Andrei Zaharescu
email: andrei.zaharescu@aimetis.com

Edmond Boyer
email: edmond.boyer@inria.fr

Radu Horaud
email: radu.horaud@inria.fr

Keypoints and Local Descriptors of Scalar Functions on 2D Manifolds

Keywords: 3D shape, Meshed surfaces, Riemannian manifold, Scale space, 3D keypoint detection, Local shape descriptors, Shape matching

  L'archive ouverte pluridisciplinaire

Introduction

The representation of visual information in terms of a structured collection of local features has been an active research topic for the last decade and it is of great importance for a variety of tasks, such as tracking, registration, recognition, retrieval, etc. Feature-based approaches were introduced in the computer vision literature three decades ago [START_REF] Bolles | Recognizing and locating partially visible objects, the Local-Feature-Focus method[END_REF][START_REF] Bolles | 3DPO: A threedimensional part orientation system[END_REF] for the purpose of recognizing and localizing partially occluded objects. Initially, features represented local geometric information. More recently, feature-based image analysis has become very popular [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF][START_REF] Mikolajczyk | A performance evaluation of local descriptors[END_REF]. The vast majority of existing methods detect and describe features using photometric information from a single image. Recently, image features were extended to 2D+t features, used for characterizing short image sequences for video analysis [START_REF] Laptev | On space-time interest points[END_REF].

Recent progress in multiple-view stereo and imagebased modelling and rendering allows the recovery of geometric and photometric information directly from images (Seitz et al, 2006). This means that one can characterize 3D shapes based on both geometric and photometric features. However, if taken separately, geometric or photometric information have limited utility, as the whole richness of the data is not fully ex-ploited. Consider, for example, 3D deformable or articulated objects. Their 2D appearances, i.e., images, are not full invariant to motions and to viewing conditions and hence, image-based features do not yield robust 3D descriptors. Similarly, geometric features are not robust to complex object motions that can considerably change the topology. Therefore, we believe that photometric and geometric information need to be handled in a consistent and simultaneous manner. To this purpose, we observe that both photometric and geometric information available with a 3D object can both be viewed as scalar (real-valued) functions defined over a 2D manifold, e.g., the surface of an object; This may well be considered as a generalization of scalar functions defined over an Euclidean domain, e.g., light intensity over an image, to non-Euclidean domains. One can thus build on existing feature-based image description paradigms, in order to investigate and propose extensions to 3D shapes.

The main contribution of this paper is a novel family of keypoint-based local surface descriptors that takes into account both the geometric properties of the surface and any scalar field defined over the surface, e.g. the photometric information available from multiplecamera setups. We develop an interest point detector and associated local scale-space descriptors that can be applied to various function defined over the surface, e.g, texture, colour, Gaussian curvature, mean curvature, geodesic integral, etc. To this end, we use discrete operators, e.g., the gradient defined for scalar functions on discrete surface domains (meshes), thus taking into account both the functions' differential properties as well as the underlying intrinsic geometry of the surface. Based on these operators, both an interest point detector and a local descriptor are introduced, namely Mesh-DOG and MeshHOG. MeshDOG is a generalization of the Difference of Gaussian (DOG) operator [START_REF] Marr | Theory of edge detection[END_REF][START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF] and it is used to build a discrete Laplacian operator on a mesh. This allows us to represent scalar functions over multiple scales, i.e., by convolution with a discrete Laplacian operator, and to detect points of interest as local extrema. MeshHOG is a generalization of the Histogram of Oriented Gradients (HOG) descriptor that was proposed for describing 2D images [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF]. The newly descriptor, MeshHOG, is defined with respect to the measurements available at each vertex of the discrete surface and it can be implemented with any scalar function.

The newly proposed surface-based interest point detector and region descriptor exhibit a number of interesting properties with respect to the corresponding image-based counter-parts: 1. A 3D detector takes more information into account than an image based one and it does not suffer from false detections due to occluding contours. 2. There are no affine (nor perspective) distortions, since the computations are performed in a local metric space. Hence, there is no need for affine-invariance; 3. Image descriptors are sensitive to 2D occlusions. This is not an issue with surface descriptors, provided that a complete reconstruction of the underlying 3D object is available, which is typically the case with reconstructions from multiple camera setups; 4. The descriptor captures both the (local) 3D geometry and the gradient information associated with the scalar function; 5. In a multiple-camera setting, photometric information can be elegantly fused together from several images. This provides a photometric descriptor that is both image-invariant and more robust to image noise.

The remainder of the paper is organized as follows. Section 2 discusses related work and it emphasizes our own contribution to the problem of extracting local surface features. Section 3 introduces the mathematical formulation used to define the gradient and hessian operators on discrete manifolds that are needed to build local features. Section 4 and 5 introduce the local feature detector and region descriptor, respectively. Section 6 presents a detailed performance evaluation and comparison with other methods, before concluding in Section 7.

Related Work

In the recent past, there has been a lot of work aimed at visually characterizing 3D objects for the purposes of modelling and recognition.

Surface Features Based on Image Processing

One possible approach is to rely on image keypoints with their associated local descriptors and to use them in order to build 3D surface features.

Detectors. Keypoints are often associated with the detection of interest points in images, such as the extrema of the Laplacian of the image intensity function. They can be detected at various scales using the difference of Gaussian (DOG) approximation of the Laplacian (see [START_REF] Mikolajczyk | A performance evaluation of local descriptors[END_REF] for a review and the references therein).

Descriptors. 2D feature descriptors are generally designed to be robust to changes in illumination and invariant to image transformations such as translation, rotation, or scale [START_REF] Matas | Robust wide-baseline stereo from maximally stable extremal regions[END_REF][START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF][START_REF] Dufournaud | Image matching with scale adjustment[END_REF][START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF][START_REF] Bay | Speededup robust features (SURF)[END_REF] and, more generally, to 2D affine transformations [START_REF] Mikolajczyk | Scale & affine invariant interest point detectors[END_REF]. This type of imagebased descriptors have been successfully used to characterize 3D objects [START_REF] Rothganger | 3D Object Modeling and Recognition Using Local Affine-Invariant Image Descriptors and Multi-View Spatial Constraints[END_REF]. We note, however, that there are inherent limitations with these approaches. First, it is required that the image descriptors are back projected onto the 3D surface of the object which may lead to a redundant and ambiguous representation, since a 3D surface point corresponds to several points belonging to different images. Second, these descriptors are well suited only for objects which are locally planar. Third, the image-based descriptors are limited to photometric information and hence, one cannot build 3D descriptors based on the geometric properties of the underlying surface, e.g., curvature. Our method is quite different because it directly exploits the local geometric properties of the object's underlying surface.

Image keypoints and associated local descriptors were extended to image sequences [START_REF] Laptev | On space-time interest points[END_REF], (Wong and Cipolla, 2007), [START_REF] Kläser | A spatiotemporal descriptor based on 3D-gradients[END_REF], by considering the 2D+t (spatio-temporal) volume defined by a short image sequence. Such image+time features can be seen as local detectors/descriptors defined over a volumetric representation, i.e., a regular 3D grid. We consider a different problem, namely the extension from features defined on regular domains to features defined over irregular non-Euclidean domains.

