## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Generalized multidimensional scaling: A framework for isometry-invariant partial surface matching

Edited by Alexandre J. Chorin, University of California, Berkeley, CA, and approved December 5, 2005 (received for review October 3, 2005)

## Abstract

An efficient algorithm for isometry-invariant matching of surfaces is presented. The key idea is computing the minimum-distortion mapping between two surfaces. For this purpose, we introduce the generalized multidimensional scaling, a computationally efficient continuous optimization algorithm for finding the least distortion embedding of one surface into another. The generalized multidimensional scaling algorithm allows for both full and partial surface matching. As an example, it is applied to the problem of expression-invariant three-dimensional face recognition.

The problem of comparison between deformable surfaces has recently gained a lot of interest in pattern recognition, especially in the fields of three-dimensional (3D) data analysis and synthesis. It arises, for example, in medical imaging and 3D face recognition (1). The fundamental question is how to efficiently yet accurately quantify the similarity between a given reference surface (“model”) and some other surface (“probe”), usually a bent version of the model. We would like to capture the distinction of the model and the probe intrinsic properties associated with the metric structure of the surface, while ignoring the extrinsic properties that describe the way the surface is immersed into the ambient space and that often change while the surface bends. A deformation that preserves the intrinsic structure of the surface is called an “isometry.” Our goal is thus to define a computable isometry-invariant measure of similarity between surfaces. Often, the problem is even more complicated, when having the probe only partially available (see Fig. 1). We refer to the latter setting as “partial surface matching” and to the case where the whole probe is available as “full surface matching.”

One of the earliest attempts of isometry-invariant surface matching is the classical “iterative closest point” (ICP) algorithm (2). It addresses a particular case of the partial matching problem where only rigid (Euclidean) isometries are allowed. An efficient method for the construction of near-isometry-invariant representations of surfaces [called “canonical forms” (CFs)] based on Euclidean embeddings was presented in ref. 3, as a generalization of ref. 4. This approach used a multidimensional scaling (MDS) algorithm (5). MDS is closely related to dimensionality reduction (6, 7) and can be performed in a computationally efficient manner. Euclidean embeddings are used in theoretical computer science for representing metric spaces usually arising from geometry of graphs (8).

Recently, Mémoli and Sapiro (9) used a discrete probabilistic approximation of the Gromov–Hausdorff (GH) distance (10) to compare surfaces up to isometries. Their method is based on evaluating all of the permutations between two sets of surface samples, which is computationally expensive. Moreover, the GH distance is inappropriate for partial matching.

One of the main contributions of this work is the “partial embedding (PE) distance,” motivated by Gromov's theory. For its efficient computation, we introduce an extension of MDS, which we call generalized MDS (GMDS). The key idea is to measure the minimum possible distortion when trying to isometrically embed one surface into another. By using GMDS, we can handle full as well as partial surface matching; this ability is one advantage over a straightforward use of the GH distance.

The paper includes six sections. *Theoretical Foundations* reviews the theoretical background of isometry-invariant surface matching. We define our PE distance and show its relation to the GH distance. The next section, *GMDS*, introduces the GMDS and addresses implementation and optimization considerations. The following section focuses on partial matching of surfaces. In *Results*, we show some experimental results. Finally, in *Conclusions*, we summarize the main results.

## Theoretical Foundations

Let 𝒮 be a smooth connected and compact Riemannian surface (2-manifold) with minimal geodesics $$mathtex$$$$mathtex$$ between every two points *s, s*′ on 𝒮. We define the geodesic distance function $$mathtex$$$$mathtex$$ as $$mathtex$$$$mathtex$$. For simplicity, we say “surface 𝒮,” implying the underlying metric space (𝒮, $$mathtex$$$$mathtex$$) induced by the Riemannian geometry of 𝒮. In our context, surfaces are perceived from a point of view of metric rather than Riemannian geometry; in other words, we do not consider the Riemannian metric tensor, but the geodesic distances that it induces.

