1 Introduction

Finding correspondence among three-dimensional shapes undergoing non-rigid transformations is one of the most fundamental recent problems in computer vision. It is actually a highly active research area, since it represents a key task in diverse applications such as motion tracking and recognition, shape interpolation and morphing, space-time reconstruction, shape retrieval and videos indexing.

Formally, matching pairs of shapes consists on establishing a map \(f: S \rightarrow T\) between two given surfaces S and T which are semantically equivalent or their geometrical properties are similar.

Finding a correspondence between shapes in the rigid transformations case has been efficiently treated. It consists on the estimation of a rotation and a translation. However, matching two shapes undergoing non-rigid deformations, still, remains a challenging problem. Unlike the rigid case, the non-rigid one involves an important number of freedom degrees increasing with the number of matched points and thus, estimating these non-rigid transformations is, always, a hard task.

The process of non-rigid shapes matching could be formulated as an optimization problem. It consists on finding the full correspondence mapping all point of one surface to their equivalent points on a second one with the minimal distortion.

The non-rigid correspondence methods can be classified on two major categories. We consider full and partial correspondence. The full or the complete one aims mapping the entire surfaces. While Partial matching is more complex, since it requires identifying optimal sub-parts or regions on shapes for which a right correspondence can be found.

According to the matched points resolution, the matching task can be also distinguished by dense or sparse. For the dense correspondence, the goal is to establish mapping between a large number of points on surface or even to find point-to-point matching of the corresponding shape. The sparse one is defined for a small number of elements or a set of feature points locally described. Besides, considering smaller sets of points permits to reduce the computational complexity comparing with the dense correspondence. Note that these categories of shapes matching may be combined, hence full as well as partial correspondence can be whether sparse or dense.

This challenging task has been widely addressed in the literature. Van Kaick et al. [20] proposed a detailed survey on 3D shapes matching methods.

For the full correspondence category, notably the sparse resolution, diverse methods have been proposed to tackle this problematic. The most common ones of these methods consist first on extracting a set of feature points and then constructing intrinsic surface descriptors.

Within this context, Ovsjanikov et al. [13] proposed an approach which relies on matching feature points in a space of a heat kernel for a given point on a surface and then the correspondence is obtained by searching the most similar heat kernel maps. Otherwise, Funkhouser et al. [17] used a Mobius voting approach. This approach consists on applying a Mobius transform for triplets of points, then generating conformal maps and voting for each couple of correspondences. The pairs with high votes are hence matched. On the other hand, Zhang et al. [21] proposed a method based on searching for correspondences while minimizing alignment and deformation error. Other alternative for the sparse correspondence searching, is presented by Bronstein et al. [2], by introducing the generalized multidimensional scaling (GMDS) which allows finding the minimum-distortion embedding of one surface into another. In the same context, Sahillio\(\breve{g}\)lu et al. [16] presented a method based on greedy optimization of an isometric distortion function.

For the dense resolution of full matching, various methods in the literature seek to find correspondences for all points on a surface. Kim et al. [8] proposed to combine multiple low-dimensional intrinsic maps to produce a blended map. They, then, associated confidence and consistency weights to each map and find the best blending to establish a final correspondence. For the same purpose, Bronstein et al. [3] replaced the geodesic distance in the GromovHausdorff framework by the diffusion distance. Jiang et al. [7] proposed to embed the shapes into a spectral domain, and, then, to find the correspondence using a non-rigid variant of the ICP (Iterative Closest Point) algorithm. Very recently, Lähner et al. [10] proposed an algorithm for non-rigid 2D-to-3D shape matching, which consists on finding the shortest circular path on the product 3-manifold of the surface and the curve. In [14], Sahillio\(\breve{g}\)lu et al. implemented an optimization mechanism of the method proposed by Sahillio\(\breve{g}\)lu et al. [16] idea improved within the EM (Expectation-Maximization) framework and coupled with a more sophisticated sampling scheme.