Surface Features Based on Local Geometric Properties

Another category of methods attempts to extract local geometric information from range images or from pointcloud data.

Detectors. [START_REF] Novatnack | Scale-dependent 3D geometric features[END_REF] defines the scale space in a planar parameterization of the surface using the normal map; 3D keypoints are detected as the extrema of this representations, based on a gradient operator defined over a planar vector field. An analysis of the scale-variability of geometric structures captured in range images is proposed in [START_REF] Novatnack | Scaledependent/invariant local 3D shape descriptors for fully automatic registration of multiple sets of range images[END_REF], while [START_REF] Bariya | Scale-hierarchical 3d object recognition in cluttered scenes[END_REF] extends this method to deal with cluttered 3D scenes in an object recognition task. The automatic identification of interest regions on surfaces, taking into account geometric features such as scale-space extrema based on the average mean curvature flow, is proposed in (Schlattmann et al, 2008). In [START_REF] Mian | On the repeatability and quality of keypoints for local featurebased 3d object retrieval from cluttered scenes[END_REF], it is proposed to detect keypoints that characterize 3D and 2.5D surfaces. Mesh saliency methods are proposed in [START_REF] Lee | Mesh saliency[END_REF][START_REF] Castellani | Sparse points matching by combining 3D mesh saliency with statistical descriptors[END_REF], based on the centre-surround operator, adapted from the visual attention literature.

Descriptors. A number of 3D descriptors were proposed. 3D spin images [START_REF] Johnson | Using spin images for efficient object recognition in cluttered 3D scenes[END_REF], one of the first proposed methods, build 2D histograms by accumulating points that fall on a rotating plane along the normal. 3D shape contexts [START_REF] Körtgen | 3D shape matching with 3D shape contexts[END_REF][START_REF] Frome | Recognizing objects in range data using regional point descriptors[END_REF] extend the idea of spin images, by accumulating 3D histograms within a spherical support region. For a detailed survey, see (Tangelder and Veltkamp, 2004;[START_REF] Bustos | Feature-based similarity search in 3D object databases[END_REF]. Intrinsic shape signatures are proposed in (Zhong, 2009), thus improving on shape contexts by using a different histogram partitioning scheme. In [START_REF] Novatnack | Scaledependent/invariant local 3D shape descriptors for fully automatic registration of multiple sets of range images[END_REF], an image-based descriptor is proposed using the local R 2 embedding of the surface normal information. Similarly, [START_REF] Hua | Geodesic distance-weighted shape vector image diffusion[END_REF] builds 2D conformal maps of a 3D shape by mapping an irregular domain, such as a mesh, to a regular grid. In the resulting shape vector image, standard 2D scale-space keypoint detection and description [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF]) can be applied to build 3D features. Such features are robust to isomorphic deformations. However, due to conformal mapping limitations, they are sensitive to even small topological changes. More recently, (Ruggeri et al, 2010) proposed to build local 3D descriptors of a meshed surface at the locations of the critical points, the maxima, minima and saddle points of the eigenfunctions of the Beltrami-Laplace operator. Similarly, [START_REF] Bronstein | Shape Google: geometric words and expressions for invariant shape retrieval[END_REF] proposed to use the critical points of the auto-diffusion function, the diagonal elements of the heat-kernel matrix of a mesh, in order to obtain an intrinsic scale-space representation of the mesh's geometry. These methods allow to describe an object in terms of its intrinsic local geometry, but they do not allow to characterize other scalar functions. Again, these methods cannot exploit photometric information.

Combining Photometric Information with 3D Features

Closer to our own methodology, a category of methods attempts to characterize the photometric information that is available when considering 3D objects obtained from multiple-view reconstructions. (Starck and Hilton, 2007) proposes a concatenated surface descriptor, encompassing local geometry (a region descriptor based on geodesic-intensity histograms) and photometric information, edge and corner descriptors that take into account the local isometric mapping to R 2 ; In (Wu et al, 2008) a SIFT-based descriptor using 3D oriented patches is proposed, namely VIP, or viewpoint-invariant patches, which was used for 3D model matching. Both (Starck and Hilton, 2007) and (Wu et al, 2008) are among the first attempts to devise a descriptor that combines geometric and photometric information. Our proposed approach is similar in spirit to (Wu et al, 2008), but instead of back-projecting an image descriptor onto the surface, we propose to detect keypoints and build an associated descriptor directly onto the surface, taking full advantage of the 3D nature of the surface.

Recently, a number of extensions were proposed to our previous work (Zaharescu et al, 2009). (Smith et al, 2011) consider the specific scenario when the 3D triangulated mesh and the scalar function are from single range-intensity image. The discrete gradient approximation, inspired from (Xu, 2004), requires the scalar function to be defined over the faces of a triangulated mesh. This may be difficult and expensive to retrieve from the initial data, or it may not be available at all. Alternatively, we propose a radically different way of approximating gradient computation on any polygonal mesh that only requires that the scalar function is defined at the mesh vertices. This allows us to compute not only gradients of functions but also gradients of directional derivatives of these functions, which are needed for keypoint detection.

Registration and Recognition

Interestingly, many applications of 3D modelling make use of local features, e.g., rigid and non-rigid registration, object recognition, shape retrieval, etc. Recent work [START_REF] Furukawa | Dense 3D motion capture from synchronized video streams[END_REF][START_REF] Ahmed | Dense correspondence finding for parametrization-free animation reconstruction from video[END_REF][START_REF] De Aguiar | Marker-less 3D feature tracking for mesh-based human motion capture[END_REF]Varanasi et al, 2008) addressed non-rigid mesh registration using observations from multiple views. The vast majority of the proposed methods (the only notable exception being [START_REF] Furukawa | Dense 3D motion capture from synchronized video streams[END_REF])) use both geometric information extracted from surfaces and photometric data available from images. The latter is first extracted using 2D image descriptors (such as SIFT), and subsequently back-projected onto the mesh. This sparse description is generally used to bootstrap dense matching. Surface descriptors may well be used for 3D object recognition, as it has been already done in (Shilane et al, 2008), using the Princeton shape benchmarking database1 . Our work contributes to these efforts by taking a different, yet complementary approach: image-feature detection and description methodologies are extended to features defined directly onto discrete 2D manifolds.

Paper Contributions