In practice, we work with sampled surfaces represented by finite metric spaces. We use the following definitions: a subset $$mathtex$$$$mathtex$$ is called an “*r*-covering” (or “*r*-net”) of 𝒮 if $$mathtex$$$$mathtex$$, where $$mathtex$$$$mathtex$$ is a ball of radius *r* around $$mathtex$$$$mathtex$$ in 𝒮. A finite *r*-covering of 𝒮 consisting of *N* points is denoted by $$mathtex$$$$mathtex$$. The metric on $$mathtex$$$$mathtex$$ is assumed to be the restricted metric $$mathtex$$$$mathtex$$ for all *s, s*′ in $$mathtex$$$$mathtex$$. We refer to *r* as the “sampling radius.” An arbitrary finite sampling of 𝒮 consisting of *N* points is denoted by $$mathtex$$$$mathtex$$, where 𝒮 is called “continuous” surface and $$mathtex$$$$mathtex$$ a “discrete” one.

Given two surfaces 𝒮 and 𝒬, a transformation $$mathtex$$$$mathtex$$ is said to have “distortion” ε if $$mathtex$$$$mathtex$$[1] If in addition ψ is “ε-surjective” (i.e., $$mathtex$$$$mathtex$$ for all $$mathtex$$$$mathtex$$), it is called an “ε-isometry,” and 𝒮 and 𝒬 are called “ε-isometric.” A transformation $$mathtex$$$$mathtex$$ with disψ = 0 is called an “isometry,” and 𝒮 and 𝒬 admitting such a transformation are called “isometric.”

**Desiderata.** Let $$mathtex$$$$mathtex$$ denote the set of all compact surfaces, and $$mathtex$$$$mathtex$$ denote the space of the equivalence classes of all isometries in $$mathtex$$$$mathtex$$, i.e., a space in which a point represents all of the isometries of a compact surface. The similarity between two surfaces can be measured by a nonnegative function $$mathtex$$$$mathtex$$. Slightly abusing common terminology, we call *d* a “distance,” even if it is not a metric (often, the terms distance and metric are used as synonyms). Ideally, *d* should have the following properties: **Desideratum 1: Metric properties.** *d* should be a metric on $$mathtex$$$$mathtex$$, i.e., satisfy the following conditions: (*i*) $$mathtex$$$$mathtex$$; (*ii*) $$mathtex$$$$mathtex$$ if and only if 𝒮 and 𝒬 are isometric; (*iii*) $$mathtex$$$$mathtex$$; and (*iv*) $$mathtex$$$$mathtex$$.

** Desideratum 2: Good similarity measure.** (

*i*) If $$mathtex$$$$mathtex$$, then 𝒮 and 𝒬 are

*c*ε-isometric, where

*c*> 0; (

*ii*) if 𝒮 and 𝒬 are ε-isometric, then $$mathtex$$$$mathtex$$. Note that this property is an extension of the second metric axiom to ε-isometries.

** Desideratum 3: Partial matching.** Given $$mathtex$$$$mathtex$$ a patch of $$mathtex$$$$mathtex$$.

** Desideratum 4: Sampling consistency.** Given $$mathtex$$$$mathtex$$ a finite

*r*-covering of 𝒮, lim

_{r→0}$$mathtex$$$$mathtex$$.

*Desideratum 5: Efficiency*.*d* should be efficiently computable, or at least efficiently approximated. Practically, given $$mathtex$$$$mathtex$$ and $$mathtex$$$$mathtex$$, $$mathtex$$$$mathtex$$ should be computable in polynomial time of max{*N, N*′}.

Unfortunately, *Desiderata 1*–*3* cannot be satisfied simultaneously. In practical applications, it appears that the ability to perform partial matching (*Desideratum 3*) is the most important. In the following, we show a distance that satisfies *Desiderata 3*–*5* and partially satisfies *Desiderata 1* and *2*.