Furthermore, for the partial category, with either dense or sparse resolution, different approaches exist for searching the correspondent parts of two given shapes. The common strategy to establish partial correspondence, for some methods, is to represent regions by descriptors and, then, to search the similarity between them. In this context, Funkhouser et al. [5] introduced a priority-driven algorithm for searching similar shapes from a large database of 3D objects. The authors used a priority queue of potential sets of partial matches sorted by a cost function representing feature dissimilarity and geometric deformation. Van Kaick et al. [19] explored the bilateral map in order to present a local shape descriptor whose region of interest is defined by two feature points instead of one unlike classical descriptors. Van Kaick et al. proved that their approach is more effective for partial matching but it’s not the case when the shapes have strong intrinsic symmetries.

Hence some other alternatives compute a partial correspondence without relying on shape descriptors, such as the method of Bronstein et al. [1] whose main contribution is a framework for regularized partial matching of shapes taking into account three criteria, which are the regularity, the similarity and the size of parts. The authors relied on the Mumford-Shah functional to formulate the regularization criterion. Moreover, Sahillio\(\breve{g}\)lu et al. [15] proposed a rank-and-vote-and-combine (RAVAC) algorithm that identifies and ranks potentially correct matches by exploring the space of all possible partial maps between shape extremities.

We propose in this paper a novel method for the sparse correspondence between 3D shapes undergoing non-rigid transformations. Our approach consists on an intrinsic local description of surfaces extremeties based on the construction of local discrete representation known by Darcyan Coordinates System. For each discrete representation around an extreme point on the surfaces, principal curvatures field are computed as well. The most similar point descriptors are therefore matched in the mean of the \(L_{2}\) distance.

Thus, the present paper is organized as follow: We present in the second section the proposed descriptor construction steps. The next section is consecrated to the representation of our 3D local matching approach. Experimentations on 3D objects from the Tosca database and results discussion are illustrated in the fourth section.

2 3D Local Surface Description

We propose to locally describe the partial part of a surface around an extreme point of the surface with curvature. Since two acquisitions of the same object can lead to two different meshes, it is necessary first of all to extract a local representation around the extreme point which is invariant under the initial parametrization of the mesh. We propose to use the Darcyan representation. We present in this part all the steps of the proposed approach construction.

2.1 Brief Recall of the Darcyan Representation

Here, we recall the construction process of the well known Darcyan Coordinates System introduced by D’Arcy Thompson [18]. The parametric surface representation based on these coordinates system relatively to a given reference point on surface is, in fact, obtained by the superposition of the geodesic level curves around the reference point and the radial lines coming from the same point.

Let S be a two dimensional differential manifold, and let consider \(U_{r}\) the geodesic potential field coming from a reference point r on S. The function \(U_{r}: S\rightarrow R^{+}\) computes for any point p on S the length of the geodesic curve joining it to the reference point r. This function is well defined, since a geodesic curve between two points of a 2D differential manifold exists [4].

Construction of the Geodesic Level Curves. A geodesic level curve of value equal to \(\lambda \) around a reference point r on the surface S can be formulated as follows:

$$\begin{aligned} L_{r}^{\lambda }=\{p\in S; U_{r}(p)=\lambda \} \end{aligned}$$
(1)

\(L_{r}^{\lambda }\) is materialized by a set of all points on S having the same geodesic distance \(\lambda \) from r. Therefore, the surface S can be approximately reconstructed by all geodesic level curves, as shown in Fig. 1, so that, \(S \approx \cup _{\lambda } L_{r}^{\lambda }\).

Fig. 1.
figure 1

Geodesic level curves analytically extraction on a sphere around a reference point

Full size image

Construction of the Radial Lines Curves. We remind as well as the process of radial lines curves construction from a reference point r of the surface S.

Like mentioned in [6], the radial curves represent a solution of the following system:

$$\begin{aligned} \left\{ \begin{array}{l l l} \frac{dP(t)}{dt}= -\nabla U_{r}(P)\\ P(0)=r\\ \frac{dP(t)}{dt} \mid _{t=0}=\alpha \end{array} \right. \end{aligned}$$
(2)

where P(t) is the geodesic path emanating from r and following the opposite gradient \( \nabla \) direction on \(U_{r}\). Radial lines curves, denoted by \(C^{\alpha }\), are therefore generated according to the angular direction \(\alpha \) which can be arbitrary taken. Similar to geodesic level curves, a reconstruction of the surface S can be approximated by \(\cup _{\alpha } C^{\alpha }\), which is illustrated in Fig. 2.

Fig. 2.
figure 2

Radial lines analytically extraction on a sphere from a reference point

Full size image

The Darcyan Representation. Here we define Darcyan representation D as the superposition of both n geodesic level curves and m radials lines curves relatively to a given point r.

$$\begin{aligned} D_{k,l}(r)=\left\{ L_{r}^{\lambda k}\bigcup C_{r}^{\alpha l}, 1\le k \le n , 1\le l \le m \right\} \end{aligned}$$

The Fig. 3 illustrates the steps of Darcyan coordinate system construction. We propose to compute the curvature values on the intersection points between the radial line curves and the geodesic level ones of the Darcyan coordinate system. These intersection points are invariant under the M(3) group. They are also ordered since each point on the surface is indexed by the geodesic level curve and the radial one to which it belongs.

Fig. 3.
figure 3

Darcyan system reconstruction: Radial curves (a), geodesic level curves (b) and the superposition of both system of curves (c)

Full size image

2.2 Brief Recall of the Curvature Computation

First, we recall some useful facts for the curvature calculation on surface.

Let \(X:(u,v) \in D\subset R^{2}\rightarrow (x(u,v),y(u,v),z(u,v)) \in S \subset R^{3} \) be a parametrization of S.

We consider \(\left( u,v\right) \) the corresponding basis of the tangent plane to S at a point \(p=X(x_{u},x_{v})\). \(N(p)=\frac{x_{u} \wedge x_{v}}{\left\| x_{u} \wedge x_{v}\right\| }\) is the normal vector to S at p, according to a chosen orientation.

Therefore, the curvature expression is given using the following coefficients:

$$E=x_{u}.x_{u}, F=x_{u}.x_{v}, G=x_{v}.x_{v}, L=x_{uu}.{{\mathop {N}\limits ^{\rightarrow }}}, M=x_{uv}.\mathop {N}\limits ^{\rightarrow } and N=x_{vv}\mathop {N}\limits ^{\rightarrow } $$
$$\begin{aligned} K_{G}=\frac{LN-M^{2}}{EG-F^{2}} \quad \quad \\ K_{M}=\frac{EN-2FM+GL}{2(EG-F^{2})} \end{aligned}$$

where E, F and G are the first fundamental coefficients, while L, M, N are the second fundamental coefficients.

The quantities \(K_{G}\) and \(K_{M}\) are the Gaussian curvature and the Mean curvature \(p=X(x_{u},x_{v})\) respectively.

Thus, the principal curvatures are derived from these expressions \(K_{G}=k_{max}.K_{min}\) and \(K_{M}=\frac{(k_{max}+k_{min})}{2}\).

\(k_{max}\) and \(k_{min}\) define the principal curvatures of the surface as, respectively, the maximal and the minimal curvature.

2.3 The Proposed Local Descriptor

We propose, here, a novel 3D shape descriptor which explores an intrinsic geometric property, principal curvatures fields on a local parametrization which is invariant under Euclidean motions.

The proposed descriptor relies on the computation of both principal maximal and minimal curvature field for a discrete point picked on a given geodesic level curve in the parametric representation, since these set of curves are invariant under motion group.

We propose to compute the mean of both principal maximal and minimal curvatures on the intersection points of each geodesic level curve. We denote \(\overline{k_{max}^{i}}\) and \(\overline{k_{min}^{i}}\) the mean of, respectively, \(k_{max}\) and \(k_{min}\) for the \( i^{th} \) geodesic level curve. Thus the novel descriptor is defined by \(\left\{ \overline{k^{1}_{max}}, \overline{k^{1}_{min}},..,\right. \left. \overline{k^{i}_{max}}, \overline{k^{i}_{min}},..., \overline{k^{m}_{max}}, \overline{k^{m}_{min}}\right\} _{1\le i \le m}\). The proposed descriptor is illustrated in Fig. 4.

Fig. 4.
figure 4

Illustration of the proposed descriptor: (a) the Darcyan representation construction, (b) the vector of curvature fields computation and the obtained intersection points (in red) (Color figure online)

Full size image

3 3D Local Matching

Once we have extracted a set of interest points from shape extremities by appling an Average Geodesic Distance function (AGD) [12], we compute the proposed descriptor around each point of interest. Thus, we obtain vectors consisting on the mean of principal curvature field values.

Fig. 5.
figure 5

3D local correspondence approach

Full size image

After executing this process, the similarity between the acquired vectors is measured in term of \(L_{2}\) distance. The full approach that we propose in shown in Fig. 5. We compute, indeed, the distance between all the pairs of vectors and thereafter, correspondent points are found when the resulting distance is the minimal one.

4 Experimentation

In order to evaluate the efficiency of our approach, we have conducted experiments on 3D shapes in different poses from the TOSCA database with non rigid deformations.

4.1 The Used Database

We have conducted experiments on the high resolution TOSCA database, widely used in a variety of 3D shape correspondence approaches, which contains 80 meshes representing people and animals in a variety of poses. The meshes are grouped in 8 groups with common topology [1]. The reference points and the correspondence ground truth are provided in the evaluation benchmark proposed by [8].

4.2 Approximation on 3D Meshes

Since 3D surfaces are usually approximated by 3D meshes, we use the Fast Marching Method [9] for the geodesic distances computation. We, as well as, use Meyer et al. [11] algorithm for estimating the principal curvatures. Moreover, All the objects of the database need a scale normalization.

Fig. 6.
figure 6

Local parametric representation around points of interest on David model and the obtained matches

Full size image

After being constructed, the sets of two curves require an interpolation step in order to increase their resolution and have more accuracy. Figure 6 illustrates local parametric representation around points of interest on the object David and the obtained matches.

Fig. 7.
figure 7

Matching results for various models from the Tosca database

Full size image
Fig. 8.
figure 8

Obtained correspondence rates for Tosca database

Full size image
Fig. 9.
figure 9

Incorrect matches (in blue) due to symmetric shapes (Color figure online)

Full size image

4.3 Obtained Results

We evaluate our matching approach on different objects of the Tosca database. The Fig. 7 illustrates the matching results obtained after applying our process for finding correspondence between some pairs of shapes in various poses. In fact, our method involves important correspondence rates varying from \(75\%\) to \(100\%\) according the used objects. The Fig. 8 shows the resulting rates.

Beside the correct matches, correspondence errors exist in some cases. One of the factors of these errors, is the symmetry which appears in the majority of the used database models. Figure 9 shows some matching errors which are due to symmetric objects.

5 Conclusion

We have presented a novel approach to find correspondence between pairs of surfaces undergoing non-rigid transformations. This procedure relies on the mean values of principal curvature fields computation on a intrinsic local parametrization around reference points extracted on the extremities of the surface.

The performance of our local descriptor is evaluated on different models from Tosca Database. The first obtained results have shown the accuracy of our approach to establish correspondence among non-rigid shapes relying on the reached rates. Nevertheless, our matching process may lead to some errors and this is due to certain limitations such as symmetric shapes.

We intend, in the future works, to handle the problem of symmetry. In other hand, we aim to optimize the resolution of our intrinsic descriptor by finding the optimal number of the geodesic levels curves and the radial lines curves. We also intend to perform the experimentation on others 3D non-rigid databases with different properties.