This paper is an extended version of (Zaharescu et al, 2009), which introduced 3D shape descriptors inspired from image descriptors. Unlike an image, which is a regular Euclidean domain, a 3D shape is often defined over an irregular non-Euclidean domain such as a mesh, which may be viewed as a discrete manifold. In (Zaharescu et al, 2009) we proposed a surface keypoint detection based on a difference of Gaussian operator (MeshDOG) and a local descriptor based on the histogram of oriented gradients (MeshHOG). Both Mesh-DOG and MeshHOG require the estimation of a gradient operator and of first-and second-order directional derivative operators. In this paper we propose an improved computational framework for estimating the gradient and the directional derivatives of real-valued functions defined on discrete manifolds, e.g., a mesh. With respect to (Zaharescu et al, 2009), we relax the constraint that the mesh vertices must correspond to a regular sampling of the underlying continuous surface; We do this by using the geodesic distance throughout the formulation. Also, the proposed gradient computation is now cast into a least square minimization problem that can be efficiently estimated using a linear solver. The gradient method is inspired from [START_REF] Barth | A 3-D least-squares upwind Euler solver for unstructured meshes[END_REF][START_REF] Mukherjee | Learning gradients on manifolds[END_REF] and it could, in principle, be applied to any kind of polygonal mesh or to a point cloud. It differs from previous approaches that are based on the eigenfunctions of a discrete Laplace-Beltrami operator, e.g., (Xu, 2004), [START_REF] Luo | Approximating gradients for meshes and point clouds via diffusion metric[END_REF]. Additionally, we introduce a new dataset as well as an in-depth evaluation of our method, thus testing both the repeatability of the keypoint detector and the robustness of the region descriptor under a large number of deformations. Finally, two other existing datasets are used for comparisons with other existing methods.

Differential Mesh Processing

In this section we introduce a computational framework required to estimate interest points and local descriptors of a scalar function defined over a manifold. To this end, we define several operators that can handle an irregular domain, including the gradient operator and the first-and second-order directional derivative operators.

Let M be a 2D closed manifold (i.e. compact and without boundaries) embedded in R 3 and let M be a discrete mesh representation of M composed of vertices on M and of convex polygons, i.e. facets. M can be viewed as a graph M (V, E), where V = {v i } N i=1 is the set of mesh vertices and E = {e ij } is the set of edges between adjacent vertices. We associate a 3D point v i ∈ R 3 with each mesh vertex v i . Note that an image can be viewed as a "flat" uniformly sampled mesh with boundaries, i.e., a grid of vertices with valence 4 and whose facets are rectangles.

Gradient

Let f : M → R be a smooth real-valued function defined on M, e.g., photometric data or curvature. In order to estimate the gradient ∇ M f (v) of f at point v, we consider the first order Taylor expansion approximating f at a manifold point v i in a neighborhood of point v j :

f (v i ) ≈ f (v j ) + ∇ M f (v i ) ⊤ (v i -v j ).
(

) 1
where the gradient ∆ belongs, by definition, to the tangent plane of M at v i . In the discrete case one can therefore write:

∇ M f (v i ) ⊤ (v i -v j ) ≈ f (v i ) -f (v j ), (2)
where ∇ M f (v) denotes the discrete gradient of f at v. This expression can be used to estimate the discrete gradient at any mesh vertex v i through an error minimization criterion (Sibson, 1981). We adopt the least square gradient construction that follows this principle [START_REF] Barth | A 3-D least-squares upwind Euler solver for unstructured meshes[END_REF] and we seek the 3D vector that minimizes the criterion:

∇ M f (v i ) = argmin g vj ∼vi w ij (f (v i ) -f (v j )- g ⊤ (v i -v j )) 2 , (3)
where the notation v j ∼ v i means that v j ∈ N (v i ), i.e., the neighbourhood of v i considered in the estimation and where the weights w ij balance the contributions of the neighboring vertices. Both N (v i ) and w ij are choosen as follows.

N (v i ) is usually the first ring of vertices around v i . However, in order to make it more robust to nonuniform sampling, N (v i ) can be defined as the set of vertices v j ∈ M residing within a geodesic ball centred in v i of radius r:

N (v i ) = {v j |d g (v i , v j ) < r}, (4)
where d g (v i , v j ) represents the geodesic distance between v i and v j .

The weight function w ij can be uniform or it can vary with respect to, e.g., areas (Sibson, 1981) or inverse distances [START_REF] Barth | A 3-D least-squares upwind Euler solver for unstructured meshes[END_REF] in the neighbourhood of v i . In [START_REF] Mavriplis | Revisiting the least-squares procedure for gradient reconstruction on unstructured meshes[END_REF] it is shown that weighted gradient estimations based on inverse distances significantly improve over unweighted estimations. In this work, the weight function is a zero-centred Gaussian function:

w ij = G σ (d g (v i , v j )) = exp -d 2 g (v i , v j )/2σ 2 . ( 5
)
Note that vector g in ( 3) can be advantageously constrained to belong to the tangent plane of M at v i , whenever the normal unit vector n i to this tangent plane is known:

∇ M f (v i ) = argmin g vj ∼vi w ij (f (v i ) -f (v j )- g ⊤ (v i -v j )) 2 + λ i (g ⊤ n i ) 2 , ( 6
)
where the positive scalar λ i is chosen such that the tangent-plane constraint is emphasized:

λ i = vj ∼vi w ij .
The constrained minimization ( 6) is a linear least-squares problem that that is efficiently solved using standard matrix factorization methods, such as singular value decomposition [START_REF] Lay | Linear algebra and its applications[END_REF].

The current gradient computation formalism is better motivated mathematically than our previously proposed approximation (Zaharescu et al, 2009). Even though the current formulation handles better particular edge cases, both methods behave numerically similar on the average, when dealing with evenly sampled manifold discretizations.

Directional Derivatives

The directional derivative D a f (v) of f at v ∈ M along vector a is then, by definition, the projection of the gradient vector ∇ M f (v) along the direction of a:

D a f (v) = ∇ M f (v) ⊤ a a , ( 7
)
where a is a vector lying in the tangent plane of M at v. The discrete directional derivative at v i along any direction in the tangent plane at v i can therefore be computed using [START_REF] Dufournaud | Image matching with scale adjustment[END_REF] with the discrete gradient ∇ M f (v i ) estimated using (6).

Let us now consider second discrete directional derivatives along unit vectors a and b lying in the tangent plane to the mesh at v i . Such a second directional derivative can be written as:

D ab f (v i ) = ∇ M ∇ M f (v i ) ⊤ b ⊤ a, = D a (D b f (v i )) , (8)
which requires an estimate of the gradient of the scalar function

D b f (v i ), namely ∇D b f (v i ).
This can be easily obtained by applying the least square criterion (6) to the directional derivative function. Similarly, one can estimate

D aa f (v i ), D bb f (v i ), D ab f (v i ) and D ba f (v i ).
By further assuming that the vectors a and b form an orthonormal basis of the tangent plane at v i , one obtains the Hessian matrix of f :