It is also important to note the following potential danger of partial matching. Consider an extreme case of $$mathtex$$$$mathtex$$ a smooth curve on 𝒮. Obviously, if *Desideratum 3* holds, then $$mathtex$$$$mathtex$$, but because a curve can be isometrically embedded into any surface, we are liable to obtain $$mathtex$$$$mathtex$$ for an arbitrary 𝒬 which is not necessarily an isometry of 𝒮. In other words, when the patch $$mathtex$$$$mathtex$$ is too small, partial matching can be meaningless.

**GH Distance.** Our starting point is the GH distance introduced by Gromov (10). For two surfaces 𝒮 and 𝒬, it is defined by $$mathtex$$$$mathtex$$[2] where $$mathtex$$$$mathtex$$ and $$mathtex$$$$mathtex$$ are isometric embeddings into the metric space 𝒵, and $$mathtex$$$$mathtex$$ is the Hausdorff distance between two subsets 𝒜 and ℬ of 𝒵. For bounded metric spaces 𝒮 and 𝒬, the GH distance can be defined alternatively (11) as $$mathtex$$$$mathtex$$[3] where $$mathtex$$$$mathtex$$.

**Proposition 1.** (*i*) *d*_{GH} *is a finite metric on* $$mathtex$$$$mathtex$$; (*iia*) $$mathtex$$$$mathtex$$, *then* 𝒮 *and* 𝒬 *are 2*ε*-isometric*; (*iib*) *if* 𝒮 *and* 𝒬 *are ε-isometric, then* $$mathtex$$$$mathtex$$; (*iii*) $$mathtex$$$$mathtex$$.

For proof, see ref. 11. The GH distance is a good candidate for our “ideal” similarity measure between surfaces; from *Proposition 1*, it follows that it satisfies *Desiderata 1, 2, and 4*. However, there are two main disadvantages. First, the definition in Eq. **2** is hard to compute (*Desideratum 5*). Secondly, *d*_{GH} does not allow partial matching between surfaces (*Desideratum 3*). Given a surface 𝒮 and a patch $$mathtex$$$$mathtex$$ such that $$mathtex$$$$mathtex$$ (i.e., $$mathtex$$$$mathtex$$ is an *R*-covering of 𝒮), we can have $$mathtex$$$$mathtex$$, although we would have liked the distance between the surface and its patch to be zero.

**PE Distance.** Our assumption was that we have a model surface 𝒮 and a probe surface 𝒬, which can be an isometrically bent patch of 𝒮. We have shown the disadvantage of the GH distance in such a setup. To allow for partial matching, we define the PE distance $$mathtex$$$$mathtex$$ Intuitively, *d*_{PE} measures the metric distortion obtained trying to embed 𝒬 into 𝒮 in the “most isometric” way. It can be considered as only part of the GH distance, hence the name “partial embedding.”

**Proposition 2.** *Given compact surfaces *𝒮*, *𝒬*, and *ℛ*, d*_{PE} *satisfies the following*: (*i*) $$mathtex$$$$mathtex$$; (*ii*) $$mathtex$$$$mathtex$$ *if and only if* 𝒬 *is isometrically embeddable into* 𝒮; *and* (*iii*) $$mathtex$$$$mathtex$$.

The proof closely follows the proof of theorem 7.3.30 in ref. 11. Still, *d _{PE}* is not symmetric, and furthermore, not a metric on $$mathtex$$$$mathtex$$. Yet, it satisfies some important metric properties, namely: isometry invariance (

*ii*) and a nonsymmetric version of the triangle inequality (

*iii*). Moreover, because $$mathtex$$$$mathtex$$, it trivially follows that if 𝒮 and 𝒬 are ε-isometric, then $$mathtex$$$$mathtex$$. The converse appears to be true under some more restrictive assumptions.

