Chaohui Wang

Michael M Bronstein

Nikos Paragios

Discrete minimum distortion correspondence problems for non-rigid shape matching

Keywords: shape recognition, shape similarity, Gromov-Hausdorff metric reconaissance de formes, similarité des formes, distance de Gromov-Hausdorff

Similarity and correspondence are two fundamental archetype problems in shape analysis, encountered in numerous application in computer vision and pattern recognition. Many methods for shape similarity and correspondence boil down to the minimum-distortion correspondence problem, in which two shapes are endowed with certain structure, and one attempts to find the matching with smallest structure distortion between them. Defining structures invariant to some class of shape transformations results in an invariant minimum-distortion correspondence or similarity. In this paper, we model shapes using local and global structures and formulate the invariant correspondence problem as binary graph labeling. We perform challenging non-rigid shape matching experiments, and show how different choice of structure results in invariance under various classes of deformations.

Introduction

Recent works in computer vision and shape analysis [START_REF] Berg | Shape matching and object recognition using low distortion correspondences[END_REF][START_REF] Mémoli | A theoretical and computational framework for isometry invariant recognition of point cloud data[END_REF][START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF][START_REF] Bronstein | Calculus of non-rigid surfaces for geometry and texture manipulation[END_REF] have shown that different approaches to shape similarity and correspondence can be considered as instances of the minimum distortion correspondence problem, in which two shapes are endowed with certain structure, and one attempts to find the best (least distorting) matching between these structures. Examples of such structures include multiscale heat kernel signatures [START_REF] Sun | A concise and provably informative multi-scale signature based on heat diffusion[END_REF][START_REF] Ovsjanikov | Shape Google: a computer vision approach to invariant shape retrieval[END_REF][START_REF] Bronstein | Shape Google: geometric words and expressions for invariant shape retrieval[END_REF][START_REF] Bronstein | Scale-invariant heat kernel signatures for shape recognition[END_REF], local photometric properties [START_REF] Zaharescu | Surface Feature Detection and Description with Applications to Mesh Matching[END_REF][START_REF] Thorstensen | Non-rigid Shape matching using Geometry and Photometry[END_REF], topological graphs [START_REF] Biasotti | Extended reeb graphs for surface understanding and description[END_REF], conformal factor [START_REF] Ben-Chen | Characterizing shape using conformal factors[END_REF], shape extremities [START_REF] Zhang | Feature-based surface parameterization and texture mapping[END_REF][START_REF] Zhang | Deformation-driven shape correspondence[END_REF], eigenfunctions of the Laplace-Beltrami operator [START_REF] Rustamov | Laplace-Beltrami eigenfunctions for deformation invariant shape representation[END_REF][START_REF] Mateus | Articulated shape matching using laplacian eigenfunctions and unsupervised point registration[END_REF][START_REF] Hu | Salient spectral geometric features for shape matching and retrieval[END_REF][START_REF] Dubrovina | Matching shapes by eigendecomposition of the Laplace-Beltrami operator[END_REF], triplets of points [START_REF] Lipman | Möbius voting for surface correspondence[END_REF][START_REF] Zeng | Dense non-rigid surface registration using high-order graph matching[END_REF], and geodesic [START_REF] Mémoli | A theoretical and computational framework for isometry invariant recognition of point cloud data[END_REF][START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF][START_REF] Elad | Bending invariant representations for surfaces[END_REF], diffusion [START_REF] Bronstein | A Gromov-Hausdorff framework with diffusion geometry for topologically-robust non-rigid shape matching[END_REF][START_REF] Mémoli | Spectral Gromov-Wasserstein distances for shape matching[END_REF] and commute time [START_REF] Rustamov | Laplace-Beltrami eigenfunctions for deformation invariant shape representation[END_REF][START_REF] Qiu | Clustering and embedding using commute times[END_REF][START_REF] Bronstein | On a relation between shape recognition algorithms based on distributions of distances[END_REF][START_REF] Bronstein | Shape recognition with spectral distances[END_REF] distances. By defining a structure invariant under certain class of transformations (e.g. non-rigid deformations), one obtains correspondence invariant under that class (in the above example, deformation invariant matching). The Gromov-Hausdorff distance [START_REF] Gromov | Structures Métriques Pour les Variétés Riemanniennes[END_REF] is an important particular case of the minimum distortion correspondence problem, in which the matched structures are metric spaces, invariant to isometries of the metric structures.

Some settings of the minimum distortion correspondence problem can be reformulated as labeling problems [START_REF] Torresani | Feature correspondence via graph matching: Models and global optimization[END_REF][START_REF] Starck | Correspondence labelling for widetimeframe free-form surface matching[END_REF], such that the objective function can be optimized efficiently using the recently developed discrete optimization algorithms. For example, the dual-decomposition strategy [START_REF] Bertsekas | Nonlinear Programming[END_REF], introduced by [START_REF] Komodakis | MRF optimization via dual decomposition: Message-passing revisited[END_REF] to perform pairwise Markov random field (MRF) inferences, provides a powerful technique to solve such labeling problems. Based on such a strategy, Torresani et al. proposed a pairwise graph matching algorithm [START_REF] Torresani | Feature correspondence via graph matching: Models and global optimization[END_REF] to compute correspondence between images using a criterion combining local features and Euclidean distances between nearby features. Such an approach showed better performance than feature-only based methods in deformable 2D object tracking, the increased performance attributed to the use of inter-feature distances as a geometric consistency constraint. However, Euclidean distances are not deformation invariant and can be applied only locally, thus limiting the usefulness of such a constraint.

Main contribution

In this paper, we study the minimum distortion correspondence problem in the context of non-rigid shape analysis. We formulate invariant correspondence as a minimizer of a distortion criterion based on structures invariant to some classes of transformations. In particular, we use local and global structures invariant to important classes of transformations such as non-rigid deformations, changes in topology, and scaling. By such an axiomatic construction of invariant structures, we obtain invariant correspondence. In particular, we show scale invariant shape matching using only singleton and pairwise interactions without higher-order terms. Compared to Torresani et al. [START_REF] Torresani | Feature correspondence via graph matching: Models and global optimization[END_REF], our use of global structures in non-rigid shapes provides a better regularization to the problem and is better motivated geometrically. Yet, it also increases the computational complexity of the optimization. To address this problem, we use hierarchical matching, in which candidate correspondences are restricted to neighborhoods of matching points from coarser levels.

While the described axiomatic approach is suitable for modeling geometric shape transformations such as bendings, it is not applicable to intra-class shape variations (e.g. different appearances of a human shape). To cope with this case, we show a probabilistic extension of our framework, in which local and global structures are replaced with respective multidimensional distributions, accounting for shape variability.

Related work

Feature-based shape matching methods for non-rigid shapes were used in numerous recent works [START_REF] Zaharescu | Surface Feature Detection and Description with Applications to Mesh Matching[END_REF][START_REF] Sun | A concise and provably informative multi-scale signature based on heat diffusion[END_REF][START_REF] Hu | Salient spectral geometric features for shape matching and retrieval[END_REF]. Tree-based [START_REF] Zhang | Deformation-driven shape correspondence[END_REF], RANSAC-type [START_REF] Tevs | Isometric Registration of Ambiguous and Partial Data[END_REF], and branch-and-bound techniques [START_REF] Raviv | Symmetries of non-rigid shapes[END_REF][START_REF] Raviv | Full and partial symmetries of non-rigid shapes[END_REF][START_REF] Raviv | Diffusion symmetries of non-rigid shapes[END_REF] were used to find the matches between the feature points. Elad and Kimmel [START_REF] Elad | Bending invariant representations for surfaces[END_REF] used multidimensional scaling (MDS) [START_REF] Borg | Modern multidimensional scaling -theory and applications[END_REF] to represent shapes in a low-dimensional Euclidean space and compare them as rigid objects. The use of an intermediate embedding space was eliminated in [START_REF] Mémoli | A theoretical and computational framework for isometry invariant recognition of point cloud data[END_REF] using the Gromov-Hausdorff formalism [START_REF] Gromov | Structures Métriques Pour les Variétés Riemanniennes[END_REF]. Bronstein et al. [START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF] proposed an MDS-like algorithm referred to as generalized MDS (GMDS) for the computation of the Gromov-Hausdorff distance and deformation invariant correspondence between shapes [START_REF] Bronstein | Calculus of non-rigid surfaces for geometry and texture manipulation[END_REF]. This framework was extended in [START_REF] Bronstein | A Gromov-Hausdorff framework with diffusion geometry for topologically-robust non-rigid shape matching[END_REF][START_REF] Mémoli | Spectral Gromov-Wasserstein distances for shape matching[END_REF] using diffusion geometry instead of the geodesic one. In [START_REF] Mémoli | On the use of Gromov-Hausdorff distances for shape comparison[END_REF], Mémoli extended [START_REF] Mémoli | A theoretical and computational framework for isometry invariant recognition of point cloud data[END_REF][START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF] by modeling shapes as metric-measure spaces. He introduced the Gromov-Wasserstein distance based on measure coupling between two metric-measure spaces, and formulated is as a quadratic assignment problem (QAP). Thorstensen and Keriven [START_REF] Thorstensen | Non-rigid Shape matching using Geometry and Photometry[END_REF] extended the GMDS framework to textured shapes introducing photometric stress as a local matching term in addition to geodesic distance distortion. Dubrovina and Kimmel [START_REF] Dubrovina | Matching shapes by eigendecomposition of the Laplace-Beltrami operator[END_REF] generalized this approach for the matching of textureless shapes using Laplace-Beltrami eigenfunctions as local geometric descriptors. Mateus et al. [START_REF] Mateus | Articulated shape matching using laplacian eigenfunctions and unsupervised point registration[END_REF] showed a non-rigid shape correspondence approach with inexact graph matching based on spectral embedding. In the image domain, Torresani et al. [START_REF] Torresani | Feature correspondence via graph matching: Models and global optimization[END_REF] used graph labeling problem to match 2D images.