H a,b (f (v i )) = D aa f (v i ) D ab f (v i ) D ba f (v i ) D bb f (v i ) . (9)
Assuming that the directional derivatives of f are continuous, the order of the differentiation does not matter, and hence (by Clairaut's theorem), one should expect the Hessian in ( 9) to be symmetric, namely that

D ab f (v i ) = D ba f (v i ).
However, in our case, the gradients ∇D a f (v i ) and ∇D b f (v i ) are obtained by numerical optimization of ( 6). Hence, the Hessian is not guaranteed to be symmetric. This means that it is not guaranteed that the two eigenvalues of ( 9) are real. Therefore, we propose to use the following real symmetric matrix:

Ha,b (f (v i )) = 1 2 H a,b (f (v i )) + H ⊤ a,b (f (v i )) ( 10
)
which corresponds to the projection of H onto the linear space of 2 × 2 symmetric matrices [START_REF] Horn | Matrix Analysis[END_REF].

Convolution

Finally, using the same notations, the normalized convolution of the function f with a Gaussian kernel G yields:

F σ (v i ) = f ⋆ G σ (v i ) = 1 K i vj ∼vi f (v j ) exp -d 2 g (v i , v j )/2σ 2 (11)
where G σ is the Gaussian function defined in ( 5) and K i is a normalization term such that:

K i = vj ∼vi G σ (d g (v i , v j )).
Using the properties of convolution, one can easily compute the first-and second-order directional derivatives of F σ , namely:

D a F σ (v i ) = D a f ⋆ G σ (v i ), ( 12
)
D ab F σ (v i ) = D ab f ⋆ G σ (v i ). ( 13
)

Numerical Approximations

Geodesics. The computation of geodesic distances on arbitrarily triangular meshes can be computed by the fast marching method [START_REF] Kimmel | Computing geodesic paths on manifolds[END_REF] or by other approximations (Surazhsky et al, 2005). In practice, in the interest of computational speed, we have used a local shortest path approximation on the edge connectivity graph. It has been observed experimentally that, for typical meshes, the variations due to the triangulation are minimal.

Normals. At first, a local normal estimation is used, using a one ring neighbourhood. In order to increase the robustness of the normal estimation, a smoothed version is then computed, using the mean estimate in a small geodesic neighbourhood.

Curvatures. Instead of using the classical curvature estimation method [START_REF] Meyer | Discrete differential geometry operators for triangulated 2-dimensional manifolds[END_REF] that employs only a one ring neighbourhood, we employ the more robust method proposed in [START_REF] Dong | Curvatures estimation on triangular mesh[END_REF], using a 3 ring neighbourhood.

Feature Detection (MeshDOG)

Feature detection comprises three steps, as illustrated in Figure 2. First, the extrema of the function's Laplacian are found across scales using a one-ring neighbourhood. The Laplacian is approximated with the standard difference of Gaussian (DOG) operator. Second, the extrema thus detected are thresholded. Finally, the unstable extrema are eliminated, only retaining the features exhibiting some degree of cornerness.

Scale-Space Construction

A scale-space representation of any scalar function f defined on a mesh is considered, built by progressive convolutions over f . The scale-space is built over s = 3 octaves, covering each octave in c = 6 steps. This is accomplished by progressive convolutions with a Gaussian kernel (11) of the original scalar function:

F 0 = f, (14)
F t = F t-1 ⋆ G σ(t) , (15)
with t = {1, 2, . . . , s•c}. The standard deviation parameter σ(t) of the Gaussian at iteration t is chosen based on the formula:

σ(t) = 2 1 c-2 ⌈ t c ⌉ e avg , (16)
where e avg represents the average edge length. As it can be observed, the value of σ(t) only changes when t starts spanning a new octave. σ(t) remains unchanged for the next c iterations, while t covers the current octave, thanks to the term ⌈ t c ⌉. This behaviour emulates in spirit the 2D grid downsampling, introduced in [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF], but without modifying the mesh geometry, which can be an expensive operation, in the case of non-uniformly sampled meshes. Therefore, an important observation is that, when building the scale space of scalar function defined over the mesh, the mesh geometry does not change. This contrasts to other approaches, such as [START_REF] Hou | Efficient computation of scalespace features for deformable shape correspondences[END_REF], that construct the scale-space by generating meshes with different samplings, thus requiring further mesh processing and simplification.

The difference of Gaussian operator is then used as an approximation of the Laplacian operator, built by subtracting adjacent convolved functions:

L t = F t -F t-1 . ( 17
)
An example can be observed in Figure 3, where the data being used correspond to frame 30 of the pop2lock sequence from the University of Surrey data set. The features being shown are colour intensity (first row) and mean curvature (second row).

Feature Detection

Feature points represent a subset of all vertices that can be detected with high repeatability. Using local gradient information is one way to detect a repeatable feature. Therefore, the feature points are selected as the local extrema over one ring neighbourhoods, in the current and in the adjacent scales. Such an example can be observed in Figure 1(b).

From the extrema of the scale space, only the top β = 5% of the maximum number of vertices are being considered, sorted by magnitude. We have chosen a percentage value, versus a hard value threshold, in order to keep the detector flexible, no matter which scalar function is being considered, i.e. colour intensity or mean curvature, without the need for normalization. However, when the threshold response is known a priori for a particular scalar function, such as it is the case in [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF] with image intensity, it can be easily used instead.

Additionally, in order to eliminate more non-stable responses, we only retain the features that exhibit corner characteristics. As proposed in [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF], this (a) color can be done by examining the eigenvalues of the 2 × 2 Hessian matrix of second directional derivates of the difference of Gaussian operator, i.e., Section 3.

(b) f 2 (c) f 18 (d) L 2 (e) L 18 (f) curvature (g) f 2 (h) f 18 (i) L 2 (j) L 18
Let us consider Difference of Gaussian operator L i , defined in (17) and the symmetric Hessian matrix approximation Hx,y (L i (v i )) from (10), where (x, y) is a pair of orthonormal vectors lying in the tangent plane T i of M at v i . The absolute value of the ratio between the two sorted eigenvalues, |µ 1 | ≥ |µ 2 | of this matrix is a good indication of a corner response. By construction, this ratio is independent of the choice of the local coordinate frame, i.e., vectors x and y. We use |µ 1 /µ 2 | = 10 as threshold value to eliminate the non-stable edge responses.

Feature Descriptor (MeshHOG)

In association with the detector presented in the previous sections, we propose building a local descriptor, named MeshHOG, similar in spirit to the histogram of gradient descriptor (HOG) [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF], but extended to 2D manifolds. A 2D image is in essence a 2D regular grid. The regularity assumption does not always hold in the case of a 2D manifold, that can exhibit non regular sampling. For this reason, the support region of the descriptor has to be chosen using a measure invariant to local triangulation, such as the geodesic distance. In addition to invariance to mesh sampling, the descriptor should also exhibit invariance to a number of other transformations, such as rotation and scale. The scale invariance is achieved by consider-ing the gradient information at the scale of the detected interest point. Rotation invariance is achieved by defining a local coordinate system using the normal at the detected interest point, the dominant gradient in the support region and their cross product. Finally, a two level histogram of gradient is computed, both spatially, at a coarse level, in order to maintain a certain highlevel spatial ordering, and using orientations, at a finer level. Trilinear interpolation is used in order to minimize the sampling effect in the histogram bins. The fact that the gradient vectors are 3D allows the computation of the histograms in 3D.

Support Region

The descriptor for vertex v i is computed within a support region N (v i ), defined using a geodesic ball of radius r, as in ( 4). The geodesic support region is chosen adaptively based on a global measure, such that the descriptor is robust to scale and to spatial sampling. The value of r is chosen such that it covers a proportion α r of the the total mesh area, where α r ∈ (0, 1). By denoting with A M the total area of the mesh M , which can be computed as the sum of all triangle areas, the radius of the circular support region is:

r = α r A M Π , (18)
where [x] denotes the integer of x and assuming that the surface covering the ring neighbourhood can be approximated with a circle. In practice, we use r such that α r = 2%.

In the current work we assumed that we are dealing with fully reconstructed objects, thus recovering the "true" global object scale. If this is the case, choosing the size of a support region based on (18) makes it scale and sampling invariant. However, when dealing with partially reconstructed objects, the support-region size should be chosen based on the scale reported by the detected interest point, as suggested in [START_REF] Novatnack | Scaledependent/invariant local 3D shape descriptors for fully automatic registration of multiple sets of range images[END_REF][START_REF] Bariya | Scale-hierarchical 3d object recognition in cluttered scenes[END_REF][START_REF] Mian | On the repeatability and quality of keypoints for local featurebased 3d object retrieval from cluttered scenes[END_REF].

Local Coordinate System

As mentioned earlier, a local coordinate system is desirable, in order to make the descriptor invariant to mesh rotations. A local coordinate system can be built using the unit vector n i orthogonal to the plane T i tangent to M at vertex v i , and a pair of orthonormal vectors residing in this plane. Given an arbitrary unit vector a i ∈ T i , a local coordinate system C is constructed as:

C = {a i , n i , a i × n i }. ( 19
)
It is therefore important to choose the unit vector a i based on some intrinsic local property of the scalar function, and hence make the choice of the local coordinate system invariant to mesh rotations. The direction corresponding to the most dominant gradient magnitude in the neighourhood exhibits such a desired behaviour. Therefore, the unit vector a i is chosen as the direction associated with the dominant bin in a polar histogram h a , with b a = 36 bins. The histogram is constructed by considering the projected participating neighbouring vertices v j ∈ N (v i ) onto T i . The vertex contribution c i (v j ) to the appropriate bin takes into account the gradient magnitude and the geodesic distance from the centre vertex v i , weighted by a Gaussian (see (5)):

c i (v j ) = ||∇ M f (v j )|| G σ (d g (v i , v j )), ( 20
)
where the standard deviation is: σ = ǫr, with r the support region radius and ǫ the spatial influence. ǫ is set to 0.5 in our experiments.

Therefore, bin h a (k, v i ) yields:

h a (k, v i ) = vj ∼vi χ ha(k,vi) (v j ) c i (v j ), (21)
where k = {1, 2, . . . , b a } is the bin index and χ ha(k,vi) (v j ) represents the indicator function for bin selection. In general, the indicator function is defined as :

χ Z (x) = 1 if x ∈ binZ, 0 otherwise. ( 22
)
In order to reduce aliasing and the boundary effects of binning, votes c i (v j ) are interpolated trilinearly between neighbouring bins during the histogram computation. In practice, that means relaxing the indicator function definition ( 22). The same interpolation technique is used in the next sub-section.

The HOG Descriptor

Instead of computing full 3D orientation histograms, as proposed in [START_REF] Kläser | A spatiotemporal descriptor based on 3D-gradients[END_REF], we project the gradient vectors onto the three planes associated with the local coordinate system (19), in order to provide a more compact representation of the descriptor. A possible drawback of the approach described below is a decrease in the discriminability of the descriptor. However, given the fact that local surface neighbourhoods and the associated gradients estimated at the mesh vertices typically span a limited subset of the possible 3D spatial/orientation bins, the current compression scheme does not incur any practical setbacks.

For each of the three planes, a two-level histogram h s,o is computed. First, the plane is divided in b s = 4 polar slices h s , starting with the origin and continuing in the direction dictated by the right hand rule, with respect to the other orthonormal vector. When projected onto the plane, the participating neighbouring vertices v j will fall within one of the spatial slices.

Second, for each spatial slice, the space is divided into b o = 8 orientation slices h o . The projected gradient vectors ∇ M f (v j ) of the vertices v j ∈ N (v i ) that projected onto spatial slice are used to determine the orientation slice, as shown in Figure 4(d).

Similar to the histogram definition (21), the histogram bin h s,o (e, l, v i ) yields:

h s,o (e, l, v i ) = vj ∼vi χ hs,o(e,l,vi) (v j ) c i (v j ), (23)
where e = {1, 2, . . . , b s } is the spatial bin index, l = {1, 2, . . . , b o } is the orientation bin index and χ hs,o(e,l,vi) represents the indicator function for bin selection, defined as:

χ hs,o(e,l,vi) (v j ) = χ hs(e,vi) (v j ) χ ho(l,vi) (∇ M f (v j )). ( 24
)
The final descriptor is obtained by concatenating the b s • b o histogram values for each of the 3 orthonormal planes. In order to make the descriptor invariant to mesh sampling, the concatenated histograms are normalized using the L2 norm. The final descriptor will have 3 × b s × b o elements. Given the previous choice of parameters, the dimensionality of the descriptor is 3 × 4 × 8 = 96.

Whenever reducing the descriptor dimensionality is a requirement, keeping only histograms computed in the tangent plane T i is a possibility, thus shrinking the descriptor to 32 elements. Also, if multiple scalar functions are available, an aggregate descriptor can be built by concatenating multiple individual descriptors.

The method has been implemented in C++. The source code has been made available under a GPL license and it can be downloaded from 2 .

Performance Evaluation

In this section an extensive evaluation will be performed, covering both the interest point detector and the region descriptor. The original paper (Zaharescu et al, 2009) introduced some preliminary empirical results, in the context of sparse mesh matching.

Results are presented in four different scenarios. In Section 6.1, the performance is evaluated on the newly proposed PHOTOMESH dataset, which consists of deformations on meshes equipped with photometric information. In section 6.2 the performance of the currently proposed method is compared with other state of the art methods on the SHREC 2010 features database [START_REF] Bronstein | Shrec 2010: robust feature detection and description benchmark[END_REF], containing only geometric deformations. Additional comparisons with other methods are presented in Section 6.3, using the database introduced in [START_REF] Kovnatsky | Photometric heat kernel signatures[END_REF], containing both geometric and photometric deformations, but in the context of mesh retrieval. Finally, additional results are presented in Section 6.3 in the context of rigid and non-rigid sparse mesh matching.

Benchmark on the PHOTOMESH Dataset

While there already exist datasets that test a number of geometrical deformations (see the next section), they do not contain photometric meshes and limit scalar functions to measures of the surface geometry alone. Therefore, a new dataset is proposed, named PHOTOMESH, aimed at testing the repeatability of the detector and the robustness of the descriptor under both photometric and geometric deformations.

Dataset

The dataset consists of 3 base shapes. also called null shapes, endowed with photometric information at each vertex. Simulated transformations are applied to them. Two of the null meshes are obtained from multi-view stereo reconstruction algorithms, whereas one is generated with a modelling program. The photometric transformations are noise and shot noise. The geometric transformations are noise, shot noise, rotation, scale, local scale, sampling, holes, micro-holes, topology and isometry. Each transformation has 5 levels of noise applied to it. Therefore, for one base shape, a total of 13•5 = 65 shapes are obtained. Hence, the database contains 135 shapes.

Noise. Generally, the noise level corresponds to the noise amplitude. For more information on how the noise is generated for the various transformations, please consult Table 1. For isometries, groundtruth is obtained either by means of non-rigid semi-elastic transformations [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF], when using real data, or manually, in a modelling program, such as Blender, when using fully synthetic data. Additional noise is inherently introduced by the multi-view image-based mesh reconstruction and mesh colour estimation process. In addition to the colour differences between views, there are regions that are not visible in any camera view, but which are still reconstructed in 3D due to the interpolation process, such as the sole of a foot.

Evaluation Methodology

The evaluation was performed for feature detection and for region description. The performance is measured by comparing the features and descriptors obtained for the null shape with the ones obtained for the different transformations.

Feature Detection. The criteria employed for quantifying the quality of feature detection is repeatabil-ity. Given that the ground-truth (one-to-one correspondence) is known for each transformed shape B of the null shape A, the repeatability is calculated as the percentage of detected feature points in B that are within a geodesic ball of radius r = 1% of the surface area, from one of the detected interest points in A.

Feature Description.The quality of the feature description is measured as the average normalized L 2 distance between descriptors corresponding to matched feature points.

Results and Discussion

Results are presented for three different scalar fields, defined over the manifolds: colour intensity (f I ), mean curvature (f M ) and Gaussian curvature (f G ). Tables 345summarize the repeatability of the detector, whereas Tables 678show the robustness of the descriptor.

In the case of colour noise and color shot noise performance slowly degrades if the scalar function utilizes the colour information (f I -Tables 3 and6). In the case of noise and shot noise, the performance slowly degrades for all functions. Notice, however, that f I leads to better results than f M and f G , because the last two scalar functions use geometric information both when computing the gradient and the actual curvature function.

Rotation and scale transformations prove experimentally that the method is indeed invariant to rigid transformations.

Holes, micro-holes and topological transformations affect the performance of the method linearly. That is partly because these transformations modify the total surface area, which in turn changes the local support area of the descriptor. Also, topological changes introduce new keypoints with the new structures. In the case of holes, some of the keypoints are simply missing.

Even though some invariance to the density of the mesh discretization is built into the method, sampling still affects the performance. One aspect is related to the fact that, even though the detection method has invariance built in, during the re-sampling process (when generating the transformation) some of the features could have moved further away than 1% and they are now detected as incorrect matches. Also, resampling will affect the computation of curvature. Lastly, when computing the descriptor histograms, as the sampling decreases, more bins can become empty, due to the increased sparsity of the participating vertices. One potential way to overcome this effect is to ensure sufficient sampling, by Gaussian Noise with σ x = {0.002t, 0.005t, 0.01t, 0.02t, 0.05t}, t = 255 Shot Noise Color Shot noise with signal to noise ratio SN R x = {0.002, 0.005, 0.01, 0.02, 0.05} and noise amplitude σ = 50. Noise Geometry Gaussian Noise with σ x = {0.1t, 0.2t, 0.3t, 0.4t, 0.5t}, where t = e avg . Shot Noise Geometry Shot noise with signal to noise ratio SN R x = {0.002, 0.005, 0.01, 0.02, 0.05} and noise amplitude σ = 20 e avg . Rotation Geometry Gaussian noise with σ x = {0.1t, 0.2t, 0.3t, 0.4t, 0.5t}, t = π.

Scale Geometry Scale factor s x = {0.5, 0.83, 1.25, 1.62, 2.0}. Local Scale Geometry The input mesh is dilated 3 * x times. At each iteration, the vertices are moved along the their normal by e avg /3. Sampling Geometry The mesh is resampled to attain a desired edge size e avg (1.0 + x), using edge split, edge collapse and edge swap operations, as described in (Zaharescu et al, 2011). Holes Geometry x random round holes are created, each having area corresponding to 5% of the total initial mesh area. Micro-Holes Geometry 3 * x random round micro-holes are created, each having an area corresponding to 3 neighbouring rings from the chosen centre vertex. Topology Geometry The input mesh is sliced with x equidistant planes into closed non-connected components, using (Zaharescu et al, 2011).

Isometry & Noise Mixed

Meshes are chosen from the captured 3D temporal sequences. x does not encode noise amplitude. Noise is inherently introduced by the multi-camera mesh reconstruction method and the multi-image colour estimation process.

Table 1 Transformations and noise levels of the PHOTOMESH dataset.

sub-sampling densely enough to guarantee that at least one sample is available for each bin.

The method also exhibit quasi-invariance to isometric transformations. In theory, pure isometric transformations do not affect Gaussian curvature, but they do affect mean curvature and the estimation of normals. In practice, however, the transformations are not purely isometric, due to errors introduced during the mesh tracking or mesh deformation process. In addition, when using real meshes obtained from the multiview stereo reconstruction, a significant amount of colour noise is implicitly introduced, since colours are interpolated in areas that are non-visible. When considering just the synthetic mesh isometry, the results are a lot more accurate, as presented separately in Table 2. Overall, the best results are obtained when using the detector / descriptor in conjunction with the photometric information (Tables 3 and6).

Measure f I f M f G Repeatability 0.

Strength Transform.

1 < 2 < 3 < 4 < 5 Color Noise 0.00 0.00 0.00 0.00 0.00 Color Shot Noise 0.00 0.00 0.00 0.00 0.00 Geometry Noise 0.26 0.29 0.31 0.33 0.34 Geometry Shot Noise 0.04 0.09 0.14 0.21 0.29 Rotation 0.01 0.01 0.01 0.01 0.01 Scale 0.01 0.01 0.01 0.01 0.00 Local Scale 0.21 0.25 0.28 0.30 0.31 Sampling 0.31 0.34 0.34 0.36 0.36 Holes 0.02 0.02 0.07 0.07 0.07 Micro-Holes 0.01 0.01 0.07 0.07 0.08 Topology 0.13 0.20 0.22 0.25 0.28 Isometry + Noise 0.23 0.24 0.22 0.25 0.25 Average 0.10 0.12 0.14 0.15 0.17

Table 8 Robustness of MeshHOG (Gaussian curvature) 2010), in order to be able to compare with other existing methods.