Given $$mathtex$$$$mathtex$$ an *r*′-covering of 𝒬 and $$mathtex$$$$mathtex$$ an *r*-covering of 𝒮, the distance *d*_{PE} satisfies the following important properties.

**Proposition 3.** (*i*) $$mathtex$$$$mathtex$$; (*ii*) $$mathtex$$$$mathtex$$.

Note that a particular case of *Proposition 3i* is $$mathtex$$$$mathtex$$. By virtue of this property [unlike the corresponding property $$mathtex$$$$mathtex$$ of the GH distance] *d*_{PE} can be used for partial matching. From *Proposition 3* it also follows that *d*_{PE} is consistent with sampling, namely, given $$mathtex$$$$mathtex$$ and $$mathtex$$$$mathtex$$, two finite coverings of 𝒮 and 𝒬, respectively, $$mathtex$$$$mathtex$$. In other words, if we have a sufficiently dense finite sampling of 𝒮 and 𝒬, $$mathtex$$$$mathtex$$ is a sufficiently good approximation of the continuous $$mathtex$$$$mathtex$$.

## GMDS

We now address a practical question of computing the PE distance *d*_{PE}, starting from the CF approach (3). As an input, we assume to be given two sampled surfaces $$mathtex$$$$mathtex$$ and $$mathtex$$$$mathtex$$ (representing the continuous surfaces 𝒮 and 𝒬) and the geodesic distances between the samples $$mathtex$$$$mathtex$$, represented by an *N* × *N* matrix $$mathtex$$$$mathtex$$ and an *N*′ × *N*′ matrix $$mathtex$$$$mathtex$$, respectively.

**CF and MDS.** In ref. 3, a three-stage algorithm for near-isometry-invariant surface matching was proposed. First, each of the sampled surfaces is mapped into an *m*-dimensional Euclidean space by a near-isometric embedding (“flattening”), obtained by minimization of the stress function $$mathtex$$$$mathtex$$[4] with respect to **X**. Here **X** denotes an *N* × *m* matrix of coordinates in $$mathtex$$$$mathtex$$ and $$mathtex$$$$mathtex$$ denotes the Euclidean metric. The stress can be thought of as an *L*_{2} measure of the distance distortion caused by such an embedding. Algorithms minimizing the stress σ are known as MDS (5). $$mathtex$$$$mathtex$$ was called the CF of 𝒮 and $$mathtex$$$$mathtex$$ the “embedding error.”

Note that both CFs $$mathtex$$$$mathtex$$ and $$mathtex$$$$mathtex$$ are defined up to a Euclidean isometry (translation, rotation, and reflection). The second stage of the CF algorithm is an alignment of the CFs in $$mathtex$$$$mathtex$$. It is obtained by diagonalizing the matrix of the second-order moments of each CF and reordering the axes such that the diagonal values are decreasing. At the third stage, a distance based on high-dimensional moments (12) is used to compare between $$mathtex$$$$mathtex$$ and $$mathtex$$$$mathtex$$ after alignment. We denote the distance computed as rigid matching between the CFs of 𝒮 and 𝒬 by $$mathtex$$$$mathtex$$.

Strictly speaking, the CF algorithm satisfies only *Desideratum 5*: it is very efficient computationally. Unlike alternative combinatorial optimization, MDS is a continuous optimization problem that can be solved by using standard optimization methods. Recently, a very fast multigrid implementation of the CF algorithm was presented (13). However, the CF algorithm lacks true isometry invariance, because, in general, a surface cannot be isometrically embedded into $$mathtex$$$$mathtex$$ (8). Therefore, this “flattening” introduces a distortion. Partial matching is also complicated using this method.

**GMDS.** The choice of a Euclidean embedding space was followed by other choices, like spherical (14, 15) and hyperbolic (16) spaces, which proved to be advantageous for certain types of surfaces. The leitmotif of these approaches was choosing an embedding space in which geodesic distances can be expressed analytically. In practice, this possibility exists only in very few spaces.

