1 Introduction
Modern computer graphics has experienced a paradigm shift: Traditional manual modeling is increasingly complemented by datadriven techniques where measured data, such as 3D scans, are used as a basis for building 3D models. An important data source are 3D scans of deformable models, such as humans or animals in varying poses. Today, there exist sophisticated scanning setups for acquiring moving geometry in realtime [1, 2, 3, 4] and there are even consumer devices on the market such as the Microsoft Kinect. This leads to new applications such as virtual cinematography [5], or the creation of datadriven generative shape models of deformable objects [6, 7, 8, 9]. Finding correspondences among such data is a fundamental problem for all of these applications: Almost any further processing, such as the registration of partial scans into a complete shape, editing of sequences, or statistical analysis, requires dense correspondences between surface points.
Matching deformable shapes is in many cases a difficult problem: If we permit rather general deformations this might require many parameters that have to be explored during matching. The size of this search space is exponential with respect to the available degrees of freedom. However, for the important special case of a single object in different poses, we can often assume that the deformation is approximately isometric, i.e., preserving the intrinsic metric structure. Concretely, the distances along surfaces of objects such as humans, animals, plants, or cloths do not change a lot without serious injury or damage. This restriction leads to a strongly constrained search space. Lipman and Funkhouser
[10] argue that isometries between topological disks are a special case of conformal mappings, thereby limiting the degrees of freedom to six (three pointtopoint correspondences are sufficient). Ovsjanikov et al. [11] show that a single point correspondence is sufficient for a special class of shapes where the spectrum of the LaplaceBeltrami operator is not degenerate, thus showing that there are only two degrees of freedom in this special case. Additional cues, such as carefully selected constellations of local features [12], can even reduce the complexity for many shapes to the point of leaving a trivial search space, eliminating all degrees of freedom. As approximate isometry makes the correspondence problem feasible while still permitting significant pose changes, many of the recent shape matching algorithms are based on this assumption [13, 14, 15, 16, 17, 10, 11, 18, 12, 19].However, the isometric matching problem is not yet solved: Because of the intrinsic view of the geometry, it is naturally sensitive to topological noise. In case of holes, missing data, or contacts, intrinsic distances become distorted and thus no longer constitute invariants that can be exploited for matching. One solution is to replace the notion of distance. For example, by using diffusion distances, or variants thereof, one can reduce the sensitivity to topological artifacts [17, 20, 21] when the pieces of geometry that cause the problem are small in relation to the overall shape. Nonetheless, these invariants still break down in case of large artifacts (wide contacts, large holes, as shown for example in Section 5 of this paper). Unfortunately, realworld 3D scanner data, one of the main practical application areas, is almost universally troubled by substantial artifacts of this kind.
Formulated more generally, we have to address the problem of partial matching
, where not the whole manifold can be brought into correspondence by a (near) isometric mapping but only an excerpt of the surface can be matched. In this case, typical invariants (geodesic paths, Laplacian eigenfunctions) become unreliable because they utilize global information. Importantly, the excerpts of the surfaces that can be matched are not known a priori (otherwise, we could just restrict a traditional method accordingly) but need to be determined along with the matching. This seems to reintroduce prohibitive complexity as we now have to choose from an exponential number of such subsets
[22]. The most important contribution of this paper is to show that with an appropriate matching model this is not the case. The search space is not much larger than in the global problem and we give an efficient algorithm for computing such matches.The core of our method is based on the observation from differential geometry that an isometric map can be fully specified by using local information up to first order only: An isometry between Riemannian manifolds is fixed by a single point correspondence and an orthogonal map of the tangent spaces (see for example [23, p. 201]). In the case of surfaces, this means that the map of a point and a local direction (plus orientation, in case of unoriented surfaces) is sufficient to determine an isometry. We will sketch a constructive proof in Section 3 that directly yields a propagation algorithm for computing matches: Starting with a single oriented point match, correspondence information is incrementally propagated to the neighborhood, thereby floodfilling a partially consistent region of isometric geometry. Being local in nature, the method handles partial matches naturally and is robust to topological noise, which is reported naturally as boundary of partiality.
From a structural point of view, we can understand partial isometries among smooth, connected manifolds as equivalence classes of mappings of a pair of points and their tangent spaces on the two surfaces involved. In the case of oriented 2manifolds, each such object has six degrees of freedom (a point and an angle around the surface normal for each for source and target surface, respectively). The mapping is captured by a equivalence class of such 6tuples. This set is redundant in that one degree of freedom (the sum of the angles) can be chosen arbitrarily and two further parameters (either the starting or the end point) are only needed to select the partial region to be mapped. This leaves three degrees of freedom that need to be actually explored densely, with false starting points being rejected in time.
In this view, we can perform a relaxation: In order to also find approximate isometries, we can cast this problem as the task of finding approximate equivalence classes. Kim et al. [18] have used this idea in the context of globally consistent isometric mappings in order to efficiently find approximate isometries. In our paper, we demonstrate how our partial matching framework can be adapted to perform the same task in a partial matching scenario. We perform agglomerative clustering [24] in the space of nearisometric mappings, which are concisely represented using the tuples introduce above.
We validate our algorithm on standard benchmark data sets and raw scanner data, and compare the results to previous work. We show a significant improvement in quality over global methods in shapes with topological noise. Our algorithm yields similar or better results as previous heuristics for partial matching, but with stronger guarantees of discovering existing isometries as outlined above.
In summary, our paper proposes a systematic framework and new algorithms for extending isometric matching to the case of partial consistency, thereby making the following specific contributions:

We characterize partial isometries of shapes by singlepoint maps up to first order, which yields a tight bound on the inherent degrees of freedom.

This leads to a novel matching algorithm that provides a systematic approach to the general setting of partial intrinsic matching, where both surfaces may be incomplete, including robustness to strong topological noise.