Dataset

The SHREC 2010 feature dataset is similar to the previously proposed benchmark dataset from Section 6.1, except for the fact that it does not contain scalar fields defined over the manifolds (i.e. textures), nor any of scalar function transformations. In addition, the null shapes in the current dataset are all synthetic, whereas 66% of null shapes from thePHOTOMESH dataset are obtained from multi-view stereo. The geometric transformations are similar to the ones introduced in the Section 6.1.1. An example can be seen in Figure 6 6.2.

Evaluation Methodology

The evaluation methodology is similar to the one presented in Section 6.1.2. The measures tested are the repeatability of feature detections and the robustness of the descriptors.

Results

The results of the proposed method are presented in Tables 9 and14. The Gaussian curvature was used as the scalar function. The results using mean curvature are very similar, but they are omitted in the interest of space.

Results for other top-performing methods, taken from [START_REF] Bronstein | Shrec 2010: robust feature detection and description benchmark[END_REF], are also included. For feature detection, the following other methods are included: two heat kernels variants based on the work of (Sun et al, 2009) -HK1 (Table 10) and HK2 (Table 11); a Harris 3D corner detection method (Sipiran and Bustos, 2010) -H1 (Table 12) and a saliency based method [START_REF] Castellani | Sparse points matching by combining 3D mesh saliency with statistical descriptors[END_REF]) -SP3 (Table 13). For region description, the following methods are included: two sparse heat kernel signature variants based on the work of (Sun et al, 2009) -SHK1 (Table 15) and SHK2 (Table 16) and the spin images method [START_REF] Johnson | Using spin images for efficient object recognition in cluttered 3D scenes[END_REF] -SI (Table 17). For more information about the competing methods and for additional results, please consult [START_REF] Bronstein | Shrec 2010: robust feature detection and description benchmark[END_REF].