So far, we have a good distance *d*_{PE}, which satisfies our Desiderata on one hand and an efficient MDS framework for computing Euclidean embeddings on the other. Next, we unite theory and practice. We consider a generalization of the CF algorithm, where surface 𝒬 is embedded into another surface 𝒮. We call this procedure GMDS.

Similarly to the Euclidean case, we define the “generalized stress” as $$mathtex$$$$mathtex$$[5] for 1 ≤ *p* < ∞, and $$mathtex$$$$mathtex$$[6] for *p* = ∞, where the matrix **U** represents the positions of *N*′ points on 𝒮 in some local or global parametric coordinates **u*** _{i}*, and

**W**= (

*w*) is a symmetric matrix of nonnegative weights. Note two major differences between the generalized stress and the stress defined in Eq.

_{ij}**4**. First, the metric in the embedding space is not given analytically and is approximated numerically. Second, we replaced the

*L*

_{2}norm with a general

*L*norm. Furthermore, the embedding space is 2D, whereas in the CF algorithm the embedding space is usually at least 3D.

_{p}Choosing *p* = ∞ and **W** = **1**_{N}_{′×}_{N}_{′} we easily obtain $$mathtex$$$$mathtex$$[7] In other words, minimization of the generalized stress allows computing *d*_{PE}. We discuss numerical optimization algorithms in the next section. Moreover, it is possible to use a similar procedure to compute *d*_{GH} (A.M.B., M.M.B., and R.K., unpublished results). In practice, it may appear that using other norms (e.g., *p* = 1 or 2) is advantageous. In *Results*, it is demonstrated that we can use σ* _{p}* instead of σ

_{∞}. We denote the distances computed in this way by $$mathtex$$$$mathtex$$.

An important issue is the choice of the parametric coordinate system. For objects acquired by a range sensor, a global parameterization $$mathtex$$$$mathtex$$ (e.g., *U* = [0, 1]^{2}), is usually readily available. Global parameterization for objects with simple topology (e.g., homeomorphic to a disk or a sphere) is also very natural and can be computed from a general mesh by a variety of efficient algorithms (17). For general triangulated meshes with more sophisticated topology, a single global coordinate system may not suffice. In this case, we propose to resort to a local coordinate system, in which each point on the polyhedral surface is represented by an index of the triangle to which it belongs and a vector in barycentric coordinates of that triangle.

**Computational Aspects.** For 1 ≤ *p* < ∞, the generalized stress (5) is minimized by iterative gradient-type optimization methods. Here, for simplicity, we show the projected gradient descent algorithm for a global parameterization $$mathtex$$$$mathtex$$; in practice, more efficient optimization methods such as conjugate gradients or quasi-Newton can be used (18). Note that because the problem is nonconvex, such methods do not necessarily converge to the global minimum. Local convergence is a common property of MDS algorithms in general (5). There are standard techniques to overcome this problem, one of which is multiscale or multigrid optimization (13, 19).

For globally parameterized surfaces, the basic projected gradient descent step has the form $$mathtex$$$$mathtex$$[8] where *k* is the iteration number. α^{(}^{k}^{)} is the step size on iteration *k*; it is computed by using line search (18). *P _{U}* is a projection operator onto the parameterization domain

*U*. ∇

*σ*

_{U}*denotes the gradient of σ*

_{p}*with respect to*

_{p}*U*; its computation involves the distance function $$mathtex$$$$mathtex$$ and its derivatives. When local parameterization is used, the optimization algorithm becomes more elaborate because it requires unfolding of triangular faces of the mesh and transformations from barycentric coordinate systems of adjacent triangles to compute the stress along a polylinear path.

For *p* =∞, we rewrite the problem as constrained optimization with 2*N*′ + 1 variables, a linear objective, and *N*′(*N*′ – 1) nonlinear inequality constraints using an artificial variable τ, $$mathtex$$$$mathtex$$[9] for all *j* > *i*. The new problem can be solved by using standard constrained optimization methods, e.g., penalty-barrier or augmented Lagrangian (18).