Problem formulation

Our shape model is an extension of the metric model used in [START_REF] Mémoli | A theoretical and computational framework for isometry invariant recognition of point cloud data[END_REF][START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF][START_REF] Elad | Bending invariant representations for surfaces[END_REF]. We assume that the shapes are endowed with local and global structure, and try to find such a correspondence between the shapes that best preserves these structures. The structures are defined having in mind certain invariance properties required in the particular problem, as discussed in Section 3. Given a shape X, modeled as a connected surface (possibly with boundary) embedded into R 3 (or R 2 in case of planar shapes), its local structure is modeled by a vector field f X : X → R m referred to as a local descriptor. For example, f X can be texture [START_REF] Torresani | Feature correspondence via graph matching: Models and global optimization[END_REF], local photometric [START_REF] Zaharescu | Surface Feature Detection and Description with Applications to Mesh Matching[END_REF] or geometric [START_REF] Mitra | Probabilistic fingerprints for shapes[END_REF][START_REF] Toldo | Visual vocabulary signature for 3D object retrieval and partial matching[END_REF][START_REF] Ben-Chen | Characterizing shape using conformal factors[END_REF][START_REF] Sun | A concise and provably informative multi-scale signature based on heat diffusion[END_REF][START_REF] Ovsjanikov | Shape Google: a computer vision approach to invariant shape retrieval[END_REF][START_REF] Bronstein | Scale-invariant heat kernel signatures for shape recognition[END_REF] descriptor. The global structure of the shape is modeled as a metric d X : X × X → R, defined as a positive-definite subadditive function between pairs of points on X.

Given two shapes X and Y with the local descriptors f X and f Y and metrics d X and d Y , respectively, we define a bijective correspondence between X and Y as C ⊂ X × Y satisfying ∀x ∈ X ∃!y ∈ Y such that (x, y) ∈ C and ∀y ∈ Y ∃!x ∈ X such that (x, y) ∈ C. A good correspondence should match similar descriptors between corresponding points and similar metrics between corresponding pairs of points. This can be quantified using firstand second-order distortion terms, dis(C) = f (C) and dis(C × C) = d(C × C) , measuring the quality of correspondence of local and global structures, respectively. (here,

f (C) is a |C| × 1 vector with elements f X (x) -f Y (y) for all (x, y) ∈ C; d(C × C) is a |C| 2 × 1 vector with elements d X (x, x ) -d Y (y, y )
INRIA for all (x, y), (x , y ) ∈ C; and • is some norm). In particular,

dis 2 (C) = (x,y)∈C f X (x) -f Y (y) 2 ; dis 2 (C × C) = (x,y),(x ,y )∈C (d X (x, x ) -d Y (y, y )) 2 .
The optimal correspondence is found by minimizing a combination of first-and second-order distortion terms,

min C dis(C) + βdis(C × C), β ≥ 0. (1)
The minimizer of problem (1) is the minimum distortion correspondence between X and Y . The minimum of problem (1) can be interpreted as the similarity of X and Y . 1A particular theoretically important case is a minimum-distortion correspondence with an L ∞ second-order distortion term, referred to as the Gromov-Hausdorff distance [START_REF] Gromov | Structures Métriques Pour les Variétés Riemanniennes[END_REF]

d GH (C × C) = 1 2 min C max (x,y),(x ,y )∈C |d X (x, x ) -d Y (y, y )|.

Invariance

The choice of the local and global structures (f X , f Y and d X , d Y ) defines the invariance properties of the correspondence. Assume that the shape Y = τ (X) is obtained from X by means of some transformation τ from a class T . If

f X • τ = f Y and d X • (τ × τ ) = d Y
for all τ ∈ T , our structures are invariant under the class of transformations T . As a result, correspondence obtained by the solution of problem ( 1) is also invariant under T . Important invariance classes can be addressed by appropriate definition of the descriptors and the metric. In particular, we are interested in inelastic deformations (bendings), changing the embedding of the shape without changing its intrinsic structure; topological transformations, resulting in local changes in the connectivity of the shape, appearing as holes or "gluing" two points on the surface; and scaling.

Choice of the metric 3.2 Geodesic metric

One of the most straightforward definitions of a metric on a surface is the geodesic metric, measuring the length of a shortest path between points x and x ,

d X (x, x ) = min γ∈Γ(x,x ) (γ),
where Γ(x, x ) denotes the set of all admissible paths between x and x , γ is some admissible path, and (γ) is its length. The geodesic metric is intrinsic, dependent only on local distance structure of the shape, and is thus invariant to inelastic deformations [START_REF] Mémoli | A theoretical and computational framework for isometry invariant recognition of point cloud data[END_REF][START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF][START_REF] Elad | Bending invariant representations for surfaces[END_REF]. A notable drawback of the geodesic distance is its sensitivity to topological transformations. Connectivity changes alter the admissible paths Γ (e.g., gluing the fingers of the hand creates new paths that have not existed before), and, since the geodesic distance takes the minimum over all path lengths, sometimes the change in the geodesic metric can be very significant.

Diffusion metric

A more robust definition of an intrinsic metric based on heat diffusion properties has been recently popularized by Lafon et al. [START_REF] Coifman | Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps[END_REF]. Heat diffusion on manifolds is governed by the heat equation ∆ X + ∂ ∂t u = 0, where u is the heat distribution and ∆ X is the positive semi-definite Laplace-Beltrami operator (LBO), which can be roughly thought of as a generalization of the Laplacian to non-Euclidean domains. The heat kernel h X,t (x, z) is the solution of the heat equation with a point heat source at point x at time t = 0. For compact manifolds, the Laplace-Beltrami operator has discrete eigendecomposition of the form ∆ X φ i = λ i φ i , where λ 0 = 0, λ 1 , ... ≥ 0 are eigenvalues and φ 0 , φ 1 , ... are eigenfunctions. Using the eigenbasis of ∆ X , the heat kernel can be presented as

h X,t (x, z) = ∞ i=0 e -λit φ i (x)φ i (z). (2)
A family of metrics

d 2 X,t (x, y) = h X,t (x, •) -h X,t (y, •) 2 L2(X) = ∞ i=1 e -2λit (φ i (x) -φ i (y)) 2 ,( 3
)
parameterized by the time scale t, is referred to as diffusion metrics. Diffusion metric is inversely related to the connectivity of points x and y by paths of length t. Unlike the geodesic distance which measures the length of the shortest path, the diffusion metric has an averaging effect over all paths connecting two points. As a result, diffusion metric is less sensitive to topology and connectivity changes [START_REF] Bronstein | A Gromov-Hausdorff framework with diffusion geometry for topologically-robust non-rigid shape matching[END_REF]. With an appropriate selection of the time scale t, the effect of topological noise can be reduced [START_REF] Bronstein | On a relation between shape recognition algorithms based on distributions of distances[END_REF][START_REF] Bronstein | Shape recognition with spectral distances[END_REF].

Commute-time metric

At the same time, the need to select the scale parameter is a disadvantage, as it depends on the shape scale. Moreover, the diffusion metric is not scale invariant, since scale change affects the eigenvalues λ i and eigenfunctions φ i . A different metric,

δ 2 X (x, y) = ∞ i=1 1 λ i (φ i (x) -φ i (y)) 2 , (4)
called the commute time (or resistance [START_REF] Grinstead | Introduction to Probability[END_REF]) distance, is similar in its spirit to the diffusion metric, while being scale-invariant. The commute time metric measures the connectivity of points by paths of any length and is related to the expected time it takes a random walk initiating at point x go through point y and return to x. It is related to the diffusion distance by the following formula [START_REF] Qiu | Clustering and embedding using commute times[END_REF], Due to integration over t, δ X does not involve scale selection and is invariant to scaling transformations [START_REF] Bronstein | On a relation between shape recognition algorithms based on distributions of distances[END_REF][START_REF] Bronstein | Shape recognition with spectral distances[END_REF].