By interpreting partial matching as a problem of finding approximate equivalent classes in our novel representation, we obtain an algorithm for approximate partial isometric matching.
Our algorithmic pipeline for approximate partial isometric matching is summarized as follows. We identify distinctive feature points on both surfaces. We compute oriented point correspondences by matching feature descriptors, with orientation determined by nearby features. From oriented point matches, we perform local metric propagation, stopping the propagation when the stretching becomes too large. This gives us a set of partial isometries, covering different, but possibly overlapping regions of both surfaces. We define a dissimilarity measure between partial isometries, and use this to cluster the partial matches into equivalence classes. The cluster (equivalence class) with the smallest total intracluster distance is merged by taking the geodesic centroid of all candidate correspondences for each point on the source manifold.
The steps in this pipeline exhibit sensitivity to certain challenges. Most prominent are the sensitivity of feature matching to surface noise, missing data and topological noise, and the sensitivity of the clustering step to undersampling of the space of isometric mappings, which can be exacerbated by problems in the feature matching stage, in addition to the lack of a proper distance metric between partial matches. The difficulty of reliable feature matching on real data can be alleviated by the metric propagation algorithm, which will typically produce smaller partial isometries for incorrect feature matches than for correct ones since we expect the local metric to be vary significantly for different parts of the surface. However, in the presence of strong continuous intrinsic symmetries this assumption breaks down. We thoroughly discuss and explore in which situations and to what extent these challenges create problems in the final result in Section 5.
2 Complexity of Isometric Shape Matching
In this section, we discuss different isometric matching models and their implications on the difficulty for finding a globally optimal solution to the matching problem. To the best of our knowledge, this has not yet been analyzed explicitly in literature. We do not aim to give an extensive review of surface registration methods; for this, we refer the reader to recent surveys such as van Kaick et al. [25] or Tam et al. [26].
We start by introducing some formal notations.
Manifolds: We consider smooth, orientable manifolds embedded in threedimensional space. In order to represent partial data (such as scans with acquisition holes), we permit boundaries, denoted by . The orientation of might (optionally) be prescribed by oriented surface normals , , pointing outwards.
Tangent space: We further use to denote the tangent space of at point
. For its representation, we choose two arbitrary but fixed orthogonal tangent vectors
, i.e.: .Distances: We use for to denote the intrinsic or geodesic distance between the two points and . A geodesic connecting and is a curve that has no geodesic curvature, which means that the derivative of the curve at a point projected to is zero. We call the shortest geodesic connecting and a shortest geodesic path in .
Mappings and isometries: Consider two manifolds and , and a mapping from an open subset to . Let and denote the partial derivatives of with respect to the tangent space directions and in of . The first fundamental form of at point is then given by:
The function is an isometry if and only if for all .
Equipped with this notation, we will now identify and analyze three classes of approaches for finding surface correspondences. The difference is in how the notion of approximate isometry is handled, leading to different complexity characteristics and algorithms.
2.1 Global Approximate Isometry
The global model assumes that geodesic distances are global invariants of the shape, being consistent at least up to an error margin that accounts only for a small fraction of the object size. This means, the energy
(1) 
must be smaller than . For two points Eq. (1) corresponds to an additive error of at most , i.e.
(2) 
The global consistency criterion is sometimes modified to allow for a multiplicative error of instead, as
(3) 
The former error model considers absolute errors, while the latter one considers relative errors.
In case of exact isometry, i.e., , the set of matching candidates becomes strongly constrained. The isometry assumption has been used to embed the intrinsic geometry of a shape in a Euclidean space using multidimensional scaling, such that embeddings of isometric shapes become identical, which facilitates shape recognition and matching [27, 28, 29]. Alternatively, the geometry of shape can be embedded into shape using generalized multidimensional scaling, thereby computing a crossparameterization directly [14].
Exact isometry results in a set of matching candidates with few degrees of freedom. Lipman and Funkhouser [10] have noticed that isometries are special cases of conformal maps, thereby having only the degrees of freedom given by the Möbius group, which are fixed by three pointwise matches on spherical topologies. Ovsjanikov et al. [11] have shown that even a single point match is sufficient to fix an isometry if the LaplaceBeltrami spectrum of is nondegenerate. In this paper, we exploit a different way to uniquely describe an isometry: fixing one point, a tangential direction, and the surface orientation is necessary and sufficient to specify an isometric mapping [23]. This provides the fewest possible degrees of freedom while still covering all cases including shapes with global intrinsic symmetries.
In case of small, global error margins, statistical triangulation algorithms can be applied that compute all correspondences from a few landmark matches [15, 16] or, similarly, by voting for several approximate solutions [10, 11]. Depending on the geometry of the shape, errors can become amplified so that a bit more than only the minimal set of initial correspondences are required [12]. Nonetheless, matching according to this model is efficient and can be considered a more or less solved problem.
The global approximate isometry criterion has also been employed to study matching partial overlaps of complete surfaces. van Kaick et al. [30] use a pair of features to define a map that captures a local geodesic region between the two features, and show how this descriptor can be used for partial matching. Here, using two features instead of one is crucial because two features encode orientation information on the surface, while one feature does not.
The drawback of the global consistency model is its sensitivity to topological noise. To make globally consistent approaches more robust w.r.t. topological changes, several approaches have been proposed to perform a bandlimited analysis in the eigenspace of the manifold’s LaplaceBeltrami operator
[17, 20, 21, 31]. Unlike prior approaches based on embedding the intrinsic geometry of the shape directly [27, 28, 29, 14], these approaches successfully handle small topological errors. However, they break down in the presence of large artifacts, such as wide contacts or missing pieces. A notable exception is the heuristic region growing approach by Sharma et al. [32]. It tries to match points with similar spectral signatures using an expectation–maximization framework, which has been shown to perform well in the presence of large contacts. Despite good practical performance, the algorithm is heuristic in nature and it remains unclear under which conditions it will find a correct solution. In particular, if local descriptors are not unique, the greedy region growing might fail and the EMalgorithm does not guarantee to recover a correct global match. In contrast, our region growing algorithm, which is based on propagation of metric information rather than potentially ambiguous descriptor matching, comes with the theoretical guarantee to find correct results for exact isometries while requiring an initialization from a small search space with very few degrees of freedom.2.2 Global Approximate Isometry in Partial Regions
The global consistency model is incompatible with the notion of partial matching, since distances have to be measured on the complete surfaces and , which might not be available. Xu et al. [22] modify the criterion in Eq. (1) by restricting paths to the partially matched region :
(4) 
Note that considers again absolute errors
(5) 
and can be modified to consider relative errors as in Eq. (3).
The drawback is that the shortest geodesic paths and thus the energy depend on the shape of the domain, which makes it difficult to optimize ; changing the domain influences which pairs of points are mapped in a geodesically consistent way. For this reason, Xu et al. restrict their method to considering geodesically convex regions , which are defined as regions where for any two points in , there exists a shortest geodesic path between and in . In this special case, both and measure shortest paths on and . The solution proposed by Xu et al. optimizes the scale of consistent regions and the consistent points separately, which leads to a rather complex algorithm. Further, the notion of a scale parameter is not canonically related to the original matching problem.
Sahilloglu and Yemez [33] consider the case where one of the surfaces is complete and the other an incomplete, deformed part of that surface. Using a coarse sampling and matching strategy between shape extremities, they can directly estimate a scale parameter between the two surfaces, which allows them to define a scaleinvariant isometric distortion measure. This results in onesided partial dense intrinsic matching up to an arbitrary scale. Their approach also allows matching of semantically similar, but nonisometric complete surfaces. Our approach allows both surfaces to be incomplete, at the cost that they must be scaled consistently beforehand. Given realworld scanners often provide data in known units, our method is compatible with the scenario of matching surfaces acquired with different modalities. While in this paper we do not explore the latter scenario, our approach would be compatible with this task given a reliable way to estimate scale change during metric propagation, possibly using shape extremities or other intrinsic features.
Bronstein et al. [34] introduced a general framework to evaluate partial similarity using Pareto optimality. In case of partial intrinsic shape matching, this method aims to find large parts of two surfaces that are similar to each other, where similarity is measured according to Eq. (2). In practice, the parts are found using generalized multidimensional scaling. Raviv et al. [35] use a similar technique to find partial intrinsic symmetries.
2.3 Local Approximate Isometry
Another common way to relax the requirement of exact isometry towards approximate matching is to maintain the metric tensor in a leastsquaressense. Again, let
be a mapping between two manifolds, where is the open subset of that should be mapped to in a distortion minimizing way. We can measure the distortions for example by minimizing a matrix norm of the Green deformation tensor (difference of the first fundamental form to identity):(6) 
This criterion is purely local and thus well suited for partial matching. It is worth noting that extrinsic elastic deformation techniques, such as [36, 37, 38] are closely related: They either include the preservation of the curvature tensor in the objective function to maintain the extrinsic shape, or apply Eq. (6) to the 3manifold of the embedding Euclidean volume [39, 40]. All of these methods are designed for partial matching.
The problem with both intrinsic and extrinsic elastic matching models is that the search space becomes very large, rendering any approach based on exhaustive search prohibitively expensive. The structure of the search space can be approximately understood by a linearized analysis. In order to understand the degrees of freedom of the local matching model, we can use the tool of modal analysis of such elastic models, first introduced by Pentland and Williams [41] to the field of computer graphics (the aim of their paper was actually to speed up the simulation of extrinsic elastic deformations of solid objects in threespace). Modal analysis represent the deformation of an object as a linear combination of the eigenmodes of the object’s vibration, which are found using a spectral analysis of a linearized deformation energy.
Typically, these energies are related to the Laplacian of the deformed domain, thus leading to eigenvalues that are only decreasing rather slowly. Many modes need to be retained to represent the space of lowenergy (i.e., plausible) deformations adequately. If the local matching is not very stiff, the size of the search space explodes. This is intuitive: Permitting local deformations creates a large variety of permissible shape variants, for example by adding different local dents and combinations of those everywhere. Due to this large search space, it is impractical to match shapes purely based on the deformation model of Eq. (
6). In order to reduce the large search space for elastic matching, existing methods use additional constraints, such as, template models [40], temporal coherence [39], or fairly large sets of coherent feature correspondences [37].Recently, Kim et al. [18] proposed a new paradigm for the elastic matching problem. First, the method computes multiple dense maps between two shapes by assuming a global nearisometric deformation model. Multiple maps are obtained by fixing different triples of corresponding points for the computation of global conformal maps [10]. Second, the method computes local weights for pairs of maps, depending on their local adherence to isometry, and performs a spectral analysis to combine multiple weighted maps into a single global map. The key idea is to find cliques of similar mappings by clustering approximately compatible maps. The method was shown to perform well in many interesting cases. However, it cannot handle partial mapping; in particular constellations with large topological noise cannot be handled: Each global conformal map can is highly distorted due to the lack of global consistency. This introduces distortions in many of the local matches, and it is not always possible to remove these distortions in the final blended result (as demonstrated in Section 5).
Our clustering method (Section 4.4) is based on the same idea, but it combines partial maps instead of global maps
, therefore avoiding the mentioned problems of artifacts due to only partial consistency. On the technical side, the main challenge is that we cannot measure the distance between all pairs of candidate maps but only between actually overlapping ones. This is a problem for the original spectral clustering, which we substitute by an alternative technique geared towards sparse pairwise constraints
[24].2.4 Previous Work on Local Isometry
In previous work, there have been a number of attempts to find local isometric mappings, similar to our approach. However, these did not consider mappings between general surfaces but only local planar parametrizations.
Different techniques have been proposed to locally parameterize the intrinsic geometry of a surface to a plane using a local approximate isometry criterion. Schmidt et al. [42] used exponential maps to transport a local coordinate system along a surface for the purpose of texture mapping. More recently, Schmidt [43] used transported exponential maps to produce a parameterization of a local surface patch to the plane that has low metric distortion. The local surface patch is provided through user interaction as an input stroke on the surface. Malvaer and Reimers [44] propose an alternative parameterization based on an extension of polar coordinates to surfaces. In our work, we crossparameterize local surface patches from to . This crossparameterization is more challenging than a parameterization to the plane due to the arbitrary geometry of . Our crossparameterization task is further complicated as in our application scenario, both and are scanned point clouds with missing surface information and scanner noise. Hence, the reviewed methods cannot be applied in a straight forward manner in our application.
3 Local Metric Matching
We saw that for general surface matching using a global approximate isometry criterion for partial matching (Eq. (4)) is difficult since the domain changes, and that using a local approximate isometry criterion (Eq. (6)) is difficult since this results in a large search space. In this section, we outline our new method: starting from a point and attached direction in with known correspondence and corresponding attached direction in , we find the largest domain such that Eq. (5) is satisfied.
The key assumption of our approach is that and can be matched using a global approximate isometry in some partial region containing , implying that Eq. (5) holds. Our goal is to find the largest region for which this assumption is satisfied. This assumption holds in scenarios where we know that and are actually related by a nearisometric map but the data does not comprise all of the original input and/or contains additional unrelated geometry or contacts. The most important practical example where this assumption holds is the acquisition of a surface that deforms nearisometrically but the scanning equipment introduces areas of missing data (shadowed to the scanner by other object parts) and cannot correctly resolve the topology in contact areas (such as the hand of a person being in contact with the body).
Our approach can be viewed as a hybrid approach between global and local approximate isometric matching. We use the assumption that and can be matched using a global approximate isometry in some partial region , and we compute by growing a region using a local approximate isometry criterion. This allows us to combine the advantage of the global methods of having a lowdimensional search space with the advantage of the local methods of being well suited to describe partial isometric matches.
Let denote the parameter domain of all partial isometric matchings between and . Our goal is to compute all partial isometries that map maximal subsets of to . The vector parametrizes the set of all such mappings.
In the following, we will discuss that:

Isometric deformations of have three degrees of freedom.

A partial isometry can be parametrized by , where denotes the unit circle.

There exists a global representation redundancy that identifies all choices of parameters and .
3.1 Three Degrees of Freedom
For a (sufficiently) smooth Riemannian manifold fixing one point on , a tangential direction in , and a surface normal at suffices to specify an isometric mapping [23]. For details on the smoothness criteria, refer to [45]. The following proof sketch aims to give the intuition behind this statement for 2dimensional surfaces embedded in by showing that we can transfer the local coordinate frame defined by and to any point on in a canonical way.
We start with two definitions. The injectivity radius at is the largest radius, such that for any point on with geodesic distance less than from , there is a unique shortest geodesic path between and . The injectivity radius of is defined as .
For closed surfaces, the injectivity radius can be bounded from below by the minimum of , where is the Gaussian curvature, and half of the length of the smallest periodic geodesic [23, Thm. 10]. It follows that holds for closed Riemannian surfaces with finite Gaussian curvature. For general surfaces with boundary, the injectivity radius may become zero. Note that in this paper, we consider Riemannian surfaces with finite Gaussian curvature with boundaries. The boundaries are present because the closed Riemannian surface of interest is only partially observed by the acquisition device. In this special case, the injectivity radius is still positive.
Imagine that is covered by overlapping regions of intrinsic radius less than . These regions are all topologically equivalent to disks. Consider a disk containing a point as shown in Figure 1. In a small neighborhood of , is arbitrarily close to the tangent plane . Note that , , and fix a local orthonormal coordinate frame in . This coordinate frame can be transported to a point along the (unique) shortest geodesic path from to in by parallel transport, which can be thought of as repeatedly projecting the direction to the tangent planes of consecutive points along in infinitesimally small steps (for details on parallel transport see for example Berger [23, Chapter 3.1]). Let denote the transported direction. The transferred direction lies in the tangent plane , and can again be used to fix a local orthonormal coordinate frame at . This transfer can be repeated by chaining together disks until every point on was assigned a fixed tangential direction by parallel transport along a path connecting and that consists of an arbitrary but fixed sequence of (unique) shortest geodesic paths within the chained disks. Note that by construction, we only use intrinsic information to propagate the direction to the entire surface. Hence, the resulting local coordinate frames are invariant with respect to isometric deformations of when encoded relative to fixed local coordinate systems on .
This implies that any isometric deformation of an oriented surface can be specified using three degrees of freedom: one point on (accounting for two degrees of freedom) and a direction in the tangent plane of .
3.2 Representation
We can use the fact that any isometric deformation of an oriented surface can be specified using three degrees of freedom to derive a representation for intrinsic mappings. Specifically, identifying one corresponding point and one corresponding tangential direction completely determines an isometric mapping between two oriented surfaces. More formally, to define an isometric mapping between (subsets of) and , it suffices to specify a point , its intrinsically corresponding point , a tangential direction in , and its intrinsically corresponding tangential direction in .
Starting from this information, and assuming that and are isometric, we can propagate the correspondence information by mapping the metric structure of onto as follows. Starting from the corresponding points and along with the corresponding directions and , we can propagate the correspondence information to a sufficiently small geodesic neighborhood of by simultaneously walking along corresponding geodesic paths starting at and , respectively, and by matching points that are reached at the same time. Here, is sufficiently small if its geodesic radius is below . Once the correspondence information is computed for , we continue propagating the correspondence information from a point close to the boundary of to its geodesic neighborhood, and iterate until every point on has a correspondence.
The assumption that and are (near) isometric can also be used to detect the boundary of the largest region containing for which the mapping is nearisometric. The reason is that the propagation algorithm allows us to measure the difference in intrinsic geometry in newly mapped parts of and directly. Hence, we can stop the region growing algorithm if a newly added correspondence would induce a stretching larger than .
Figure 2 illustrates the nearisometric region growing process. The plane on the left () and the plane with a hill on the right () are isometric except for the hill. Beginning with a point and direction match in Figure 2a, the isometric region grows outward in all directions where the isometry condition is locally satisfied. This way, the largest nearisometric partial match is identified, as shown in Figure 2d. Note that this region can have a complex topology and geometry.
3.3 Representation Redundancy
The representation discussed above contains redundant information: Infinitely many may represent the same nearisometric mapping.
The first degree of redundancy is the choice of the direction in the tangent plane . Changing this direction merely rotates the field of directions in (and similarly, the field of corresponding directions in ). Hence, the choice of this direction has no influence on the final result. To remove this redundancy, we start from an arbitrary but fixed direction and precompute and fix for all on .
The second redundancy is the choice of the start point . Let be the maximal set within which an isometry can be constructed. We can replace in by any other point . See Fig. 3. This requires an update of the remaining parameters. The direction is set to the precomputed value (see above), and the remaining parameters can be updated using the computed mapping . More specifically, and the direction of is set to the direction of , where is set small enough that is in . Thus, in the case of a global isometry, or when we know beforehand the isometric region , the mapping has three degrees of freedom. In the general partial isometry case, however, where is unknown, the starting point is not fully redundant; it still selects the equivalence class that represents a certain partial map; computationally, this is the patch to be matched by the propagation algorithm. This does not increase the complexity strongly as we just have to restart the matching in case multiple partial matches exist. Ideally, we would sample one starting point per partial isometric region, which in practice will be far fewer than there are samples on . While we do not determine this lowerbound beforehand, we maintain low complexity by using features to identify potential starting points, and marking a starting point as redundant if another starting point produces an isometric region containing . In case of a mismatch, it can be discovered quickly if the target area does not match, which will become evident within constant time.
In summary, all mappings represented by form an equivalence class in the parameter space . Since we can remove one degree of redundancy by fixing the tangential directions on , for each mapping we have two degrees of freedom that vary among equivalent maps (the choice of the start point on ), which along with the manifold structure of implies that each equivalence class forms a manifold in , which can be computed directly using the propagation algorithm introduced in Section 3.2.
In practice, we can take advantage of this representational redundancy as follows. Since and are discretized and corrupted by noise, the error of the correspondence information computed using the propagation algorithm increases with increasing distance from the start point of the propagation. Hence, it is possible that the propagation stops prematurely due to discretization artifacts and the influence of noise, thereby identifying a region that is smaller than the correct solution. Thanks to the redundancy in the representation, we can start the propagation algorithm from multiple oriented point pairs, identify a set of equivalent mapping functions , and them into a single mapping function that covers a larger area of and is less influenced by noise than the individual . The following section discusses a direct implementation of this theoretically motivated method, which we will use to compute correspondences of noisy scanner data.
4 Pairwise Intrinsic Matching
From the preceding analysis, we derive an algorithm for computing partial nearisometries between two surfaces and . We start our discussion by outlining how surfaces are represented and how basic steps of the algorithm are implemented (Section 4.1). Our direct, nonoptimized implementation is based on enumerating the nonequivalent choices by selecting different oriented point matches (Section 4.2), growing the isometric region locally from there until no more points locally satisfy (Section 4.3), and finally clustering the partial maps into equivalence classes (Section 4.4). Figure 4 gives a graphical overview of our matching pipeline.
4.1 Surface Representation
In our implementation, the surfaces and are either represented as oriented point clouds or meshes. While most parts of the algorithm can be extended to discretized surfaces in a straight forward manner, for some parts of the algorithm we require a continuous surface representation. We obtain a continuous surface representation as implicit movingleastsquares (MLS) manifolds using the robust method of Öztireli et al. [46]. Using this method, a discretized surface is represented continuously as the zero levelset of an implicit function derived from the oriented vertices of .
Using a MLS representation allows a point to be projected to a continuous representation of in the case where is given as oriented point cloud or mesh. In the following, whenever we refer to projecting a point onto a discretized surface , this projection is implemented as projecting to the MLS representation derived from .
One basic operation needed by our algorithm is the computation of geodesic distances and paths on . In some parts of the algorithms, rough estimates of geodesic distances and paths suffice, and these are computed using Dijkstra’s algorithm. In other parts of the algorithm, it is crucial to have accurate estimates of geodesic distances and paths. In these cases, we initialize a geodesic path to the Dijkstra path and refine by minimizing the length of using the constraint that must not leave . This optimization is carried out iteratively using a conjugate gradient method. After each step, is projected back to . In the following, we refer to these refined geodesic distances and paths as smoothed geodesic distances and paths.
A core part of our approach is to compare distances measure on the target surface to distances measured on the source for corresponding points. Comparing geodesics between all pairs of correspondences quickly becomes prohibitively expensive and would limit the practical applicability of our algorithm. To remedy this, we construct a topology hierarchy on similar to [47] as follows. We define level of the hierarchy to be the original set of vertices and their connectivity–either the original triangle mesh or the nn graph for a point cloud. (In all our experiments involving point clouds, we set .) The sample spacing is defined as the average edge length. Level of the hierarchy is constructed by selecting an evenly spaced subset of the level vertices and connecting vertices in a topology preserving way. Subsequent levels are constructed in the same way as a subset of level . At each level the subset for the next level is determined by doubling the desired sample spacing, . At a coarser level of this hierarchy fewer vertices are connected to any others, but at a greater distance to each other. In our algorithm, we only consider geodesic distances between vertices that share an edge in at least one level of the hierarchy, which results in a sparse set of distances to be optimized, even for a dense set of correspondences.
This hierarchical structure ensures that locality is respected, as only geodesics between points that are neighbors in some level of the hierarchy are considered. Having a structure that respects locality is crucial to allow for partial matching. The number of levels in the hierarchy determines the tradeoff between local and global isometry constraints. We use levels in all experiments in this paper.
Note that during the execution of the matching algorithm, the vertices on are fixed, as we aim to find a correspondence for every vertex of on (if such a correspondence exists). Hence, we can precompute all geodesic distances on for vertices connected in some level of the hierarchy. This way, only distances on need to be updated during the metric optimization.
4.2 Finding Oriented Point Matches
In principle, by exhaustively trying all possible starting points and directions, our algorithm will recover all partial isometries. However, since this is infeasible, we reduce the search space using sparse feature matches. This section outlines an algorithm to find features matches. However, note that this part is not a novel contribution of this paper and only given for completeness; in principle, any feature matches can be used to reduce the search space, as demonstrated in Section 5, where we show a result using imagebased feature matches from a multiview capture setup. We denote the feature descriptor at a point with a vector , and similarly for points on .
We identify features on and using the geodesic fingerprint descriptor [48], which compares how an isocontour of the geodesic distance from a point deviates in length from the isocontour of the same 2D Euclidean distance. We use isocontours in all experiments, for geodesic radii between and . The values of and need to be varied slightly depending on the amount of noise in the input data, and are discussed in Section 5. We define the distinctiveness of a point as the sum of distances to the rest of the vertices of in descriptor space. Features are points that locally maximize distinctiveness. The leftmost box in Figure 4 shows features colorcoded by distinctiveness (red for most distinct, blue for least distinct).
To find feature matches, we begin by computing the Cartesian product of descriptor distances between features on and . The vast majority of these are not correct mappings between the surfaces. We filter these potential feature matches using both the distances between descriptors and the distinctiveness of the features. More precisely, for each feature on , we only consider the features of that have the closest descriptor matches to . (We use in all our experiments.) The initial dissimilarity between two features and is defined as
(7) 
where is a weight. (We use in all our experiments.) This dissimilarity is minimized for features with high distinctiveness that are similar in descriptor space.
This procedure often produces a good set of sorted feature matches on clean data. For noisy data from real scanners, however, the descriptors will be less discriminative, and considering spatial relations between features in addition to descriptors and distinctiveness leads to more reliable feature matches.
To do this, we iteratively build (possibly overlapping) clusters of consistent feature matches. In each iteration, the feature match with the next lowest is chosen as a starting point for . The cluster is built by repeatedly adding the closeby feature match that has the lowest dissimilarity to . More precisely, let . In the next step, all features that are neighbors of in the topology hierarchy are considered, along with their potential matches . The dissimilarity between a feature match and cluster is defined as
(8) 
where is a weight. (In all our experiments, we set to one over the square of the average edge length on .) We repeatedly add the match with smallest dissimilarity to as long as is below a threshold. (We use threshold in all our experiments.) We stop adding new clusters once the sum of the cardinalities of the clusters exceeds the initial number of feature matches.
Note that the above clustering scheme is equivalent to modeling both the descriptor distances and the summed stretches of geodesic distances of a feature match to a cluster as normally distributed. Hence, the above stopping criteria correspond to stopping the clustering once the joint probability of a match belonging to a cluster becomes small.
After the clustering, we place the feature matches in a minpriority queue, so that we start isometric region growing first from the matches we expect to be most reliable. A feature match is assigned priority if is not part of any cluster and priority otherwise. We repeatedly take the minimum element from queue and use it to generate partial isometric mappings using the region growing of Section 4.3.
However, so far we have only established a positional correspondence between features, and we need a directional correspondence as well to fix the partial isometry we wish to grow. To establish direction, we build a minpriority queue as outlined above, but only with the subset of features in the neighborhood of the current positional match in the topology hierarchy. Matches of neighboring features allow us to find corresponding tangent plane directions from corresponding smoothed geodesic paths between feature points.
Once a partial isometry has been computed, we increase the priority of feature matches that are redundant given the already computed partial match. If the same source and target points are already matched, by our model, the result of running the growing again from those same points will be equivalent.
To keep the runtime of our method bounded, we stop after a fixed number of oriented point matches have been tried. (We use in all our experiments.) A more general stopping criterion could be devised based on determining the maximal coverage of and subject to a consistent mapping. The second box from the left of Figure 4 shows three oriented feature matches found using this method.
4.3 Isometric Region Growing
Starting from an oriented point match , we grow the region by adding matches incrementally in the local neighborhood of the boundary of . For a new point near the boundary of such that , let denote the neighborhood of in level of our topology hierarchy. We use parallel transport along corresponding directions in and emanating from an oriented point match , where is in a local neighborhood of . We know the full path between and on , and from and we know the corresponding start direction on . Hence, we can transport the start direction along corresponding paths on and until we have traveled . In practice, we implement parallel transport by moving in small steps in the tangent plane, and by reprojecting the resulting point to the surface and updating the tangent to the path. We use smoothed geodesic paths to transport the position of the new match on to reduce discretization errors for differently sampled surfaces. For robustness, we use all available oriented point matches in and take the Riemannian center of mass [49] of the transported points, where all transported points have equal mass.
The newly matched source point is added to only if it locally respects the stretch factor as given in Eq. (5). This is necessary for two reasons. First, subsequent matches will be estimated based on the assumption that the existing matches are nearisometric. Second, as discussed next, we use a nonlinear leastsquares optimization to refine the positions of the matched points on , which effectively distributes the error evenly over the matched region. To avoid introducing errors, it is therefore important to verify that each newly added point match locally respects the stretch factor .
To reduce the effect of quantization errors and noise, we optimize the metric matching of using a nonlinear optimization technique [14] every time the area of has doubled. This means more frequent optimizations at the start of the growing process. This optimization reduces the amount of drift as it realigns the matched regions while taking into account long geodesics, between neighbors in the top level of our topology hierarchy, in . For increased efficiency, we only consider edges in the topology hierarchy of (explained in Section 4.1) during the optimization.
An important difference to the global approach used by Bronstein et al. [14] is that we optimize the metric only using geodesic paths which are entirely within and . This models the isometry criterion in partial regions and is crucial to handling topology changes and missing data. Note that such defects may cause large differences between and , as well as and , respectively. The second box from the right in Figure 4 shows the first three partial isometries found by growing isometrically from oriented feature matches in our priority queue.
The proper value for the stretching threshold depends on a number of factors: material properties, the resolution at which the surface is sampled (as it affects the accuracy of the Dijkstra paths), and the noise of the acquisition process. In our experiments, we do not consider material properties, and we assume that the acquisition noise has an equal influence on both source and target. Hence, we set to to account for quantization effects in the computation of geodesics.
4.4 Combining Equivalent Partial Maps
It remains to identify and merge a set of partial mappings that represent the same mapping function . If the surfaces were related by exact isometries and noise was negligible, the following step would not be required. However, for realworld data, the identification of functions that are approximately equivalent improves the results substantially.
The problem at this point similar to the blending problem by Kim et al. [18]. Recall that their approach uses a spectral method to find blending weights for different maps. This is a good approach in their case as blending weights are given as the solution to a quadratic energy function. Note that since in our model, equivalent mappings form a 2manifold in parameter space , and the parameter space is nonlinear, it is not appropriate to use a spectral approach.
However, we can take advantage of the property that equivalent mappings form a 2manifold in . We employ the agglomerative clustering algorithm of Zhang et al. [24]
for discovering manifold structures in highdimensional data based on the indegree and outdegree of the nearestneighbor graph of points in highdimensional space. In our case, we consider the nearest neighbor graph of partial isometric mappings, where the dissimilarity between these points in
is measured as follows.We compare different maps and using a dissimilarity measure based on their domains and :
(9)  
where , , , and denotes the surface area. In practice, we compute a discrete version of this dissimilarity by replacing integrals over regions by sums over vertices in the region. We cluster different mappings together until the maximum affinity between any two clusters is greater than a threshold . Affinity is computed from the weighted graph degree between mappings, where the weights have a doubleexponential falloff as increases. See Zhang et al. [24] for details. While the direct relation of to the allowed stretch is not easily defined, it should be lower when we want to allow for greater stretching between partial maps. We only have to adjust this value in a few cases in our experiments, as discussed in Section 5. Following the clustering, we select for our final mapping the cluster with the highest intracluster affinity, or connectivity, as proposed by Zhang et al. [24]. This most often correlates with the cluster that covers the largest portion of the surface.
We merge clustered maps by computing a weighted geodesic average on for each source point in the union of the mapped regions as
(10) 
where are the clustered partial isometries, which we want to merge into
. The weights are computed as an exponential distribution in the geodesic distance from the starting point match
as , reflecting that we expect errors to accumulate as the growing proceeds because of discretization artifacts, data noise, and deviations from isometry. We set , where is the sample spacing of the level of the topology hierarchy described in Section 4.1. Note that this is equivalent to finding the Riemannian center of mass [49] of the estimates . The rightmost box of Figure 4 shows the final clustered and merged result.5 Experiments
To validate our theoretical analysis, we evaluate a direct implementation of our algorithm and compare our results to state of the art approaches. In the following, we demonstrate that our method achieves results that are either comparable or superior to the state of the art, which demonstrates that our algorithm not only has theoretical advantages, but is applicable in practice as well.
5.1 Implementation Details
We implemented the algorithm described in Section 4 in C++, and conducted our evaluations on a standard laptop PC. During the evaluation, all but two types of parameters are fixed. The first type of parameters that is varied is , which controls the size of the neighborhood used to compute surface descriptors. If the radii are set higher, then the method is more robust with respect to noise at the cost of potentially missing features of small scale. In our experiments, we only use two settings for these parameters. The first setting, and , is for relatively clean data, and the second setting, and , is for noisy data. Here, is set to of the diameter of . The second parameter that is varied is the threshold used to control the clustering. Lower values of allow for less isometric patches to be clustered together. We vary in our experiments.
For the examples discussed below, our algorithm takes between 30 minutes and 8 hours to compute the final result. To give an idea of the distribution of the time, we discuss the running time for one pair of models (the spacecarved samba models shown in Fig. 9) in more detail. For this pair, finding oriented feature matches takes about 2 minutes, growing partial mappings takes about 1.5 hours, clustering the mappings takes about 9 minutes, and merging the patches takes about 44 minutes. Hence, the total time to compute the results is about 2.4 hours. Note however, that the running time of our method depends significantly on the distinctiveness of the intrinsic geometry of the surfaces, relative to the noise level. For example, the templatefitted samba models shown in Fig. 6 take about 1.5 hours in total, despite having more than twice as many vertices as the spacecarved versions, because the approximate isometry criterion is more discriminative (geodesic distances are less perturbed by surface noise).
5.2 Comparison to State of the Art
We compare the performance of our algorithm against four existing methods, namely heat kernel maps (HKM) [11], blended intrinsic maps (BIM) [18], the method of Sharma et al. [32], and the method of Tevs et al. [50]. We select these methods for comparison because they represent the state of the art for matching between surfaces. More specifically, HKM is derived from the theoretical complexity of isometric mappings, BIM can match surfaces that exhibit local deviations from isometry, and the methods of Sharma et al. and Tevs et al. are heuristics that have been specifically designed for matching partial data with topological noise. For HKM we use our own implementation, which uses two feature correspondences to initialize the mapping, for BIM we use the authors’ implementation, for comparisons with Sharma et al., we run our code on their data, and for the method of Tevs et al., the authors were kind enough to run their algorithm on our data.
We show comparative evaluations on a variety of types of data. The first type is a synthetic data set that helps in understanding the major difference between our approach and HKM, and illustrates the partial isometric matching model. The second type is data acquired using either a laser scanner or an imagebased reconstruction system that was processed by fitting a template to the data. In this case, we treat the result of the template fitting as ground truth. The third, and most challenging, type is unprocessed realworld data acquired using different acquisition systems.
To compare our approach to previous methods, we use two evaluation methodologies. In cases where ground truth is available, we compare quantitatively by evaluating the accuracy of different results with respect to the ground truth. For data that has no ground truth, we rely on visual evaluation. Our visualization scheme is as follows. A texture on is mapped to corresponding points on , and regions of that have no correspondence in are colored red. As a texture we combine constant coloring of semantically distinct parts (where applicable) with a checkerboard pattern. This type of texture simultaneously shows both global semantic accuracy and finescale distortion.
5.3 Synthetic Data
We start by comparing our algorithm against HKM using the synthetic example shown in Fig. 5, where we map from a plane to a plane with three peaks. This experiment illustrates the difference between a global isometric matching model and partial one: the two surfaces are globally nonisometric, but the planar part is isometric. In this example, we use two ground truth correspondences to initialize both algorithms to remove differences due to different feature matching approaches. Since HKM aims to map the shapes using a global isometry, the result maps planar parts to the peaks, while our method successfully detects the largest part of the surfaces that can be isometrically mapped.
5.4 TemplateFitted Scan Data
Next, we consider acquisitions of realworld data that was processed by fitting a template shape to the raw data. For all experiments in this section, we use and .
We first use two frames of the samba sequence by Vlasic et al. [51]. These frames are locally very close to isometric, but globally have high nonisometric distortion due to the dress being connected to the legs. We therefore set in the clustering step. Vlasic et al. provide a processed version of the frames, where a template was fitted to the data. This processing ensures that the models have the same topology, and the processed data can be used as ground truth correspondence.
We compare our method to HKM and BIM using the processed frames of the samba sequence. Figure 6 shows the results. Note that while HKM leads to a result with visual artifacts, the results of BIM and our method are visually pleasing. Furthermore, since we have ground truth correspondences, we show the cumulative error distributions for all three methods in Figure 7. Here, geodesic error is measured as a fraction of the square root of the surface area of . Note that our method numerically outperforms the two other methods.
The experiments conducted so far have shown that HKM, while being based on a solid theoretical foundation, leads to results of low quality when the aim is to find dense correspondences in datasets that contain noise and nonisometric distortion. Hence, in the following, we exclude this method from our comparisons.
Second, we test our algorithm on the SCAPE dataset [52] consisting of 71 scans of a male scanned in different postures. In our experiment, we match the neutral posture to all 70 remaining postures using BIM and our approach. The cumulative error distributions for all 70 mappings are shown in Figure 8. This data differs from the samba frames in that it is globally nearisometric, but contains local areas of highdistortion (at joints for example). For this reason it provides a different kind of nearisometric test, and poses a greater challenge for our method, which does not exploit global assumptions. This is reflected in our error curve being below that of BIM. We set in this experiment to reduce the influence of partial isometric patches that were thrown off by high local distortions.
5.5 Raw Scan Data
Finally, we consider raw scan data acquired using different acquisition systems. This type of data is noisy and incomplete, and is therefore significantly more challenging to match than the data used in the previous experiments. By using data from a variety of acquisition systems, we show our methods robustness to different types of acquisition noise. We also show our method’s robustness to other acquisition artifacts that violate global isometry: large holes and contacts.
First, we use the same two frames of the samba sequence by Vlasic et al. [51] that were used above. However, this time, we consider the geometry reconstructed by a spacecarving algorithm rather than template fitting. In this case, the frames we choose have different topology (as the hands merge with the body at the hips in one of the models) and are severely corrupted by noise. For this reason, we set , , and . We use these models to compare our method to BIM, Tevs et al.’s method, and Sharma et al.’s method. Figure 9 shows the results. Note that the results using BIM and the method by Tevs et al. match parts of the body of to the arms and legs of , thereby leading to visual artifacts. (See enlarged areas in Figure 9.) The method of Sharma et al., which is especially well suited to the scenario of handling contacts, leads to a result that covers the surfaces well. For this example, our method is run using two feature sets: first, using the standard features of our method and second, using the same imagebased features used by Sharma et al. Our method produces a visually accurate mapping in both cases. However, when using standard features, the result of our method does not cover the right foot of the target surface, while the entire target surface is covered when using imagebased features. Note that both the result by Sharma et al. and our results detect the contacts correctly and stop the growing in these regions. Hence, all areas of the surface are matched well. The improved performance with imagebased features illustrates some of the technical challenges. The high level of surface noise makes matching geometric features difficult, which results in poorer sampling of the space of isometries, which in turn results in a poorer clustering result. We note however, that using the same features as Sharma et al., we obtain equal coverage and accuracy, and that when using purely geometric information, we outperform the other purely geometric methods tested.
Second, we compare our method to BIM and Tevs et al.’s method using two models of the BU3DFE face database [53]
. The two models contain numerous small holes and outliers. Furthermore, the models have different topology because the mouth is closed in one model and open in the other one. For these reasons, we set
, , and . Figure 10 shows the results. Note that our method obtains a mapping with higher visual accuracy than both BIM and the method of Tevs et al. In the case of BIM, this is to be expected, since the topological change breaks the global isometry. (Note the lips mapped to the side of the face.) The method of Tevs et al. produces a much more accurate result, but with still significant overall distortion and outliers (specklelike effect). Our method leverages local information to fix partial isometries, and produces a mapping that largely preserves the semantic coloring with significantly lower distortion and without outliers.Third, we compare our method to that of Tevs et al. on two point clouds acquired using a laser scanner. The two models contain numerous holes and outliers. As the models are point clouds and BIM requires input meshes, we do not compare our result to BIM for this experiment. For these models, we set , , and . The results are shown in Figure 11. As can be seen, the method of Tevs et al. obtains better coverage, however our method has fewer outliers–islands of incorrectly mapped points within larger smoothly mapped regions (see the enlarged parts of Figure 11). The suboptimal coverage of our method is due to the difficulty in matching features on surfaces plagued by missing data. A feature descriptor less sensitive to holes could improve this.
5.6 Limitations
In the previous sections, we have demonstrated that our method not only has theoretical advantages, but also computes results that improve upon the state of the art results in challenging cases where the input data is a pair of raw scans with topological noise.
We now discuss some limitations of our algorithm. Our algorithm is based on growing nearisometric mappings between partial regions of two surfaces and then clustering consistent mappings together. In this way, our algorithm enumerates a set of nearisometric mappings. For models that exhibit a large number of partial intrinsic symmetries, this technique enumerates a large set of nearisometric mappings where many of the mappings are inconsistent with each other. Since we stop growing new nearisometric regions after the first 200 oriented feature point matches have been considered, for shapes with a large number of partial nearisometric symmetries, it may happen that there is no cluster of consistent mappings covering a large area of .
An example where the clustering step fails to identify a cluster of consistent mappings covering a large area of is shown in Figure 12. The two frames of the flashkick dataset [54] contain a large topological change due to a merge of the subject’s pants in one of the models. Note that our algorithm correctly stops the growing of single partial mapping in this area, as shown in Figure 12. However, since the legs and core of the target body are intrinsically symmetric (similar to cylindrical), many inconsistent partial mappings are found by the algorithm, and they cannot be clustered in a consistent way.
We should also note that the computational costs are quite high. This is partially due to unoptimized code, but an algorithmic shortcoming is the rather simple featurematching algorithm for finding start positions and directions. Optimizing this was not the focus of this work; our method should rather be understood as an alternative for the dense matching step where it can replace previous approaches based on conformal maps ([10, 18] and followups), heatkernel maps ([11] and followups), or geodesic triangulation with landmarks ([15, 16] and followups) in order to handle partial matching.
For future work, to address this limitation, we plan to combine our approach with an approach that detects continuous intrinsic symmetries [55], thereby reducing the search space for the initial feature matching and allowing us to efficiently enumerate all partial matches. Similarly, our framework could be extended to find partial intrinsic symmetries within a single object.
6 Conclusion
We have analyzed the complexity of the isometric matching problem under global and local isometry assumptions and based on this analysis we have introduced a new approach to solve the partial isometric matching problem using a representation for partial isometries that is both lowdimensional and redundant. Underpinning this is the fundamental observation that isometric mappings can be determined using purely local information and have only three degrees of freedom on 2manifolds. The local metric propagation algorithm we derived from this observation is designed to handle topological noise that could affect large portions of the model, including both large holes and contacts. The redundancy in the representation can be exploited to increase robustness of and to combine partial matches.We have shown how a direct implementation of this theoretical framework can be used to match challenging surfaces with different types of topological noise.
The insights gained by studying the partial isometric matching problem have the potential to impact other shape processing tasks. For instance, the representation for partial isometric matches introduced in this paper can be used to derive new algorithms to detect partial symmetries of shapes. For future work, we plan to further investigate this option.
Acknowledgements
The authors thank Art Tevs, Aurela Shehu, Waqar Khan and Avinash Sharma for their help in conducting the comparative evaluation, Vladimir Kim for making his code available, and Art Tevs, Silke Jansen, Alexander Berner, QiXing Huang, and Leonidas Guibas for discussions. We also thank the anonymous reviewers for their valuable comments that we used to improve the paper.
References
 [1] J. Davis, D. Nehab, R. Ramamoorthi, and S. Rusinkiewicz, “Spacetime stereo: A unifying framework for depth from triangulation,” IEEE PAMI, vol. 27, no. 2, pp. 296–302, 2005.
 [2] S. König and S. Gumhold, “Imagebased motion compensation for structured light scanning of dynamic scenes,” Int. J. of Int. Sys. Tech. App., vol. 5, no. 3/4, pp. 434 – 441, 2008.
 [3] T. Weise, B. Leibe, and L. V. Gool, “Fast 3d scanning with automatic motion compensation,” in Proc. CVPR, 2007, pp. 1–8.
 [4] D. Vlasic, P. Peers, I. Baran, P. Debevec, J. Popović, S. Rusinkiewicz, and W. Matusik, “Dynamic shape capture using multiview photometric stereo,” TOG, vol. 28, no. 5, pp. 174:1–174:11, 2009.
 [5] P. Debevec, “Virtual cinematography: Relighting through computation,” IEEE Computer, vol. 39, no. 8, pp. 57–65, August 2006.
 [6] V. Blanz and T. Vetter, “A morphable model for the synthesis of 3d faces,” in Proc. SIGGRAPH, 1999, pp. 187–194.
 [7] B. Allen, B. Curless, and Z. Popović, “The space of human body shapes: Reconstruction and parameterization from range scans,” TOG, vol. 22, no. 3, pp. 587–594, 2003.
 [8] N. Hasler, C. Stoll, M. Sunkel, B. Rosenhahn, and H.P. Seidel, “A statistical model of human pose and body shape,” CGF (Proc. Eurographics), vol. 28, no. 2, 2009.
 [9] T. Weise, S. Bouaziz, H. Li, and M. Pauly, “Realtime performancebased facial animation,” TOG, vol. 30, no. 4, pp. 77:1–77:10, 2011.
 [10] Y. Lipman and T. Funkhouser, “Möbius voting for surface correspondence,” TOG, vol. 28, no. 3, pp. 72:1–72:12, 2009.
 [11] M. Ovsjanikov, Q. Mérigot, F. Mémoli, and L. J. Guibas, “One point isometric matching with the heat kernel,” CGF, vol. 29, no. 5, pp. 1555–1564, 2010.
 [12] A. Tevs, A. Berner, M. Wand, I. Ihrke, and H.P. Seidel, “Intrinsic shape matching by planned landmark sampling,” CGF, vol. 29, no. 2, pp. 543–552, 2011.
 [13] D. Anguelov, P. Srinivasan, H.C. Pang, D. Koller, S. Thrun, and J. Davis, “The correlated correspondence algorithm for unsupervised registration of nonrigid surfaces,” in Proc. NIPS, 2005, pp. 33–40.
 [14] A. M. Bronstein, M. M. Bronstein, and R. Kimmel, “Generalized multidimensional scaling: a framework for isometryinvariant partial surface matching,” PNAS, vol. 103, no. 5, pp. 1168–1172, 2006.
 [15] Q.X. Huang, B. Adams, M. Wicke, and L. J. Guibas, “Nonrigid registration under isometric deformations,” CGF, vol. 27, no. 5, pp. 1449–1457, 2008.
 [16] A. Tevs, M. Bokeloh, M. Wand, A. Schilling, and H.P. Seidel, “Isometric registration of ambiguous and partial data,” in Proc. CVPR, 2009, pp. 1185–1192.
 [17] J. Sun, M. Ovsjanikov, and L. Guibas, “A concise and provably informative multiscale signature based on heat diffusion,” CGF, vol. 28, no. 5, pp. 1383–1392, 2009.
 [18] V. G. Kim, Y. Lipman, and T. Funkhouser, “Blended intrinsic maps,” TOG, vol. 30, no. 4, pp. 79:1–79:12, 2011.
 [19] M. Ovsjanikov, M. BenChen, J. Solomon, A. Butscher, and L. Guibas, “Functional maps: a flexible representation of maps between shapes,” TOG, vol. 31, no. 4, pp. 30:1–30:11, 2012.
 [20] A. M. Bronstein, M. M. Bronstein, R. Kimmel, M. Mahmoudi, and G. Sapiro, “A GromovHausdorff framework with diffusion geometry for topologicallyrobust nonrigid shape matching,” IJCV, vol. 89, no. 23, pp. 266–286, 2010.
 [21] Y. Lipman, R. M. Rustamov, and T. A. Funkhouser, “Biharmonic distance,” TOG, vol. 29, no. 3, pp. 27:1–27:11, 2010.
 [22] K. Xu, H. Zhang, W. Jiang, R. Dyer, Z. Cheng, L. Liu, and B. Chen, “Multiscale partial intrinsic symmetry detection,” TOG, vol. 31, no. 6, pp. 181:1–181:11, 2012.
 [23] M. Berger, A Panoramic View of Riemannian Geometry. Springer, 2002.
 [24] W. Zhang, X. Wang, D. Zhao, and X. Tang, “Graph degree linkage: Agglomerative clustering on a directed graph,” in Proc. ECCV, 2012, pp. 428–441.
 [25] O. van Kaick, H. Zhang, G. Hamarneh, and D. CohenOr, “A survey on shape correspondence,” CGF, vol. 30, no. 6, pp. 1681–1707, 2011.
 [26] G. Tam, Z.Q. Cheng, Y.K. Lai, F. Langbein, Y. Liu, D. Marshall, R. Martin, X.F. Sun, and P. Rosin, “Registration of 3d point clouds and meshes: A survey from rigid to nonrigid,” TVCG, vol. 19, no. 7, pp. 1199–1217, 2013.
 [27] A. Elad and R. Kimmel, “On bending invariant signatures for surfaces,” TPAMI, vol. 25, no. 10, pp. 1285–1295, 2003.
 [28] V. Jain, H. Zhang, and O. van Kaick, “Nonrigid spectral correspondence of triangle meshes,” IJSM, vol. 13, no. 1, pp. 101–124, 2007.
 [29] S. Wuhrer, Z. B. Azouz, C. Shu, and P. Bose, “Posture invariant correspondence of incomplete triangular manifolds,” IJSM, vol. 13, no. 2, pp. 139–157, 2007.
 [30] O. van Kaick, H. Zhang, and G. Hamarneh, “Bilateral maps for partial matching,” CGF, vol. 32, no. 6, pp. 189–200, 2013.
 [31] M. Aubry, U. Schlickewei, and D. Cremers, “The wave kernel signature: A quantum mechanical approach to shape analysis,” in Proc. ICCV Workshops, 2011, pp. 1626–1633.
 [32] A. Sharma, R. P. Horaud, J. Cech, and E. Boyer, “Topologicallyrobust 3d shape matching based on diffusion geometry and seed growing,” in Proc. CVPR, 2011, pp. 2481–2488.
 [33] Y. Sahillioglu and Y. Yemez, “Scale normalization for isometric shape matching,” Computer Graphics Forum (Proc. Pacific Graphics), vol. 31, no. 7, 2012.
 [34] A. Bronstein, M. Bronstein, A. Bruckstein, and R. Kimmel, “Partial similarity of objects, or how to compare a centaur to a horse,” IJCV, vol. 84, no. 2, pp. 163–183, 2009.
 [35] D. Raviv, A. Bronstein, M. Bronstein, and R. Kimmel, “Full and partial symmetries of nonrigid shapes,” IJCV, vol. 89, no. 1, pp. 18–39, 2010.
 [36] D. Hähnel, S. Thrun, and W. Burgard, “An extension of the ICP algorithm for modeling nonrigid objects with mobile robots,” in Proc. IJCAI, 2003, pp. 915–920.
 [37] H. Zhang, A. Sheffer, D. CohenOr, Q. Zhou, O. van Kaick, and A. Tagliasacchi, “Deformationdriven shape correspondence,” CGF, vol. 27, no. 5, pp. 1431–1439, 2008.
 [38] H. Li, R. W. Sumner, and M. Pauly, “Global correspondence optimization for nonrigid registration of depth scans,” CGF, vol. 27, no. 5, pp. 1421–1430, 2008.
 [39] M. Wand, B. Adams, M. Ovsjanikov, A. Berner, M. Bokeloh, P. Jenke, L. Guibas, H.P. Seidel, and A. Schilling, “Efficient reconstruction of nonrigid shape and motion from realtime 3d scanner data,” TOG, vol. 28, no. 2, pp. 1–15, 2009.
 [40] H. Li, B. Adams, L. J. Guibas, and M. Pauly, “Robust singleview geometry and motion reconstruction,” TOG, vol. 28, no. 5, pp. 175:1–175:10, 2009.
 [41] A. Pentland and J. Williams, “Good vibrations: modal dynamics for graphics and animation,” in Proc. SIGGRAPH, 1989, pp. 215–222.
 [42] R. Schmidt, C. Grimm, and B. Wyvill, “Interactive decal compositing with discrete exponential maps,” TOG, vol. 25, no. 3, pp. 605–613, 2006.
 [43] R. Schmidt, “Stroke parametrization,” CGF, vol. 32, no. 3, 2013.
 [44] E. L. Malvaer and M. Reimers, “Geodesic polar coordinates on polugonal meshes,” CGF, vol. 31, no. 8, pp. 2423–2435, 2012.
 [45] R. Palais, “On the differentiability of isometries,” Proc. of the AMS, vol. 8, no. 4, pp. 805–807, 1957.
 [46] A. Öztireli, G. Guennebaud, and M. Gross, “Feature preserving point set surfaces based on nonlinear kernel regression,” CGF, vol. 28, no. 2, 2009.
 [47] F. Mémoli and G. Sapiro, “A theoretical and computational framework for isometry invariant recognition of point cloud data,” FoCM, vol. 5, no. 3, pp. 313–347, 2005.
 [48] Y. Sun and M. Abidi, “Surface matching by 3d point’s fingerprint,” in Proc. ICCV, 2001, pp. 263–269.
 [49] R. Rustamov, “Barycentric coordinates on surfaces,” CGF, vol. 29, no. 5, 2010.
 [50] A. Tevs, A. Berner, M. Wand, I. Ihrke, M. Bokeloh, J. Kerber, and H.P. Seidel, “Animation cartography  intrinsic reconstruction of shape and motion,” TOG, vol. 31, no. 2, p. 15, April 2012.
 [51] D. Vlasic, I. Baran, W. Matusik, and J. Popović, “Articulated mesh animation from multiview silhouettes,” TOG, vol. 27, no. 3, pp. 97:1–97:9, 2008.
 [52] D. Anguelov, P. Srinivasan, D. Koller, S. Thrun, J. Rodgers, and J. Davis, “Scape: Shape completion and animation of people,” TOG, vol. 24, no. 3, pp. 408–416, 2005.
 [53] L. Yin, X. Wei, Y. Sun, J. Wang, and M. J. Rosato, “A 3d facial expression database for facial behavior research,” in Proc. FG, 2006, pp. 211–216.
 [54] J. Starck and A. Hilton, “Surface capture for performance based animation,” CG&A, vol. 27, no. 3, pp. 21–31, 2007.
 [55] M. BenChen, A. Butscher, J. Solomon, and L. Guibas, “On discrete killing vector fields and patterns on surfaces,” CGF, vol. 29, no. 5, 2010.