In the context of feature detection, the proposed MeshDOG method performs very well: it has the top average results for noise levels 3 to 5 and scores slightly worse (< 2%) than HK1 and HK2 for noise levels 1 and 2. In the context of region description, MeshHOG's performance is affected by sampling and noise, which pulls down the average performance. As mentioned earlier, the method is not fully robust to sampling errors, mostly due to the fact that for very sparse samplings, a large number of bins are left empty. Proper dense resampling of the input mesh would alleviate this problem.

Strength

Benchmarking Using the SHREC Photometric Dataset

A new dataset of photometric meshes is proposed in [START_REF] Kovnatsky | Photometric heat kernel signatures[END_REF] in the context of shape retrieval.

Dataset

The query set consists of 270 real-world human shapes from 5 classes, obtained from multi-view stereo reconstruction, to which a number of transformations have been applied. Geometric transformations are divided into isometry+topology (real articulations and topological changes due to acquisition imperfections), and partiality (occlusions and addition of clutter). Photometric transformations include contrast, brightness, hue, saturation and color noise. Mixed transformations include isometry+topology transformations in combination with two randomly selected photometric transformations. For each class, there are 5 different transformation strength levels, adding up to 54 instances per shape.

The null shape of each of the 5 classes was added to the queried corpus, in addition to the other 75 shapes used as clutter. Retrieval was performed by matching 270 transformed queries to the 75 + 5 null shapes. Each query had exactly one correct corresponding null shape in the dataset.