δ 2 X (x, y) = 2 ∞ 0 d 2 X,

Choice of the descriptor

Similarly to our motivation in the selection of the metric, the choice of the local descriptor is also dictated by the desired invariance properties. Due to their locality, many types of descriptors are usually less susceptible to changes as a result of nonrigid deformations. However, some descriptors have explicit invariance properties by construction.

Several works (e.g. [START_REF] Mitra | Probabilistic fingerprints for shapes[END_REF][START_REF] Toldo | Visual vocabulary signature for 3D object retrieval and partial matching[END_REF]) successfully used rigid shape patches to compare non-rigid shapes. Raviv et al. [START_REF] Raviv | Symmetries of non-rigid shapes[END_REF] used local histograms of geodesic distances as local shape descriptors. Such descriptors are deformation invariant, however, suffer from all the problems of geodesic metrics (sensitivity to topological noise and scale). Hu and Hua [START_REF] Hu | Salient spectral geometric features for shape matching and retrieval[END_REF] and Dubrovina and Kimmel [START_REF] Dubrovina | Matching shapes by eigendecomposition of the Laplace-Beltrami operator[END_REF] used the pointwise values of the Laplace-Beltrami eigenfunctions as local descriptors. Such descriptors are deformation-invariant, but sensitive to changes in topology.

Heat kernel signature

Sun et al. [START_REF] Sun | A concise and provably informative multi-scale signature based on heat diffusion[END_REF] introduced intrinsic descriptors based on multi-scale heat kernels, referred to as heat kernel signatures (HKS). The HKS is constructed at every point of the shape by considering the values of the heat kernel diagonal at multiple time scales, f X (x) = (h X,t1 (x, x), . . . , h X,tn (x, x)), where t 1 , . . . , t n are some time scale. The HKS is invariant to inelastic deformations and was also shown to be insensitive to topological transformations [START_REF] Ovsjanikov | Shape Google: a computer vision approach to invariant shape retrieval[END_REF].

Scale-invariant heat kernel signature

The disadvantage of HKS is the lack of scale invariance. In a follow-up work, Bronstein and Kokkinos [START_REF] Bronstein | Scale-invariant heat kernel signatures for shape recognition[END_REF] introduced a scale-invariant modification of HKS, referred to as SI-HKS. The main idea is to sample the time scales logarithmically (t = α τ ) such that shape scaling corresponds to a scale-space shift. Such a shift is then undone by taking the magnitude of the Fourier transform w.r.t. τ . The SI-HKS enjoys the invariance properties of HKS, while in addition also being scale-invariant.

Photometric descriptors

If in addition the shape is endowed with some photometric properties or texture, they can be used as a local descriptor. The advantage of photometric descriptor is its insensitivity to any geometric transformations. Recent works of Thorstensten and Keriven [START_REF] Thorstensen | Non-rigid Shape matching using Geometry and Photometry[END_REF] used a photometric stress in an extension of the GMDS framework [START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF]. Zaharescu et al. [START_REF] Zaharescu | Surface Feature Detection and Description with Applications to Mesh Matching[END_REF] constructed SIFT-like [START_REF] Lowe | Distinctive image features from scale-invariant keypoint[END_REF] texture descriptors on shapes. Such descriptors were shown to be very robust in shape matching with significant topological noise.

Correspondence as a graph labeling problem

Our minimum-distortion correspondence problem can be formulated as a binary labeling problem with uniqueness constraints [START_REF] Torresani | Feature correspondence via graph matching: Models and global optimization[END_REF] in a graph with vertices defined as pairs of points and edges defined as quadruplets. More formally, let V = {(x, y) :

x ∈ X, y ∈ Y } = X × Y be the set of pairs of points from X and Y , and let E = {((x, y), (x , y )) ∈ V × V and (x, y) = (x , y )}. Let L = {0, 1} further denote the set of binary labels. We can represent a correspondence C ⊆ V as binary labeling u ∈ L V of the graph (V, E), as follows: u x,y = 1 iff (x, y) ∈ C and 0 otherwise. When using L 2 distortions, the correspondence problem (1) can be reformulated as

min u∈L V (x,y)∈V u x,y ( f X (x) -f Y (y) -γ) + β ((x,y),(x ,y ))∈E u x,y u x ,y |d X (x, x ) -d Y (y, y )| 2 s.t. y u x,y ≤ 1 ∀x ∈ X; x u x,y ≤ 1 ∀y ∈ Y. ( 5
)
where γ > 0 is an occlusion term [START_REF] Torresani | Feature correspondence via graph matching: Models and global optimization[END_REF] to penalize unmatched points. We can choose a sufficiently large γ to ensure the bijective correspondence and the equivalence of the two problems 1 and 5.

In general, optimization of this energy is NP-hard [START_REF] Gold | A graduated assignment algorithm for graph matching[END_REF]. Here, we adopt the graph matching algorithm [START_REF] Torresani | Feature correspondence via graph matching: Models and global optimization[END_REF] based on dual decomposition to perform the optimization of (5). The key idea of this approach is, instead of minimizing directly the energy (5) of the original problem, to maximize a lower bound on it by solving the dual to the linear programming (LP) relaxation of ( 5). This approaches demonstrate good global convergence behavior [START_REF] Komodakis | MRF optimization via dual decomposition: Message-passing revisited[END_REF]. We first decompose the original problem, which is too complex to solve directly, into a series of sub-problems, each of which is smaller and solvable. After getting the solution of the sub-problems, we combine them using a projectedsubgradient scheme to get the solution of the original problem. In the numerical experiments, following [START_REF] Torresani | Feature correspondence via graph matching: Models and global optimization[END_REF], we decompose problem (5) into a linear subproblem, a maxflow subproblem and a set of local subproblems. of problem ( 5) is O(|V| 2 |E|) multiplied by the number of iterations, i.e., O(N 8 ). This complexity can be reduced by adopting a hierarchical matching strategy: after finding a coarse correspondence between a small number of points, correspondence between nearby points only is looked for. This allows to significantly reduce the graph size.

Hierarchical matching

Let x 1 , x 2 , . . . denote a progressive sampling of the shape X, such that X n = {x 1 , . . . , x n } constitutes an r n -covering of X (i.e., d X (X, X n ) ≤ r n , where d X is some metric on X). Such a sequence of points can be found using e.g. farthest point sampling (FPS) strategy [START_REF] Hochbaum | A best possible heuristic for the k-center problem[END_REF], in which x 1 is selected arbitrarily and the next point is selected as x k+1 = arg max x∈X min i=1,...,k d X (x, x i ). Same way, Y n = {y 1 , . . . , y n } will denote an r n -covering of Y .