A crucial part of the GMDS problem is the computation of geodesic distances. We use the “fast marching method” (FMM), which computes geodesic distances on surfaces by solving the eikonal equation on general triangulated meshes (20). For surfaces with global parameterization, an efficient parametric version of FMM can be used (21).

The geodesic distances $$mathtex$$$$mathtex$$ between the fixed sample points of 𝒬 can be precomputed by using FMM. However, the **u*** _{i}* coordinates representing ψ(

*q*) on 𝒮 change during the iterations of the numerical minimization algorithm. Thus, the distances $$mathtex$$$$mathtex$$ have to be reevaluated at every iteration. This computation is critical for the GMDS. For that goal, we developed a numerical procedure we termed “three-point geodesic distance approximation” (detailed in A.M.B., M.M.B., and R.K., unpublished results). The idea is to produce a $$mathtex$$$$mathtex$$ for $$mathtex$$$$mathtex$$ and its derivatives, interpolating their values from the matrix $$mathtex$$$$mathtex$$.

_{i}## Partial Matching of Surfaces

Let $$mathtex$$$$mathtex$$ be a patch of the probe surface 𝒬 and $$mathtex$$$$mathtex$$ be a finite sampling of $$mathtex$$$$mathtex$$. We consider the distance $$mathtex$$$$mathtex$$. By our assumption, the metric on $$mathtex$$$$mathtex$$ is given by restricting $$mathtex$$$$mathtex$$ to $$mathtex$$$$mathtex$$, i.e., $$mathtex$$$$mathtex$$. However, $$mathtex$$$$mathtex$$ is computed numerically on $$mathtex$$$$mathtex$$ and can be inconsistent with $$mathtex$$$$mathtex$$. Specifically, geodesics $$mathtex$$$$mathtex$$ that pass through points not in $$mathtex$$$$mathtex$$ give rise to inconsistent distances, for which $$mathtex$$$$mathtex$$ (see Fig. 2 *Left*). In pathological cases, the inconsistency can be arbitrarily large, such that the partial matching makes little sense (see example in Fig. 2 *Right*). To overcome this problem, the inconsistent distances $$mathtex$$$$mathtex$$ are excluded defining the corresponding weights *w _{i}*

_{′}

_{j}_{′}= 0 in the generalized stress (5).

A practical problem is how to locate the inconsistent distances. Here, we address two possible scenarios: (*i*) the original surface 𝒬 or its sampled version is available, and (*ii*) only the patch $$mathtex$$$$mathtex$$ is available. The first scenario occurs when one intentionally wishes to exclude parts of 𝒮. In this case, the FMM is used first to compute the distances $$mathtex$$$$mathtex$$ on 𝒬 and then to compute the distances $$mathtex$$$$mathtex$$ on $$mathtex$$$$mathtex$$. The weights are defined by $$mathtex$$$$mathtex$$[10] for all *i, j* = 1,... *N*′, where δ stands for the FMM accuracy.

The second scenario is more general, yet, also a more challenging one. When the original surface 𝒬 is not given, we need to rely on the geodesics of $$mathtex$$$$mathtex$$ [in practice, the geodesics are found by backtracking (22)]. By using the above argument, we define the weights by $$mathtex$$$$mathtex$$[11] for all *i, j* = 1,... *N*′, i.e., excluding the geodesics that touch the boundary $$mathtex$$$$mathtex$$.

## Results

**Embedding of Spherical Surfaces.** In the first experiment, we measured the distance between a unit 2D sphere sampled at 3,200 points, and spheres with radii in the range 0.5 ÷ 2.5 sampled at a smaller number of points according to the “farthest point sampling” strategy (23) with random seed. Ten random samplings were used for each radius. Two distance measures were compared as follows: $$mathtex$$$$mathtex$$ (obtained by the minimization of σ_{2}) with 100 points and *d*_{CF} with 100, 250, and 500 points. $$mathtex$$$$mathtex$$ was computed by using local representation of points in barycentric coordinates.