Evaluation Methodology

Performance was evaluated using the precision-recall characteristic for the shape retrieval task. Precision P (r) is defined as the percentage of relevant shapes in the first r top-ranked retrieved shapes. Mean average precision (mAP ) was used as a single measure of performance, where mAP = r P (r) • rel(r) and rel(r) is the relevance of a given rank.

Results

Results for the proposed method are summarized in Table 18. Photometric information was used in the scalar field.

Results for other top-performing methods are also presented in Tables 192021222324. For more information, please consult [START_REF] Kovnatsky | Photometric heat kernel signatures[END_REF]. The methods using bag of features with heat kernel signatures (Tables 19) and the spectral distance (Table 20) are purely geometric, which makes them automatically invariant to the photometric noise (contrast, brightness, hue, saturation, noise). Conversely, the colour histogram method (Table 21) uses purely photometric information, thus making it invariant to purely geometric transformations (isometry + topology and partial). The methods presented in Tables 22 -24 use spectral decomposition methods in conjunction with both photometric and geometric information.

Comparing the overall results, the proposed method (Table 18) performs very well. It provides the best results for the partial geometric transformations (occlusions and the addition of clutter), as well as the best overall results for noise strength 1. While overall top ranking results are obtained by the method described in Table 23, the proposed method comes to a very close second. Table 22 Performance (mAP in %) of BoFs with cHKS descriptors using fixed colorspace scaling factor (w = 0.1).

Strength

Shape Matching

The current subsection introduces results of the currently proposed approach in a shape matching application. Finally, the candidate match is only accepted if the second best match is significantly worse (γ = 0.7 or less from the best match score). This method is not meant to fully solve the matching problem, as would a global approach (Starck and Hilton, 2007). It is intended to allow further validation and evaluation of the proposed detector and descriptor.

Datasets

In the evaluation, the following scenarios are considered: (i) the two meshes are representations of the same rigid object, which can thus be aligned using a rigid transformation; (ii) the two shapes are representations of the same non-rigid object, i.e. a moving person. In this context, the following datasets are introduced:

-Rigid Objects: Reconstructions of the same object are considered, using different camera sets. In particular, different mesh reconstructions are obtained using (Zaharescu et al, 2011) Photometric information. The colour of each vertex of the surface is computed by considering the median colour in the visible images. It is assumed that the colours of a vertex follows a non-Gaussian distribution, due to errors that can occur around occluding contours. In the Synth-Dance dataset the vertices are randomly coloured.

Evaluation Methodology

For each of the cases, a number of matches are produced by the above mentioned matching algorithm. In the case of rigid objects and synthetic non-rigid objects, the match groundtruth is readily available. Therefore, error distributions of the matches can be computed. They accumulate binned error distance from the groundtruth match, with the size of a bin being the average edge length.

Comparison with 2D SIFT. In addition, for the non-rigid motion cases, a comparison is presented between the proposed mesh matching framework using the MeshHOG descriptor and a similar matching framework, based on 2D image detectors and descriptors backprojected onto the mesh. In the image-based framework, the matching is performed in images and only then the 2D matches are back-projected onto the surface. The SIFT [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF] keypoint detector and image descriptor was used in 2D. When matching the two surfaces, only matches from the same cameras are considered. In order to carry such a comparison for the Synth-Dance dataset, 16 virtual cameras and the associated images have been generated. The virtual cameras are distributed in a circular pattern around the object.

Results

Rigid Matching Results. Figure 7 presents results on the Dino and Temple datasets, with different possible combinations of scalar functions for both detection and description: photometric information, mean curvature and Gaussian curvature. As the results show, the biggest benefit is obtained from using the photometric information for region description, irrespective of which scalar function is used for detection, e.g. the red curves in the graphs from Figure 7. However, the best overall results are obtained when photometric information is used for both keypoint detection and region description. Mean curvature seems to be the second most informative measure when used for region description (the green curves in the graphs from Figure 7),whereas the Gaussian curvature (the blue curves) is the worst performer.

Given the datasets and the matching procedure, it is to be expected that photometric information provides the best choice for region description. Geometry alone does not provide sufficient unique regions. Both the Dino and Temple meshes exhibit a large number of repetitive geometric features, that can only be disambiguated due to the slightly different photometric information.

Non-Rigid Matching Results. Synthetic comparative results are presented in Figure 8. The mesh in the first frame was matched with all the other 199 meshes across the sequence. Observe that the Mesh-HOG descriptor generates very few false positives in comparison with the SIFT equivalent, clearly demonstrating the advantages of the proposed approach. In addition, empirical results are presented in Figure 9 for the INRIA Dance-1 sequence. There are only 54 matches found using the SIFT back-projected method between frame 525 and 526, whereas MeshHOG finds 174 matches. Even when matching across distant frames (530 and 550), our proposed method finds 27 correct matches, whereas the SIFT back-projected method fails completely. It is to be expected, since most of the interframe matches are due to local creases formed by the clothes.The head is the only unique feature that can be robustly matched across time.

Conclusion