At the first stage of hierarchical matching, correspondence is found between X N1 and Y N1 , where N 1 is some small number (in our experiments, it varied between 4 and 10), solving the labeling problem (5) on the full graph (

V 1 = X N1 × Y N1 , E 1 = {((x, y), (x , y )) ∈ V 1 × V 1 and (x, y) = (x , y )}. The solution provides a coarse correspondence C 1 ⊂ X N1 × Y N1 .
At the (k + 1)st level, correspondence is found between X N k+1 and Y N k+1 (the number of points is increased by a factor typically 2 ≤ q = N k+1 /N k ≤ 4), restricting the correspondence candidates for points within a certain radius around x to points within a certain radius around y, where (x, y) ∈ C k . This way, the (k + 1)st level labeling problem is solved on the graph with vertices

V k+1 = {(x i , y i ) ∈ X N k+1 × Y N k+1 : ∃(x, y) ∈ C k s.t. d X (x, x i ) < ρr k , d Y (y, y i ) < ρr k },
where ρ > 1, and

E k+1 = {((x, y), (x , y )) ∈ V k+1 × V k+1 and (x, y) = (x , y )}. For ρ ≈ 1, the size of the ρr k -neighborhood in X N k+1 of a point from X N k contains O(q) points. Thus, |V k+1 | = O(q 2 N k ), and |E k+1 | = O(q 4 N 2 k ) points, a significant reduction compared to O(N 2
k+1 ) vertices and O(N 4 k+1 ) edges in a full graph, as shown Figure 1. As a result, the complexity of the optimization becomes O(N 4 ).

Probabilistic matching and shape prototypes

While invariance to geometric transformations such as bending can be accounted by the selection of local (descriptor) and global (metric) structures, many types of shape variability cannot be accounted for in this way. For example, variability within the shape class (e.g. fat or thin man) results in different local and global structures that cannot be modeled explicitly. At the same time, such a variability can be modeled statistically. Instead of a descriptor and metric defining the shape X, we now have distributions f X ∼ F X and d X ∼ D X , e.g., Gaussian mixture model for distances,

p xx (d) = K k=1 π xx k 1 √ 2πσ xx k exp - (d -µ xx k ) 2 2σ 2 xx k , K k=1 π xx k = 1; ∀x = x ∈ X,
and descriptors,

p x (f ) = K k=1 π xk exp -1 2 (f -µ xk ) T Σ -1 xk (f -µ xk ) (2π) m/2 (det Σ xk ) 1/2 , K i=1 π xk = 1; ∀x ∈ X,
where p denotes probability density. The distance distribution between points x and x is parameterized by

D xx = {µ xx k , σ 2 xx k , π xx k } K k=1 ; the distribution of descriptors at each points x is parameterized by F x = {µ xk , Σ xk , π xk } K i=1
, where µ xk are m × 1 vectors and Σ xk are m×m matrices. We call X = ((F x ) x∈X , (D xx ) x =x ∈X ) a shape prototype.

In this probabilistic setting, given a shape Y , we determine the correspondence between Y and the prototype X by solving a problem similar to [START_REF] Borg | Modern multidimensional scaling -theory and applications[END_REF], with the distortion terms replaced by negative log-likelihood functions,

min u∈L V - (x,y)∈V u x,y (log p x (f Y (y)) + γ) - β ((x,y),(x ,y ))∈E u x,y u x ,y log p xx (d Y (y, y )) s.t. y u x,y ≤ 1 ∀x ∈ X; x u x,y ≤ 1 ∀y ∈ Y. ( 6
)

Results

To assess the performance of the presented approach, we performed multiple experiments of shape correspondence and similarity computation under a variety of transformations. Shapes from the TOSCA [START_REF] Bronstein | Numerical geometry of nonrigid shapes[END_REF] and Princeton [START_REF] Shilane | The Princeton shape benchmark[END_REF] datasets were used in our experiments. Textured shapes acquired with a multicamera system were taken from the INRIA Grenoble dataset [START_REF] Zaharescu | Surface Feature Detection and Description with Applications to Mesh Matching[END_REF]. The shapes were represented as triangular meshes with 2000-10000 vertices. Textures were given as RGB values for each vertex.Geodesic distances were computed using fast marching [START_REF] Kimmel | Computing geodesic paths on manifolds[END_REF]. Diffusion and commute time metrics were computed using the spectral formulae (3) and (4) taking the first 100 eigenvalues. The Laplace-Beltrami operator was approximated using cotangent weights [START_REF] Meyer | Discrete differentialgeometry operators for triangulated 2-manifolds[END_REF] and its eigenvalues and eigenfunctions were computed by solving the generalized eigensystem as described in [START_REF] Levy | Laplace-Beltrami eigenfunctions towards an algorithm that "understands" geometry[END_REF]. The heat kernel was approximated using formula (2). Hierarchical matching was implemented in MATLAB with discrete optimization module in C++.

INRIA

Typical running times for pairwise shape matching in the following experiments were about 10 -20 sec.

Invariance and the choice of the metric/descriptor

In the first experiment, matching was performed between eight points with equal weight given to the local and global distortion terms in the optimization problem. Three combinations of first-and second-order structures were used: geodesic metric/HKS descriptor, diffusion metric/HKS descriptor, and commute time metric/SI-HKS descriptor. Figure 2 shows the result of correspondence computation between shapes with different transformations for different choice of metric/descriptor. All three methods are invariant to bendings (first row; note that correspondence is defined up to an intrinsic symmetry). The combination geodesic metric/HKS descriptor is sensitive to topology (a human with hands glued to legs, second row) and scale. The combination diffusion metric/HKS descriptor is insensitive to topology but sensitive to scale. Finally, commute time metric with SI-HKS descriptor are invariant to all of the above.

Geodesic / HKS Diffusion / HKS Commute time / SI-HKS

Comparison to GMDS

In the second experiment, we compared our method to generalized multidimensional scaling (GMDS)2 , introduced in [8] and, as of today, a state-of-the-art method for deformable shape matching. GMDS computes correspondence between two shapes by trying to embed one shape into another with minimum distortion of geodesic distances (referred to as stress).

Given a fixed set of points x 1 , . . . , x N ∈ X, GMDS attempts to find a set of corresponding points on Y in barycentric coordinates y i = (t i , u i ) (where t i ∈ T (Y ) is a triangle index, and u ij ∈ [0, 1], j u ij = 1 is a vector of barycentric weights) by minimizing the stress

min {ti;ui} N i=1 i>j (d X (x i , x j ) -d Y ((t i ; u i ), (t j ; u j ))) 2 , (7)
where d X (x i , x j ) is a pre-computed geodesic distance between x i and x j on X, and d Y is interpolated from pre-computed geodesic distances between points on Y . While the stress function is highly non-convex, GMDS optimization is performed in a multi-resolution manner and in practice shows good convergence if initialized sufficiently close to the global minimizer [START_REF] Bronstein | Efficient computation of isometry-invariant distances between surfaces[END_REF]. A branch-and-bound technique with local histogram descriptors proposed in [START_REF] Raviv | Symmetries of non-rigid shapes[END_REF] to initialize the GMDS optimization was used in this experiment. For a fair comparison, we applied our method with geodesic distances in the second-order distortion term and no first-order term. The number of matched points was eight in both methods. Figure 3 shows the correspondence obtained using GMDS (left) and the proposed method (right). The first row shows a case when GMDS converges incorrectly. In the second row, the results of GMDS and the proposed method are similar.

Figure 4 shows similarity computed on a set of six shape classes, with four nearisometric deformations in each class (total 24 shapes) using GMDS and the proposed method. Similarity was computed as the distortion of correspondence between six corresponding points. Using similarity to rank matches in a retrieval experiment, in which the deformed shapes are used as queries, the mean average precision (mAP) achieved by our method is 100% (meaning that the first ranked matched is always from the correct shape class). For comparison, GMDS achieves only 40.1% mAP.

Matching of textured shapes

In the third experiment, to show the performance of our approach on the problem of textured shape matching from the INRIA Grenoble dataset [START_REF] Zaharescu | Surface Feature Detection and Description with Applications to Mesh Matching[END_REF]. These shapes contained various acquisition artifacts, such as topological noise and irregular triangulation. As local descriptors, we used the texture RGB values averaged in a small neighborhood of the point.

Shape prototypes

In the fourth experiment, a shape prototype was created based on 64 examples of a human shape, in which the length of the hands and legs and the size of the head was varied. Distance and descriptor distributions were represented using Gaussian mixtures with 5 components. Figure 6 shows a comparison of deterministic and probabilistic matching. Using deterministic matching, the shape of a humanoid alien from the Princeton database [START_REF] Shilane | The Princeton shape benchmark[END_REF]. is matched to the human shape from TOSCA dataset incorrectly (second column from left), because of different proportions of the head, legs, and hands. On the other hand, matching to the human shape prototype using probabilistic matching produced correct symmetric correspondence (third column). Figure 6 (columns four and five) shows additional examples of shape prototype matching. These results show that the probabilistic matching framework allows to address shape variability that cannot be simply accommodated into the metric model by choosing the metric.

Conclusions

We presented a generic framework for invariant matching between shapes, in which matching is performed by minimizing the distortion of local and global geometric structures under the correspondence. Using structures invariant to pre-defined classes of transformations (or, using their statistical distributions if such transformations cannot be modeled explicitly) allows obtaining invariant matching between shapes. Our approach generalizes many previous works in the field, in particular, methods based on metric distortion minimization [START_REF] Mémoli | A theoretical and computational framework for isometry invariant recognition of point cloud data[END_REF][START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF][START_REF] Bronstein | A Gromov-Hausdorff framework with diffusion geometry for topologically-robust non-rigid shape matching[END_REF] and global and local features [START_REF] Thorstensen | Non-rigid Shape matching using Geometry and Photometry[END_REF][START_REF] Dubrovina | Matching shapes by eigendecomposition of the Laplace-Beltrami operator[END_REF][START_REF] Torresani | Feature correspondence via graph matching: Models and global optimization[END_REF], allowing incorporating many existing geometries and local descriptors [START_REF] Sun | A concise and provably informative multi-scale signature based on heat diffusion[END_REF][START_REF] Bronstein | Scale-invariant heat kernel signatures for shape recognition[END_REF][START_REF] Zaharescu | Surface Feature Detection and Description with Applications to Mesh Matching[END_REF][START_REF] Ben-Chen | Characterizing shape using conformal factors[END_REF].

In particular, it extends the Gromov-Hausdorff framework [START_REF] Mémoli | A theoretical and computational framework for isometry invariant recognition of point cloud data[END_REF][START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF][START_REF] Gromov | Structures Métriques Pour les Variétés Riemanniennes[END_REF]. Formulating the problem as graph labeling, we use powerful optimization method recently developed for this class of problems which are known to have favorable convergence properties. Our approach is especially appropriate for the challenging problems of finding similarity and correspondence between non-rigid shapes.

Limitations and extensions

The problem of symmetric correspondences, inherent to all approaches based on intrinsic structures, cannot be resolved without resorting to some side information. There are a few potential cures to this problem. First, providing some initial correspondence between the shapes could be used to restrict the vertex set, ruling out symmetric correspondences. Second, exploiting shape orientation could be used to find orientationconsistent matches. Finally, using higher-order distortions (in particular, third-order between triplets of points) can resolve the symmetry problem [START_REF] Zeng | Dense non-rigid surface registration using high-order graph matching[END_REF].

INRIA

Though the presented approach shows better global convergence compared to GMDS, the advantage of GMDS is that it works with barycentric coordinates, thus finding correspondences at the sub-vertex resolution. This effect is especially significant in coarse meshes, where vertex-to-vertex correspondence can be too coarse. The proposed algorithm, can be used as an initialization to GMDS (which requires good initialization), e.g., replacing the branch-and-bound approach used in [START_REF] Raviv | Symmetries of non-rigid shapes[END_REF]. This way, GMDS would acts as sub-vertex "refinement" stage.

Figure 1 :

 1 Figure 1: Toy example of hierarchical matching. Left: at the first level, N1 = 3, V1 = X1 ×Y1 (all possible pairs of points). The graph size is |V1| = 9, |E1| = 36, counting only (N 2 1 -N1)/2 distances. Right: at the second level, N2 = 6 (3 points are added; q = 2). The graph size is |V2| = 12, |E2| = 66 (compare to full graph with |V2| = 36, |E2| = 630).

Figure 2 :

 2 Figure 2: Invariance to different types of transformations and the choice of the metric/descriptor. Shown is matching between isometric deformations (first row), shapes with different topology (second row), and shapes with different scale (third row), using geodesic metric and HKS descriptors (left), diffusion metric and HKS descriptors (middle), and commute time metric and scale-invariant HKS descriptors (right).

Figure 3 :

 3 Figure 3: Comparison to GMDS.

Figure 4 :

 4 Figure 4: Dissimilarity between deformations of six classes of shapes computed using GMDS (left) the proposed method (right). Colormap encodes small dissimilarity values by darker colors. Ideal dissimilarity matrix has six 4 × 4 blocks of zeros on the diagonal and large values outside the diagonal.

Figure 5 :

 5 Figure 5: Matching of textured surfaces using commute time metric and RGB texture.

Figure 6 :

 6 Figure 6: Matching of an alien shape to the human shape (first column from left) using deterministic (second column) and probabilistic (third column) approaches. Columns four and five: additional probabilistic matching examples.

1 1 . 2 1 3 . 2 3 . 3 8 4 4 . 1 8 5 1

 1121323384181 Main contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Choice of the metric . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Geodesic metric . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Diffusion metric . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3.4 Commute-time metric . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3.5 Choice of the descriptor . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.6 Heat kernel signature . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.7 Scale-invariant heat kernel signature . . . . . . . . . . . . . . . . . . 7 3.8 Photometric descriptors . . . . . . . . . . . . . . . . . . . . . . . . . Correspondence as a graph labeling problem 8 Hierarchical matching . . . . . . . . . . . . . . . . . . . . . . . . . . Probabilistic Invariance and the choice of the metric/descriptor . . . . . . . . . . . 11 6.2 Comparison to GMDS . . . . . . . . . . . . . . . . . . . . . . . . . 12 6.3 Matching of textured shapes . . . . . . . . . . . . . . . . . . . . . . 13 6.4 Shape prototypes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 7 Conclusions 14 7.1 Limitations and extensions . . . . . . . . . . . . . . . . . . . . . . . 14 RR n°7333

Table 1 :

 1 Invariance properties of local (top rows) and global (bottom rows) structures.

	INRIA

t (x, y)dt. * Assuming uniform illumination. * * Only point-wise, assuming no averaging is done on texture.

RR n°7333

The minimizer of problem (1) is not necessarily unique, i.e., there may be two different correspondences C = C with the same distortion. Such situations are typical when the shapes have intrinsic symmetries. Intrinsic symmetry is manifested by the existence of a self-isometry of X with respect to the metric d X , i.e., an automorphism g : X → X satisfying d X = d X • (g × g)[START_REF] Raviv | Symmetries of non-rigid shapes[END_REF][START_REF] Ovsjanikov | Global intrinsic symmetries of shapes[END_REF][START_REF] Raviv | Full and partial symmetries of non-rigid shapes[END_REF][START_REF] Raviv | Diffusion symmetries of non-rigid shapes[END_REF]. If g and h are intrinsic symmetries of X and Y , respectively, the distortion of the correspondences C and C • (g × h) is identical and the minimum-distortion correspondence problem has multiple minimizers[START_REF] Bronstein | Calculus of non-rigid surfaces for geometry and texture manipulation[END_REF][START_REF] Raviv | Symmetries of non-rigid shapes[END_REF].RR n°7333

Code taken from tosca.cs.technion.ac.il. INRIA


Julien Rabin
email: rabin@ceremade.dauphine.fr

Gabriel Peyré
email: peyre@ceremade.dauphine.fr

Laurent D Cohen
email: cohen@ceremade.dauphine.fr

Geodesic Shape Retrieval via Optimal Mass Transport ⋆

This paper presents a new method for 2-D and 3-D shape retrieval based on geodesic signatures. These signatures are high dimensional statistical distributions computed by extracting several features from the set of geodesic distance maps to each point. The resulting high dimensional distributions are matched to perform retrieval using a fast approximate Wasserstein metric. This allows to propose a unifying framework for the compact description of planar shapes and 3-D surfaces.

Introduction

Content based 2-D and 3-D shape retrieval is an important problem in computer vision. It requires to design both representations and similarity measures to discriminate shapes from different classes, while being invariant to some deformations.

Feature-based Shape Retrieval

There is a large amount of literature on content-based retrieval using similarity measures between descriptors. In this section, a brief review is given, focusing on bending and isometric deformations (i.e. preserving the topology). We refer the reader to the following review papers devoted to planar shapes [START_REF] Veltkamp | Properties and performance of shape similarity measures[END_REF][START_REF] Zhang | Review of shape representation and description techniques[END_REF] and 3-D surfaces [START_REF] Bustos | Feature-based similarity search in 3D object databases[END_REF][START_REF] Tangelder | A survey of content based 3D shape retrieval methods[END_REF] retrieval for a complete review.

Global descriptors. Simple global features are computed using polynomial moments [START_REF] Teague | Image analysis via the general theory of moments[END_REF][START_REF] Teh | On image analysis by the methods of moments[END_REF][START_REF] Liao | On image analysis by moments[END_REF], or Fourier transform [START_REF] Zahn | Fourier descriptors for plane closed curves[END_REF] (see [START_REF] Prokop | A survey of moment-based techniques for unoccluded object representation and recognition[END_REF] for review).

The spectrum of the Laplace Beltrami operator defines a descriptor invariant to rigid motion and to simple bendings [START_REF] Reuter | Laplace-spectra as fingerprints for shape matching[END_REF]. Shape distributions [START_REF] Osada | Shape distributions[END_REF] compute descriptors as histograms of the distribution of Euclidean distance between points on the surface. This is extended to bending invariant descriptors in [START_REF] Ben Hamza | Geodesic matching of triangulated surfaces[END_REF][START_REF] Ion | Shape matching using the geodesic eccentricity transform -a study[END_REF][START_REF] Ion | 3D shape matching by geodesic eccentricity[END_REF] using geodesic distances. It is possible to replace the geodesic distance by a diffusion distance [START_REF] Gorelick | Shape representation and classification using the Poisson equation[END_REF] computed by solving a linear Poisson PDE.

Local descriptors. Many other shape representations do not make use of a single descriptor. They rather compute similarities by matching points of interest for which local descriptors are defined. Shape context features [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF] are local 2-D histograms of contours around points of interest. Geodesic shape context makes use of geodesic curves to gain bending invariance [START_REF] Ling | Shape classification using the inner-distance[END_REF]. Local tomographic projections on tangent planes (spin images) [START_REF] Johnson | Using spin images for efficient object recognition in cluttered 3D scenes[END_REF] define a set of local descriptors.

Similarity measure. Most of the previous approaches make use of Euclidean metric, Kullback-Leibler or χ 2 distance to compare low-dimensional histogrambased descriptors in linear time. When considering high-dimensional descriptors (either histograms or discrete point clouds), another possibility is to use the Wasserstein distance [START_REF] Villani | Topics in Optimal Transportation[END_REF], see e.g. [START_REF] Rubner | The Earth Mover's Distance as a Metric for Image Retrieval[END_REF][START_REF] Grauman | Fast contour matching using approximate earth mover's distance[END_REF][START_REF] Ling | An Efficient Earth Mover's Distance Algorithm for Robust Histogram Comparison[END_REF].

Contributions

This paper introduces a novel framework for bending invariant recognition of shapes. We use the setting of geodesic distances on Riemannian manifolds, which unifies both planar shape and 3-D surface retrieval problems. This novel framework builds on several already known statistical descriptors, and encompasses them into a single high-dimensional descriptor. This allows us to take advantage of the richness of information available in each separate statistical measure to enhance the retrieval performance. The retrieval method is based on an approximation of the Wasserstein distance, that works directly over discrete point clouds, and can be computed with an iterative algorithm.

Geodesic Distances

In the following, we consider shapes as compact 2-D manifolds Ω ⊂ IR s , where s = 2 (planar shapes) or s = 3 (surfaces). Note however that our approach is generic and accommodates for domains of arbitrary dimension.

Geodesic Distance Definition

The length of a curve γ : [0, 1] → Ω traced within the domain is defined as

L(γ) = 1 0 ||γ ′ (t)||dt.
The geodesic distance between two points x s , x e ∈ Ω is the length of the shortest piecewise smooth curve joining the two points

d Ω (x s , x e ) = min γ(0)=xs,γ(1)=xe L(γ). (1)
The geodesic map d Ω (x i , x) differs significantly from the the Euclidean distance map ||x i -x|| when the shapes are non convex, as it is illustrated by Fig. 1.

A curve γ ⋆ satisfying d Ω (x s , x e ) = L(γ ⋆ ) is called a shortest path, sometimes also referred to as a (globally minimizing) geodesic.

Geodesic Distance Computation

Geodesic distance within a planar shape. Given some starting point x s , the geodesic distance map U xs (x) = d Ω (x s , x) can be shown to be the unique viscosity solution of the following non-linear PDE,

∀ x ∈ Ω, ||∇U xs (x)|| = 1 and U xs (x s ) = 0. ( 2
)
where the derivative should be understood in a weak sense at points along the medial axis of x s where U xs is not smooth.

The PDE (2) can be discretized with upwind finite difference. The resulting discrete equation can be and solved in O N log(N ) operations using the Fast Marching algorithm [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF][START_REF] Tsitsiklis | Efficient Algorithms for Globally Optimal Trajectories[END_REF]. This algorithm performs a front propagation within the shape, as displayed on Fig. 2. Geodesic distance on a 3-D surface. If the surface Ω is parametrized using ψ : V ⊂ [0, 1] 2 → Ω, then one can prove that the distance map

∀ x ∈ V, U xs (x) = d Ω (x s , ψ(x)) satisfies an anisotropic Eikonal PDE ∀ x ∈ V, ||∇U xs (x)|| T -1 x = 1, and U xs (ψ(x s )) = 0, (3)
where

T x = ∂ψ ∂x i , ∂ψ ∂x j 0 i,j 1 and ||x|| A = 0 i,j 1 A i,j x i x j .
This equation (3) extends to surfaces of arbitrary topology using several charts that parametrize locally the surface. The Eikonal equation ( 3) can be discretized on 3-D meshes. In the case of mesh with no obtuse angle, the discrete equation can be solved in O N log(N ) operations [START_REF] Kimmel | Computing Geodesic Paths on Manifolds[END_REF]. For general meshes, the resolution requires more advanced schemes, see for instance [START_REF] Bornemann | Finite-element discretization of static Hamilton-Jacobi equations based on a local variational principle[END_REF]. An example of a Fast Marching propagation from a set of starting points on a 3-D shape is given on Fig. 3.

Geodesic Descriptors

Similarity measures between shapes are computed by extracting global or local features ϕ(Ω), and then performing some comparison between the resulting descriptors.

An important goal in designing a similarity measure is to achieve invariance to some class R of deformations. This requires that the descriptors are invariant, so that ϕ(RΩ) = ϕ(Ω) for any R ∈ R.

This section details a class of geodesic descriptors that are invariant under geodesic isometries, and quasi-invariant to shape articulations and bendings. This is especially relevant to perform robust retrieval on articulated shapes, such as animal or human with varying poses.

Local Descriptors

Geodesic distance distributions. To design features invariant to bendings and articulations, we consider, for each point x ∈ S ⊂ Ω, the set {d Ω (x, y)} y∈E⊂Ω ⊂ IR + of distances to a subset E ⊂ Ω. This set of distances should be thought as being a 1-D distribution of values in IR + .

For numerical applications, the set S is a discrete sub-sampling of the manifold computed as described in Sect. 3.2. The set E used to compute the distributions can be defined depending on the application. In our numerical examples, we choose E = ∂Ω to be the boundary of the manifold for 2-D shapes, and E = Ω for 3-D surfaces.

Figure 4 shows examples of geodesic distance distribution, conveniently displayed using 1-D histograms. Geodesic quantile measures. The whole set of distances {d Ω (x, y), y ∈ E, x ∈ S} is too large to be used for retrieval applications. To achieve dimensionality reduction, we retain only a few statistical measures out of this distribution of values. This article considers quantiles statistical measures Q x (α) defined as, for all α ∈ [0, 1],

∀x ∈ S, Q x (α) = F -1 x (α) = max{δ ∈ IR + , F x (δ) α} (4)
where F x is the cumulative distribution function of the set {d Ω (x, y), y ∈ E}, and F -1

x is its pseudo-inverse. Observe that Q x (0) is the minimum geodesic distance between x and E, while Q x (1/2) is the median distance. The maximum distance Q x (1) is also known as the eccentricity, and has been used for 2-D shapes [START_REF] Ion | Shape matching using the geodesic eccentricity transform -a study[END_REF] and surfaces [START_REF] Ion | 3D shape matching by geodesic eccentricity[END_REF] retrieval. Other statistical measures can be retained as well. For instance, the mean of the distance E d Ω (x, y)dy is used in [START_REF] Ben Hamza | Geodesic matching of triangulated surfaces[END_REF] to perform surface retrieval.

Geodesic local descriptors. At each location x ∈ S, the local descriptor p x ∈ IR d is a vector of d quantiles

p x = (Q x (α ℓ )) 1 ℓ d ∈ IR d ,
where 0 α ℓ 1 are equi-spaced values. Figure 5 displays each of the d = 3 components Q x (α ℓ ) of p x as a function of x ∈ Ω, with α ℓ ∈ {0, 1 2 , 1}. The key feature of our approach, that makes it significantly different from these previous works is that it uses several statistical measures, and thus builds high dimensional descriptors.

Global Descriptors

The local descriptors p x are sampled on a set S ⊂ Ω to obtain a global descriptor that characterizes the shape.

Farthest Point Sampling. Estimating the full set {p x } x∈Ω of descriptors is computationally intractable, and one thus needs to compute a sub-sampling S = {x i } i∈I of n points on the manifold, where I = {0, . . . , n -1}. To perform a uniform sampling of the manifold, we use the farthest point sampling strategy. It corresponds to a greedy scheme, originally introduced in [START_REF] Gonzalez | Clustering to minimize the maximum intercluster distance[END_REF], and extended to geodesic distances on manifolds for surface remeshing in [START_REF] Peyré | Geodesic remeshing using front propagation[END_REF].

The initial point x 0 ∈ Ω is sampled at random. Given a set of k points {x 0 , . . . , x k-1 }, the next point is computed as

x k = argmax x∈Ω min 0 i<k d Ω (x i , x).
Once this new point x k has been computed, the set {d Ω (x k , x)} x∈E of geodesic distances is computed in O n log(n) operations, and the geodesic descriptor p x k ∈ IR d is obtained by computing the quantiles (4) from these distances. This process of iteratively adding the furthest point to the set S is continued until a given number of points n is reached. Examples of this farthest point sampling method on a 2-D shape and a surface are shown in Fig. 6.

Global descriptor as a point cloud. The global descriptor is then defined as a uniform sampling of the local descriptors

ϕ(Ω) = {p xi } i∈I ⊂ IR d ,
and is thus a cloud of n points in IR d . This point cloud ϕ(Ω) should be thought as being drawn from a probability distribution. Each shape has its own distribution, that is invariant under isometric bending of the shape. An alternative representation is to compute the d-dimensional histogram of the distribution. We prefer in this paper to work directly using discretized point cloud, because it offers a more precise matching.

The resulting global descriptor ϕ(Ω) is invariant to isometric deformations, in the sense that ϕ(Ω) = ϕ(RΩ) if R is a deformation of Ω that maintains geodesic distances. More generally, if R does not modify too much the distances, meaning

∀ (x, y) ∈ Ω 2 , d Ω (x, y) ≈ d RΩ (Rx, Ry), (5)
then ϕ(Ω) ≈ ϕ(RΩ). This is the case for bending deformations and articulations, see [START_REF] Ling | Shape classification using the inner-distance[END_REF]. Observe that geodesic distances {d Ω (x, y)} x,y ∈Ω have to be normalized (according to maximum distance) to achieve invariance to scaling.

Optimal Transport Retrieval

Our shape retrieval method uses a similarity measure that compares the geodesic descriptors ϕ(Ω) using the Wasserstein metric related to the Monge-Kantorovich optimal transport problem (see [START_REF] Villani | Topics in Optimal Transportation[END_REF] for an in-depth study).

Similarity Measure

In classical settings, shapes are generally represented by histogram-based descriptors and thus compared with L p distances or the Kullback-Leibler divergence. In our setting, the descriptors are high dimensional discrete distributions, and we propose to use the Wasserstein distance [START_REF] Villani | Topics in Optimal Transportation[END_REF] which is well adapted to compare statistical distributions [START_REF] Rubner | The Earth Mover's Distance as a Metric for Image Retrieval[END_REF], and is know to be more robust that Hausdorff distance [START_REF] Veltkamp | State of the art in shape matching[END_REF]. Our similarity measure is thus defined as

∆(Ω 0 , Ω 1 ) = W (ϕ(Ω 0 ), ϕ(Ω 1 )),
where W (X, Y ) is the Wasserstein distance between two point clouds X, Y ⊂ IR d , defined in the next section. Since our geodesic descriptors satisfy the approximate invariance (5) for bendings and articulations R ∈ R, our similarity measure satisfies ∆(Ω 0 , R(Ω 0 )) ≈ 0.

Wasserstein Distance

Given two point clouds X, Y ⊂ IR d of n elements, the L 2 -Wasserstein distance is defined as

W (X, Y ) 2 = min σ∈Σn i∈I ||X i -Y σ(i) || 2 , (6)
where Σ n is the set of all permutations of n elements and I = {0, . . . , n -1}.

Our framework extends to arbitrary strictly convex cost such as L p -metric for p > 1. Computing this distance boils down to estimate the optimal assignment i → σ ⋆ (i) minimizing Formula ( 6). This can be computed exactly using linear programming or other dedicated algorithms in O n 2.5 log(n) operations [START_REF] Burkard | Assignment Problems[END_REF].

One-dimensional case. It is well known that the Wasserstein assignment problem in 1-D can be easily solved in O n log(n) operations by sorting the values [START_REF] Villani | Topics in Optimal Transportation[END_REF]. Indeed, σ X and σ Y being two permutations such that {X σ X (i) } i and {Y σ Y (i) } i are sorted in increasing order, the optimal assignment is

σ ⋆ = σ Y • σ -1 X . (7)

Approximate Wasserstein Distance

For large point clouds in high dimension (d 2), computing exactly ( 6) is too demanding. Following an idea recently introduced in [START_REF] Rabin | Wasserstein Barycenter and its Application to Texture Mixing[END_REF], we propose to use the an approximate transport cost W (X, Y ) defined as

W (X, Y ) = || X -X (∞) || , (8)
where X (∞) is computed using an iterative algorithm described in the following paragraph. Starting from X (0) = X, this algorithm computes points clouds {X (k) } k that progressively evolves X (k) toward Y , minimizing at each iteration an energy E Y which is a sum of 1-D Wasserstein distances on the unit sphere

S d-1 in IR d : E Y (U ) = 1 2 θ∈S d i∈I U i -Y σ ⋆ θ (i) , θ 2 dθ , (9)
where . , . is the L 2 -scalar product, and where σ ⋆ θ is the optimal 1-D assignment according to direction θ of { Y i , θ } i with { U i , θ } i following Formula [START_REF] Liao | On image analysis by moments[END_REF].

Algorithm. Finding a minimum of energy ( 9) can be done using a classical gradient descent strategy. For numerical considerations, this energy is estimated at each iteration k from a restricted set of directions, thus resulting in a stochastic gradient descent scheme (see e.g. [START_REF] Bottou | Online algorithms and stochastic approximations[END_REF]), which relies on three steps: Y is therefore

⊲ Step 1 Define the direction set Ψ (k) ⊂ S d-1 ,
E (k) Y U = 1 2 θ∈Ψ (k) i∈I U i -Y σ ⋆ θ (i) , θ 2 . ( 10
) ⊲ Step 2 Compute, for each direction θ ∈ Ψ (k) , the optimal 1-D assignment σ ⋆ θ of 1-D distribution { Y i , θ } i∈I with { X (k)
i , θ } i∈I using Formula ( 7). ⊲ Step 3 The set {σ ⋆ θ } being computed, update the point cloud using a Newton step with parameter λ ∈ ]0, 2[ to minimize the energy

E (k) Y , such that ∀ i ∈ I X (k+1) i = X (k) i -λ ∇ 2 E (k) Y X (k) -1 ∇E (k) Y X (k) i , = (1 -λ) • X (k) i + λ H -1 θ∈Ψ (k) Y σ ⋆ θ (i) , θ θ , (11)
where

H = θ∈Ψ (k) θθ T ∈ IR d×d is the Hessian matrix of E (k) Y at point X (k) i .
Convergence. Stochastic gradient descent is known to converge if one uses a properly chosen step size λ = λ k that decays through the iterations, see [START_REF] Bottou | Online algorithms and stochastic approximations[END_REF]. In numerical simulations, we always observed convergence of X (k) to some X (∞) using the fixed step size λ = 1/|Ψ (k) |. Furthermore, X (∞) is actually equal (up to a permutation) to Y , so that the algorithm computes an assignment between the two distributions -which is not necessarily optimal in the sense of Formula (6).

Implementation. In the numerical simulations of Sect. 5, we used a fixed number of K = 100 iterations and |Ψ | = d directions, and we noticed that using more iterations does not improve significantly the retrieval results. The complexity of the proposed algorithm is therefore O |Ψ |Kn log(n) . This method extends the algorithm proposed by [START_REF] Pitié | Automated colour grading using colour distribution transfer[END_REF] that makes use of an orthogonal set of direction Ψ k , and a descent step size λ = 1. Using a smaller step, e.g. λ = 1/|Ψ k | is important to ensure the convergence of the method. Using a larger set of directions is useful to obtain a assignment that is closer to the optimal one, see [START_REF] Rabin | Wasserstein Barycenter and its Application to Texture Mixing[END_REF] for more details. Observe that other approximation methods has been previously proposed in the literature, e.g. making use of metric embedding [START_REF] Grauman | Fast contour matching using approximate earth mover's distance[END_REF] or wavelet approximations [START_REF] Shirdhonkar | Approximate Earth Mover's Distance in linear time[END_REF].

Numerical Examples

Given a database of manifolds {Ω j } j∈J , our algorithm for shape retrieval begins by computing the global signature ϕ(Ω j ) for each manifold in the database. This is performed in parallel to the farthest point sampling algorithm described in Sect. 3.2. When an input manifold Ω is queried in the database, its global signature ϕ(Ω) is computed, and the shape in the database are ordered according to the Wasserstein distance approximation W (ϕ(Ω), ϕ(Ω j )).

To evaluate the retrieval performance of the proposed descriptor, two classical performances curves are displayed:

the average recall curve shows the average number of correct shapes (or "true-positives") retrieved per query among the r most similar ones in the dataset. This curve, plotted depending of the rank r, is obtained by querying every shape in the dataset (the query itself is no used to compute the score);

the average precision-recall curve plots the rate of correct shapes retrieved among r as a function of the average recall rate.

In order to show the interest of considering high dimensional geodesic statistics, each performance curves are shown for two different descriptors: the aforementioned multi-dimensional descriptor and also a simple 1-D descriptor corresponding to the distribution of eccentricity (maximal geodesic distances).

2-D Shape Retrieval

In this setting, the domain E of ending points {y j } j∈J is the boundary ∂Ω of the manifold. 4-dimensional distributions of n = 500 points are used as global descriptors, considering 3 quantiles Q x (α j ) (minimum, median and maximum) in addition with the mean values of geodesic distances.

We consider first the "Articulated Shapes" dataset of [START_REF] Ling | Shape classification using the inner-distance[END_REF] (see Fig. 7(a)), a small dataset being designed to evaluate the robustness of retrieval to bending deformations. Performance curves are shown in Figs. 7(b) and 7(c) for both 4-D and 1-D descriptors. Comparison with state-of-the-art methods [START_REF] Ling | Shape classification using the inner-distance[END_REF][START_REF] Ling | An Efficient Earth Mover's Distance Algorithm for Robust Histogram Comparison[END_REF] is also provided in Table 1. One can see that considering several geodesic statistics at the same time enables to catch more sophisticated information on the shape while being more robust to bending deformations. Note that it is not the case when using Euclidean metric with descriptors made of 4 1-D histograms instead of the approximate Wasserstein metric with 4-D discrete distributions.

Table 1. Retrieval results on the articulated shapes dataset [START_REF] Ling | Shape classification using the inner-distance[END_REF]. Scores correspond to the number of correct shapes retrieved among 40 depending on their rank. In order to evaluate the robustness of our approach for 2-D shape retrieval on a larger database, we consider now the MPEG-7 dataset of 1400 shapes (see Fig. 8(a)). Results are shown in Figs. 8(b) and 8(c). Again, one can see that using multi-dimensional statistics on geodesic distances yields far better results. The state-of-the art method of [START_REF] Ling | An Efficient Earth Mover's Distance Algorithm for Robust Histogram Comparison[END_REF] yields a bullseye score (average recall rate for rank r = 40) of 86.56%, whereas we obtain 59.7%. This can be explained by the strongly non-isometric variations of the objects in the MPEG-7 classes, which makes our representation quite inefficient for this retrieval. An important avenue for future work is to design a large benchmark of planar shapes undergoing bendings and articulation deformations, to explore the performance of our algorithm and related methods.

Method

3-D Shape Retrieval

For surface, the domain E = Ω is the whole manifold. Hence, the first order quantiles Q x (0) are discarded since they are zero, so that we handle d = 3 -dimensional geodesic statistics of n = 500 points, considering from now max and median in addition with mean values of geodesic distances.

To evaluate the robustness of such global descriptor for 3-D shapes, we used the McGill dataset of 3-D articulated objects [START_REF]McGill 3-D shapes dataset[END_REF] (see Fig. 9(a) for an overview). Retrieval results are shown in Figs. 9(b) and 9(c). Again, it is clear that the combination of several geodesic distance characteristics achieves better retrieval results than considering only one. A comparison with state-of-the-art approach described in [START_REF] Siddiqi | Retrieving articulated 3-D models using medial surfaces[END_REF] is also provided in Figure 9, where we obtain similar results. Following the same protocol as in [START_REF] Siddiqi | Retrieving articulated 3-D models using medial surfaces[END_REF], retrieval results are given for a subset of the complete McGill Database (14 classes out of 19).

Conclusion

The first contribution of this paper is a generic framework to represent manifolds with statistical signatures based on geodesic distances, which are robust to bendings. The second contribution of this paper is an algorithm to compute a similarity measure between multi-dimensional joint distributions, which yields a fast approximation of the Wasserstein metric. This algorithm is applied to perform shape retrieval using our geodesic framework.

Our framework extends naturally to include additional information such as texture. One can indeed use a non-constant Riemannian metric that takes into account this information. [START_REF] Siddiqi | Retrieving articulated 3-D models using medial surfaces[END_REF]) and 9(g)): comparison of our approach with stateof-the-art method described in [START_REF] Siddiqi | Retrieving articulated 3-D models using medial surfaces[END_REF]