Fig. 3 presents the normalized distances as a function of the sphere radius. The PE distance appears to be extremely sensitive to the geometry of the surfaces. A change as small as 0.1% in the spherical patch radius (from 1.000 to 0.999) increases $$mathtex$$$$mathtex$$ by a value exceeding the variance due to sampling. Similar results are achieved by using $$mathtex$$$$mathtex$$ and slightly inferior with *d*_{PE}. For comparison, with the same number of points (100), *d*_{CF} is unable to discern between spheres differing in radius by >10%. Increasing the number of points to 500 makes *d*_{CF} sensitive to radius differences of ≈2.64%, which is still one order of magnitude below the sensitivity of $$mathtex$$$$mathtex$$.

**Face Recognition Experiment.** In ref. 1, we argued that human facial expressions can be approximated as isometric deformations of the face. According to this model, an isometry-invariant surface matching results in nearly expression-invariant face recognition.

In the second experiment, we applied our GMDS approach to a toy 3D face recognition problem of 20 faces taken from the Notre Dame database (24, 25). This experiment is a proof of concept rather than a real face recognition system, in which tune-up plays a crucial role. As models, four different faces were used. Each model contained ≈4,000 points and was represented using global parameterization. As probes, faces of the same subjects with different facial expressions were used (four probes per subject). The probes were cropped by using geodesic mask (1) (leaving the nose and the forehead region) and subsampled by using farthest point sampling to 53 points. The distances were weighted as described in *Partial Matching of Surfaces*. All distances were computed numerically by a version of parametric FMM (21). We tested two PE distances: *d*_{PE} and $$mathtex$$$$mathtex$$. The generalized stress function σ* _{p}* and its gradient were implemented in C. For 53 points, the computation time of σ

*and ∇σ*

_{p}*is ≈20 msec on a Pentium IV PC.*

_{p}Fig. 4 presents the similarity matrix using *d*_{PE} (second row) and $$mathtex$$$$mathtex$$ (third row). Note that in both cases a perfect separation between different subjects is possible; the recognition error (measured as “equal error rate,” or EER; see ref. 26 for definition) in this experiment is zero. It is worthwhile noting the small number of points used in our approach, compared with the thousands of points required in the CF algorithm (1, 3) for similar matching accuracy.

An example of embedding probe surfaces of two different subjects into a model surface are shown in Fig. 5. The figure depicts the “local stress” value $$mathtex$$$$mathtex$$ at each point *k*, i.e., the contribution of the *k*th point to σ_{∞}. We observe that the stress is significantly higher when a probe of a different subject is embedded into the model surface.

## Conclusions

We introduced a previously undescribed framework for isometry-invariant surface matching. The core of our computational framework is GMDS, a generalization of standard MDS algorithms. We exemplified the use of GMDS on the problem of expression-invariant 3D face recognition.

Our approach favorably compares with previous attempts to perform isometry-invariant surface matching. First, our PE distance naturally allows for isometry-invariant matching of partially missing surfaces. Secondly, the properties of our distance and its computation are completely deterministic. Thirdly, GMDS used for the PE distance computation is a continuous optimization problem and can be solved very efficiently by using standard optimization methods. In light of recent results (13), we expect a significant performance boosting by using multigrid methods.

Finally, we note that the same numerical framework can be applied for computing the exact GH distance between surfaces (A.M.B., M.M.B., and R.K., unpublished results).

## Acknowledgments

We thank Alexander Brook and Roman Goldenberg for thorough reading of the manuscript, Facundo Mémoli and Vladimir Lin for fruitful discussions, and our reviewers for valuable comments. This work was supported in part by United States-Israel Science Foundation Grant 2004274, Israeli Science Foundation Grant 738/04, and the Charles Krown Research Fund.