We have introduced MeshDOG and MeshHOG, a new 3D interest point detector and a new 3D descriptor, defined on triangular meshes endowed with a scalar function. The descriptor is able to capture the properties of both the local geometry and of the scalar function in a succinct fashion. It is robust to changes in orientation, rotation, translation and scale. The performance of both the interest point detector and the feature descriptor was tested extensively, achieving very competitive results, comparable with the state of the art.

As a future direction, we plan to investigate how to further extend the descriptor to take into account the temporal dimension, considering the context of densely tracked surfaces.

Fig. 1

 1 Fig.1The feature detection method described in this paper can be applied to any scalar function defined over a 2D manifold such as the meshed surface shown here: photometric data (a) and associated points of interest (b); mean surface curvature (c) and the detected features (d).

Fig. 2

 2 Fig. 2 Feature detection shown with photometric data. (a) Original mesh (27240 vertices); (b) Scale-space extrema (5760 vertices left); (c) Thresholding (1360 vertices left); (d) Corner detection (650 vertices left).

Fig. 3

 3 Fig. 3 An example of processing (a) color intensity and (f) mean curvature. Scale-space representation of color (b)-(e) and of mean curvature (g)-(j). The mesh corresponds to frame 30 of the pop2lock sequence from the University of Surrey.

Fig. 4

 4 Fig. 4 a) 3D Histogram -polar mapping used for creating histograms via binning of 3D vectors; b) Choosing three orthogonal planes onto which to project the 3D histogram. c) Polar coordinate system used for creating histograms via binning of 2D vectors, shown in this example with eight polar slices. d) Example of typical spatial and orientation histograms, using four spatial polar slices and eight orientation slices.

Fig. 5

 5 Fig. 5 Examples of possible transformations of the null shape (shown in strength 3 out of 5) for the PHOTOMESH dataset.

Fig. 6

 6 Fig. 6 Examples of possible transformations of the null shape (shown in strength 5 out of 5) for the SHREC 2010 Features dataset. Image taken from[START_REF] Bronstein | Shrec 2010: robust feature detection and description benchmark[END_REF]

Fig. 9

 9 Fig. 9 Non Rigid matching using real data -Dance-1 sequence. (a) Matches between frames 525 and 526 using Mesh-HOG (174 matches); (b) Matches between frames 530 and 550 using MeshHOG (27 matches); (c) Matches between frames 525 and 526 using SIFT (54 matches); (d) Matches between frames 530 and 550 using SIFT (0 matches).

Fig. 7

 7 Fig. 7 Rigid matching results -Dino and Temple datasets. (a) (c) Results when using photometric information for both detection and description; (b) (d) Error distribution when using different combinations of scalar functions.

Fig. 8

 8 Fig. 8 Non Rigid matching using synthetic datadancer-synth dataset. Comparison between MeshHOG and SIFT matching results. Matches between frames 1 and 50 are visually depicted in (a),(d). There are 364matches for MeshHOG and 343 matches for SIFT. They are also quantified in the error histograms (b),(e). The histogram bins are of size equal to e avg . The last bin groups all the errors greater than 20 * e avg . Additionally, the average histogram errors are shown in (c),(f) for matching frame 1 with x, where x ∈ [2..200].

Table 2

 2 Performance of different scalar functions under purely isometric transformations.

	99 0.98 0.98
	Robustness 0.07 0.15 0.16

Table 9

 9 Repeatability of MeshDOG: feature detection algorithm using Gaussian curvature as the scalar field. Average number of detected points: 129. Topology 97.44 96.10 92.26 91.22 88.64 Holes 91.48 90.60 86.78 83.73 81.86 Micro holes 98.08 96.69 96.00 95.52 94.87 Scale 99.36 99.36 98.50 97.90 97.68 Local scale 98.08 94.83 90.09 83.05 78.31 Sampling 97.05 97.88 97.39 96.27 92.35 Noise 95.30 92.78 91.67 89.24 87.62 Shot noise 98.08 96.22 93.39 90.45 87.32 Average 96.99 95.91 93.79 91.70 89.63 Table 10 Repeatability of HK1: heat kernel based feature detection algorithm. Average number of detected points: 23.

	Transform.	1	≤2	≤3	≤4	≤5
	Isometry	97.44 98.72 98.03 98.49 98.49
	Topology	97.44 97.44 97.44 97.40 97.41
	Holes	96.50 96.50 96.26 95.91 95.55
	Micro holes	97.31 97.24 97.22 97.08 96.95
	Scale	97.44 97.44 97.35 97.24 97.18
	Local scale	94.62 91.67 89.27 85.99 82.62
	Sampling	88.08 84.94 81.20 77.82 72.92
	Noise	91.92 91.92 90.09 88.59 87.10
	Shot noise	97.44 97.50 97.44 97.40 97.38
	Average	95.35 94.82 93.81 92.88 91.73

Table 14

 14 Robustness of MeshHOG feature description algorithm based on features detected by MeshDOG (average L 2 distance between descriptors at corresponding points). Average number of points: 129.

				Strength
	Transform.	1	≤2	≤3	≤4	≤5
	Isometry	0.05 0.04 0.04 0.04 0.04
	Topology	0.05 0.06 0.12 0.14 0.19
	Holes	0.07 0.07 0.07 0.08 0.09
	Micro holes	0.05 0.05 0.06 0.06 0.06
	Scale	0.05 0.05 0.05 0.05 0.05
	Local scale	0.07 0.09 0.10 0.12 0.14
	Sampling	0.06 0.06 0.06 0.08 0.13
	Noise	0.08 0.09 0.11 0.12 0.13
	Shot noise	0.05 0.06 0.10 0.16 0.25
	Average	0.06 0.06 0.08 0.09 0.12
	Table 15 Robustness of SHK1: heat kernel signature feature
	description algorithm based on featured detected by HK1
	(average L2 distance between descriptors at corresponding
	points). Average number of points: 23.

http://shape.cs.princeton.edu/benchmark/

http://mvviewer.gforge.inria.fr

http://tosca.cs.technion.ac.il/book/shrec_feat. html

http://4drepository.inrialpes.fr/

Acknowledgements. We would like to thank Cedric Cagniart and Artiom Kovnatsky for their help with the datasets.

Strength

Transform.

1 < 2 < 3 < 4 < 5 Color Noise 1.00 0.99 0.99 0.97 0.93 Color Shot Noise 0.98 0.96 0.91 0.86 0.76 Geometry Noise 1.00 1.00 1.00 0.99 0.99 Geometry Shot Noise 1.00 0.99 0.99 0.99 0.98 Rotation 1.00 1.00 1.00 1.00 1.00 Scale 1.00 1.00 1.00 1.00 1.00 Local Scale 1.00 1.00 0.99 0.99 0.99 Sampling 0.96 0.96 0.95 0.90 0.94 Holes 1.00 1.00 0.99 0.99 0.97 Micro-Holes 1.00 1.00 0.99 0.99 0.99 Topology 0.93 0.86 0.82 0.82 0.78 Isometry + Noise 0.95 0.97 0.97 0.93 0.96 Average 0.98 0.98 0.97 0.95 0.94