Figure 1 (Fig. 1 .

 11 Fig. 1. Left and center: comparison of Euclidean and geodesic distances inside a 2-D shape. Right: display of geodesic curves.

Fig. 2 .

 2 Fig. 2. Fast Marching propagation inside a 2-D shape.

Fig. 3 .

 3 Fig. 3. Example of Fast Marching propagation on a triangulated mesh.

Fig. 4 .

 4 Fig. 4. Histogram of the distribution of the geodesic distance to several points.

Fig. 5 .

 5 Fig. 5. Display of x → Qx(α ℓ ) for several α ℓ ∈ {0, 1 2 , 1} and of the corresponding 3-D distribution {px = (Qx(α ℓ )) 1 ℓ 3 }x∈S ⊂ IR 3 .

Fig. 6 .

 6 Fig. 6. Illustration of the farthest point sampling strategy for a 2-D shape and a 3-D surface. For each case |S| = n = 100 points are sampled (red dots). The geodesic distances between each x ∈ Ω and these points are plotted as a colormap.

  a collection of vectors randomly and uniformly sampled on S d-1 . The corresponding energy E (k)

( a )

 a Articulated shapes dataset of[START_REF] Ling | Shape classification using the inner-distance[END_REF]. Pairs of shapes from different classes. The complete dataset is composed of 8 classes of 5 elements.

  Recall vs Image Rank.

  Average Precision-Recall.

Fig. 7 .

 7 Fig. 7. Retrieval results on the articulated shapes database [17].Figure 7(a): database overview.Figure 7(b): Average recall rate depending on the shape rank threshold for each query shape in the dataset (the score does not include the query itself).Figure 7(c): Average Precision versus Recall.

  Fig. 7. Retrieval results on the articulated shapes database [17].Figure 7(a): database overview.Figure 7(b): Average recall rate depending on the shape rank threshold for each query shape in the dataset (the score does not include the query itself).Figure 7(c): Average Precision versus Recall.

  Fig. 7. Retrieval results on the articulated shapes database [17].Figure 7(a): database overview.Figure 7(b): Average recall rate depending on the shape rank threshold for each query shape in the dataset (the score does not include the query itself).Figure 7(c): Average Precision versus Recall.

Figure 7 (

 7 Fig. 7. Retrieval results on the articulated shapes database [17].Figure 7(a): database overview.Figure 7(b): Average recall rate depending on the shape rank threshold for each query shape in the dataset (the score does not include the query itself).Figure 7(c): Average Precision versus Recall.

Fig. 8 .

 8 Fig. 8. Retrieval results on MPEG7 database.Figure 8(a): database overview.Figure 8(b): Average recall rate. Figure 8(c): Average Precision versus Recall.

Figure 8 (

 8 Fig. 8. Retrieval results on MPEG7 database.Figure 8(a): database overview.Figure 8(b): Average recall rate. Figure 8(c): Average Precision versus Recall.

Figure 8 (

 8 Fig. 8. Retrieval results on MPEG7 database.Figure 8(a): database overview.Figure 8(b): Average recall rate. Figure 8(c): Average Precision versus Recall.

( a )

 a McGill dataset of articulated and non-articulated objects[START_REF]McGill 3-D shapes dataset[END_REF] (respectively composed of 9 classes of 202 elements and 10 classes of 255 elements).

  Recall vs Image Rank on articulated dataset.

  Recall vs Image Rank on complete dataset.

  Average Precision-Recall on complete dataset. (f) Average Recall on partial McGill database (14 classes) with [36] method.

  partial datase (g) Average Recall with our approach.

Fig. 9 .

 9 Fig. 9. Retrieval results on McGill database [35].Figure 9(a): overview of the database. Figure 9(b) and 9(d): Average recall rate on articulated and complete dataset.Figure 9(c) and 9(e): Average Precision versus Recall on articulated and complete dataset. Figures 9(f) (from[START_REF] Siddiqi | Retrieving articulated 3-D models using medial surfaces[END_REF]) and 9(g)): comparison of our approach with stateof-the-art method described in[START_REF] Siddiqi | Retrieving articulated 3-D models using medial surfaces[END_REF] on partial McGill Dataset (14 classes out of 19). Results are shown as average recall rate curve (plotted in black) along with the average intra-class standard deviation (in red).

Figure 9 (

 9 Fig. 9. Retrieval results on McGill database [35].Figure 9(a): overview of the database. Figure 9(b) and 9(d): Average recall rate on articulated and complete dataset.Figure 9(c) and 9(e): Average Precision versus Recall on articulated and complete dataset. Figures 9(f) (from[START_REF] Siddiqi | Retrieving articulated 3-D models using medial surfaces[END_REF]) and 9(g)): comparison of our approach with stateof-the-art method described in[START_REF] Siddiqi | Retrieving articulated 3-D models using medial surfaces[END_REF] on partial McGill Dataset (14 classes out of 19). Results are shown as average recall rate curve (plotted in black) along with the average intra-class standard deviation (in red).

Figure 9 (

 9 Fig. 9. Retrieval results on McGill database [35].Figure 9(a): overview of the database. Figure 9(b) and 9(d): Average recall rate on articulated and complete dataset.Figure 9(c) and 9(e): Average Precision versus Recall on articulated and complete dataset. Figures 9(f) (from[START_REF] Siddiqi | Retrieving articulated 3-D models using medial surfaces[END_REF]) and 9(g)): comparison of our approach with stateof-the-art method described in[START_REF] Siddiqi | Retrieving articulated 3-D models using medial surfaces[END_REF] on partial McGill Dataset (14 classes out of 19). Results are shown as average recall rate curve (plotted in black) along with the average intra-class standard deviation (in red).

  Fig. 9. Retrieval results on McGill database [35].Figure 9(a): overview of the database. Figure 9(b) and 9(d): Average recall rate on articulated and complete dataset.Figure 9(c) and 9(e): Average Precision versus Recall on articulated and complete dataset. Figures 9(f) (from[START_REF] Siddiqi | Retrieving articulated 3-D models using medial surfaces[END_REF]) and 9(g)): comparison of our approach with stateof-the-art method described in[START_REF] Siddiqi | Retrieving articulated 3-D models using medial surfaces[END_REF] on partial McGill Dataset (14 classes out of 19). Results are shown as average recall rate curve (plotted in black) along with the average intra-class standard deviation (in red).

⋆ This work has been done with the support of the French "Agence Nationale de la Recherche" (ANR), under grant NatImages (ANR-08-EMER-009).