## Footnotes

↵† To whom correspondence should be addressed. E-mail: ron{at}cs.technion.ac.il.

Author contributions: A.M.B., M.M.B., and R.K. designed research, performed research, contributed new reagents/analytic tools, analyzed data, and wrote the paper.

Conflict of interest statement: No conflicts declared.

This paper was submitted directly (Track II) to the PNAS office.

Abbreviations: CF, canonical form; FMM, fast marching method; GH, Gromov–Hausdorff; MDS, multidimensional scaling; GMDS, generalized MDS; PE, partial embedding.

- Received October 3, 2005.

- Copyright © 2006, The National Academy of Sciences

## References

- ↵
- ↵
- ↵
- ↵
- ↵Borg, I. & Groenen, P. (1997) Modern Multidimensional Scaling: Theory and Applications (Springer, Berlin)..
- ↵Roweis, S. T. & Saul, L .K. (2000) Science 290.
**,**2323–2326.pmid:11125150 - ↵Donoho, D. & Grimes, C. (2003) Proc. Natl. Acad. Sci. USA 100.
**,**5591–5596. - ↵
- ↵
- ↵Gromov, M. (1981) Structures Métriques pour les Variétés Riemanniennes, Textes Mathématiques (CEDIC, Paris), Vol. 1..
- ↵Burago, D., Burago, Y. & Ivanov, S. (2001) A Course in Metric Geometry, Graduate Studies in Mathematics (Am. Math. Soc., Providence, RI), Vol. 33..
- ↵Tal, A., Elad, M. & Ar, S. (2001) in Proceedings of the 6th Eurographics Workshop on Multimedia, eds. Jorge, J. A., Correia, N. M., Jones, H. & Kamegai, M. B. (Springer, New York), pp. 97–108..
- ↵Bronstein, M. M., Bronstein, A. M., Kimmel, R. & Yavneh, I. (2006) Numerical Linear Algebra Applications, in press..
- ↵Elad, A. & Kimmel, R. (2002) Geometric Methods in Bio-medical Image Processing (Springer, Berlin), Vol. 2191, pp. 77–89..
- ↵
- ↵Walter, J. & Ritter, H. (2002) in Proceedings of the Eighth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, eds. Hand, D., Keim, D. & Ng, R. (Assoc. Computing Machinery, New York), pp. 123–131..
- ↵Floater, M. S. & Hormann, K. (2004) in Advances in Multiresolution for Geometric Modelling, eds. Dodgson, N. A., Floater, M. S. & Sabin, M. A. (Springer, Berlin), pp. 157–186..
- ↵Bertsekas, D. (1999) Nonlinear Programming (Atlanta Scientific, Atlanta), 2nd Ed..
- ↵
- ↵Kimmel, R. & Sethian, J. A. (1998) Proc. Natl. Acad. Sci. USA 95.
**,**8431–8435.pmid:9671694 - ↵Spira, A. & Kimmel, R. (2004) Interfaces Free Boundaries 6.
**,**315–327. - ↵Kimmel, R. (2003) Numerical Geometry of Images (Springer, Berlin)..
- ↵
- ↵
- ↵Flynn, P. J., Bowyer, K. W. & Phillips, P. J. (2003) in Audio- and Video-Based Biometric Person Authentication, eds. Kittler, J. & Nixon, M. S. (Springer, Berlin), pp. 44–51..
- ↵Ashbourn, J. (2002) Biometrics: Advanced Identity Verification (Springer, Berlin)..

## Citation Manager Formats

### More Articles of This Classification

### Physical Sciences

### Related Content

- No related articles found.

### Cited by...

- From understanding of color perception to dynamical systems by manifold learning
- Minimum action principle and shape dynamics
- Landmark-free geometric methods in biological shape analysis
- On convex relaxation of graph isomorphism
- Algorithms to automatically quantify the geometric similarity of anatomical surfaces