Shape-Aware Matching of Implicit Surfaces Based on Thin Shell Energies

A shape sensitive, variational approach for the matching of surfaces considered as thin elastic shells is investigated. The elasticity functional to be minimized takes into account two different types of nonlinear energies: a membrane energy measuring the rate of tangential distortion when deforming the reference shell into the template shell, and a bending energy measuring the bending under the deformation in terms of the change of the shape operators from the undeformed into the deformed configuration. The variational method applies to surfaces described as level sets. It is mathematically well-posed, and an existence proof of an optimal matching deformation is given. The variational model is implemented using a finite element discretization combined with a narrow band approach on an efficient hierarchical grid structure. For the optimization, a regularized nonlinear conjugate gradient scheme and a cascadic multilevel strategy are used. The features of the proposed approach are studied for synthetic test cases and a collection of geometry processing examples.


INTRODUCTION
We present a variational model for the matching of surfaces implicitly represented as level sets. The approach is inspired by the mathematical theory of nonlinear elasticity of thin shells. The model consists in an energy functional, which is to be minimized among deformations of a computational domain in which two given surfaces are embedded. A minimizer of this functional is a deformation that closely maps one (reference) surface onto the other (template) surface. As the underlying model we consider the reference surface as a thin elastic shell, i.e., a layer of an elastic material embedded in a volume of another several orders of magnitude softer isotropic elastic material. Subject to matching forces the volume is deformed in such a way that the thin shell is mapped onto the template surface. The functional reflects desired phenomena like resistance to compression and expansion of the surface, resistance to bending, and rotational invariance, while solely involving the deformation and the Jacobian of the deformation. The model is formulated in terms of projected derivatives from the tangent space of the reference surface onto the expected tangent space of the template surface. Taking into account a suitable factorization of the natural pullback under a deformation of shape operators enables us to formulate a model with appropriate convexity properties. The actual surface matching constraint is handled through a penalty, allowing for efficient numerical computation.
Through arguments of compensated compactness, we are able to show weak lower semicontinuity of the energy and consequently existence of minimizing deformations. We present a numerical approach based on a multilinear finite element ansatz for the deformation implemented on adaptive octree grids. The resulting discrete energy is minimized in a multiscale fashion applying a regularized gradient descent.
In the conference article [32] a preliminary version of this approach was presented. For the functional in that paper lower semicontinuity could not be ensured for either the membrane or bending energies. This lack of lower semicontinuity manifests itself in applications, where compression of the surface is expected, and leads to undesired oscillations in almost-minimizing deformations, which we explore in the present work through explicit examples and computations. Additionally, to increase the efficiency the computational meshes are in the present paper adapted to the surfaces. Consequently the number of degrees of freedom scales asymptotically almost like that of a surface problem.
The main pillar of our modelling is the use of polyconvex energy densities, first introduced in [3]. Energies of this type allow for geometric consistency properties like rotation invariance and the ability to measure area and volume changes. The core insight of this theory is that integrands consisting of convex functions of subdeterminants of the Jacobian give rise to integral functionals that are weakly lower semicontinuous in suitable Sobolev spaces. Indeed, this can be seen as an instance of compensated compactness [50]. A generic polyconvex isotropic energy density of the type used in this work is for Cof A:= det A A −T the cofactor matrix of A. Here, the coefficients and the function Γ are such that (1.1) attains a local minimum for A ∈ SO(n), that is, for rigid motions. Such an example is provided below in (3.8). Often in the modelling of nonlinear elasticity the condition is added, to reflect the non-interpenetration of matter [4]. In our model we make use of densities both with and without this property. Related work. Linear elasticity has been extensively used in computer vision and in graphics. Prominent applications are image registration [48,38,55,33,34], optical flow extraction [35], and shape modeling [29]. Recently, theories of nonlinear elasticity have been applied in many computer vision and graphics problems such as mesh deformation [13], shape averaging [56], registration of medical images [12]. The advantage of nonlinear models is that they allow for intuitive deformations when the displacements are large.
In this paper, we present a model for nonlinear elastic matching of thin shells. A finite element method for the discretization of bending energies of biological membranes has been introduced in [5]. Their approach uses quadratic isoparametric finite elements to approximate the interface on which the gradient flow of an elastic energy of Helfrich type is considered. The papers [9,10] discuss accurate convex relaxation of higher order variational problems on curves described as jump sets of functions of bounded variation. In particular, it enables the numerical treatment of elastic energies on such curves.
One challenge in polyhedral surface processing is to provide consistent notions of curvatures and second fundamental forms, i.e., notions that converge (in an appropriate topology or in a measure theoretic sense) to their smooth counterparts, given a smooth limit surface. One computationally popular model for discretizing the second fundamental form is Grinspun's et al. discrete shells model [30]. Another efficient, and robust method for nonlinear surface deformation and shape matching is PriMo [6]. This approach is based on replacing the triangles of a polyhedral surfaces by thin prisms. During a deformation, these prisms are required to stay rigid, while nonlinear elastic forces are acting between neighboring prisms to account for bending, twisting, and stretching of the surface. We refer to Botsch and Sorkine [7] for a discussion of pros and cons for various such methods. In comparison with methods based on polyhedral surfaces, level set approaches like ours are not dependent on specific triangulations of the shapes.
The matching of surfaces with elastic energies has recently been studied in [61]. Their energy contains a membrane energy depending on the Cauchy-Green strain tensor and a bending-type energy comparing the mean curvatures on the surfaces. The matching problem is formulated in terms of a binary linear program in the product space of sets of surface patches. For computations, a relaxation approach is used.
A different direction is the use of parametric approaches to reduce shape matching problems to the matching of functions on a fixed domain. For example, the methods presented in [62] and [60] are based on conformal maps from the unit disk. A more general variant using conformal maps on surfaces with arbitrary topology is presented in [42]. Within the family of parametric methods, a surface matching approach related to ours is presented in [44], where nonlinear elastic energies are used for matching parametrized surface patches. In comparison to all these methods, our level set approach is non-parametric and allows surfaces of any topology, which does not need to be fixed in advance.
In [59], face matching based on a matching of corresponding level set curves on the facial surfaces is investigated. To match pairs of curves an optimal deformation between them is computed using an elastic shape analysis of curves. Compared to our approach, this model does not take into account bending dissipation of the curves.
A different direction in shape recognition and matching is exploiting the intrinsic geometry of the surfaces only, thereby producing isometry-invariant methods based on the first fundamental form, like those in [24,11]. In comparison, bending is penalized in our model and we use all curvatures of the surfaces and their directions to be able to better match regions of edges and creases correctly.
A method for matching and blending of curves represented by level sets has been presented in [49]. Thereby, a level set evolution generates an interpolating family of curves, where the associated propagation speed of the level sets depends on differences of level set curvatures. In this class of approaches, geometric evolution problems are formulated, whereas here we focus on variational models for matching deformations. Variational registration of implicit surfaces was also considered in [40], but only through volume elasticity, in contrast to our shell terms.
To summarize, the main novelty of our contribution is the combination of independence of mesh topologies arising from the use of level sets, penalization of tangential distortion in a rotationallyinvariant framework, and awareness both of curvatures and curvature directions of the surfaces in the matching. We are not aware of any other methods possessing all of these features simultaneously.
Our approach is inspired by the articles [21,22] in which surface PDE models are derived in terms of the signed distance function. Shape warping based on the framework of [21] has been discussed from a geometric perspective in [14].
Outline. The paper is organized as follows. In Section 2, we review the required preliminaries about distance functions and formulate the geometric non-distortion and matching conditions that inspire our model. In Section 3, we present the different contributions to our energy. Section 4 is devoted to proving the existence of minimizing deformations under suitable Dirichlet and Neumann boundary conditions. Furthermore, the strong convergence of solutions for vanishing matching penalty parameter is discussed and counterexamples showing the lack of lower semicontinuity of related simpler models are given. In Section 5 a numerical strategy for minimizing the energy on adaptive octree grids is presented. Finally, Section 6 contains a range of numerical examples demonstrating the behaviour of solutions corresponding to our design criteria, and presents several potential applications. Some useful notation. For later usage and the purpose of reference let us collect some useful notation, mostly introduced in detail in later sections: • |B| stands for the Lebesgue measure of B ⊂ R n , and diam B = sup x,y∈B |x − y| for its diameter. • Generic matrices are denoted by A, B, M, N . We use 1 for the identity matrix. The set of rotations is denoted by O(n) and SO(n) is the set of orientation-preserving rotations. The set of all symmetric and positive definite matrices is SPD(n). • Components of vectors are denoted with subindices. For v ∈ R n , |v| denotes its Euclidean norm. The (n − 1)-dimensional sphere is S n−1 . For a matrix M , |M | is the Frobenius norm. • For two column vectors v, w ∈ R n , v ⊗ w is the tensor product of v and w, that is, the square matrix vw T . In particular, if |w| = 1 we have the identity (v ⊗ w)w = v. • P(e) = 1 − e ⊗ e is the projection onto vectors orthogonal to e ∈ S n−1 .
• Deformations on R n are denoted by φ, and deformations defined on a hypersurface M ⊂ R n by ϕ. The identity deformation is denoted by id.
• Ω ⊂ R n denotes the computational domain. Every relevant deformation φ maps Ω into R n . Ω has to contain all computationally relevant manifolds M. Ω has Lipschitz boundary, is open and bounded. • We use the notation ∂ i for partial derivatives, ∇ for the gradient of a scalar function, D for the Jacobian matrix of a vector function and D 2 for the Hessian matrix of a scalar function. • M 1 , M 2 ⊂ Ω are C 2,1 compact hypersurfaces. The inside and outside components of Ω \ M i are well defined by the Jordan-Brouwer separation theorem ( [31], Chapter 2, Section 5).
The signed distance function to M 1 , M 2 is denoted by d 1 , d 2 . The sign convention is that where the distance functions dist(·, M i ), i = 1, 2 are the unique viscosity solutions of 1 − |∇ dist(·, M i )| = 0 and dist(·, M i ) = 0 on M i . The normal fields to the offsets of M i at a point x are denoted by n i (x) := ∇d i (x). A superscript next to M i (i = 1, 2), as in M c i , denotes that we are talking about a level set of d i with value different from zero, so that M c i : is given by n i (x), and the set of points where d i is not differentiable is denoted by sing d i .
We use S i = D 2 d i for the Hessian of d i , which coincides with an extended shape operator of M i . • λ, µ are the Lamé coefficients of an isotropic material in linearized elasticity.
• C 0 (Ω; R n ) is the space of continuous functions from the domain Ω to the range R n , C k,α the Hölder spaces in which the k-th derivative is α-Hölder continuous, including the Lipschitz case α = 1. The range of the spaces is specified unless it is R. Sobolev spaces are denoted by W 1,p and the closure of compactly supported smooth functions in them by W 1,p 0 . • The letter C is reserved for a generic positive constant that may have different values in each appearance. Sequence indexing is usually denoted by a superscript k, and limits by an overline, e.g., φ k → φ.

DEFORMATION AND MATCHING OF LEVEL SET HYPERSURFACES
We are given two compact, connected embedded hypersurfaces M 1 , M 2 of class C 2,1 , which are diffeomorphic to each other, and both of which are contained in a bounded Lipschitz domain Ω ⊂ R n . In this section we deal with the tangential distortion and the change of the shape operator under a deformation φ : Ω → R n .
For any c ∈ R, we denote the c-offsets to the hypersurface M i by M c i := {x ∈ Ω | d i (x) = c} . Furthermore, we define the singularity set sing d i as the set of points where d i is not twice differentiable. With the regularity of M i that we have assumed, it is well known (e.g., Theorem 1. The gradient of the signed distance function ∇d i (x) is the outward-pointing unit normal n i (x) to M , consists of all vectors orthogonal to n i (x). Then, the corresponding projection matrices onto the tangent spaces are defined by P i (x) := P(n i (x)) = 1 − n i (x) ⊗ n i (x).
at a point x. In fact, from |n i (x)| 2 = 1 we deduce by differentiation that n T i (x)S i (x) = 0. This, together with the fact that n i ⊗ n i is the projection onto the normal of the hypersurface S i shows that With our choice of signs for d i , the symmetric matrices S i are positive semidefinite for convex hypersurfaces M i . Further information on tangential calculus for level set functions may be found in Chapter 9 of [23].
2.1. Tangential derivative and area and length distortion. First, let us assume that φ exactly maps M c 1 onto M c 2 , for all c > 0. Then, = im P 2 (φ(x)) and we define the tangential derivative induced by the deformation φ as capturing the tangential variation of φ(x) on M 2 along tangential directions on M 1 . In the variational model we consider below an energy term depending on D tg φ(x) will reflect the tangential distortion of the deformation in the context of a matching of the two hypersurfaces M 1 and M 2 even though φ(M 1 ) does not necessarily equal M 2 . Indeed, in the case M 2 = φ(M 1 ) the variation along a tangent direction on M 1 is still projected via D tg φ(x) onto the tangent space and not onto the tangent space of the deformed hypersurface φ(M 1 ) (cf. Fig. 1). Therefore there may exist tangential , such that D tg φ(x)v = 0 even though Dφv = 0. Thus D tg φ(x) can only be considered a measure of tangential distortion if φ(M 1 ) is sufficiently close to M 2 in the sense of closeness of tangent bundles. For a general deformation ψ : R n → R n the Cauchy-Green strain tensor Dψ T Dψ describes (up to first order) the deformation in a frame invariant (with respect to rigid body motions) way. Since we are interested in the effect of such a deformation between two hypersurfaces, for a suitably extended tangential gradient D tg φ + n 2 • φ ⊗ n 1 we define the extended tangential part of the Cauchy-Green strain tensor, measuring only tangential distortion: The term n 2 (φ(x)) ⊗ n 1 (x) is used to complement directions that are removed by the projections in the definition of the tangential distortion D tg φ and can be seen to realize a nonlinear Kirchhoff-Love assumption [17,Page 336], which postulates that lines normal to the middle surface of a shell remain normal after the deformation, without stretching. Next, we investigate the area and length distortion due to the tangential derivative D tg φ. For a given vector e ∈ R n we denote by Q(e) any proper rotation such that Q(e)e n = e, where e n denotes the n-th element of the canonical basis of R n . Note that this condition does not specify a unique Q(e). Then, for every B ∈ R n×n satisfying w ∈ ker B and im B ⊆ v ⊥ for some unit vectors v, w ∈ S n−1 , we have Hence, for φ(M 1 ) = M 2 and v = n 2 (φ(x)), w = n 1 (x) the area distortion under the hypersurface matching deformation φ at some position x is described by det(D tg φ(x) + n 2 (φ(x)) ⊗ n 1 (x)), which equals the positive square root of the determinant of the above Cauchy-Green strain tensor D tg φ T D tg φ + n 1 ⊗ n 1 . The squared tangential length distortion (in the sense of summing all squared distortions with respect to an orthogonal basis) is described by |D tg φ(x) + n 2 (φ(x)) ⊗ n 1 (x))| 2 and equals the trace of the Cauchy-Green strain tensor.

2.2.
Bending and curvature mismatch. Now, we quantify the change of curvature directions and magnitudes under the deformation φ. Our approach is motivated by models describing bending of elastic shells, because in our application the hypersurfaces are considered as thin shells.
In order to quantify the changes of curvature we first assume that φ(M 1 ) = M 2 , and compute the difference of the pull back of the shape operator S 2 on M 2 onto M 1 under the deformation φ and the shape operator S 1 on M 1 , which, for two arbitrary directions v, w ∈ R n , is given by If v, w are tangent vectors in T x M 1 , this difference describes the relative shape operator.
We define the extended relative shape operator For n = 3 and when φ is an isometric deformation between M 1 and M 2 (that is, Dφ(x) is an orthogonal mapping on T x M 1 for all x ∈ M 1 ), S rel appears in physical models for thin elastic shells in the context of the Γ−limit of 3D hyperelasticity [28]. Even though we do not necessarily expect our deformations to be tangentially isometric, we use this ansatz to compare curvatures of level sets in deformed and undeformed configuration, respectively. The following calculations shed some light on the properties of S rel : The assumption that φ(M 1 ) = M 2 can be rewritten as is again a distance function, and since d 2 • φ = 0 it follows that the left hand side of (2.5) is the shape operator of the hypersurface M 1 . The first term in the right hand side is the pullback of S 2 .
Let us remark that the appearance of a second fundamental form is consistent with Koiter's nonlinear thin shell theory [36], [17, Section 11.1]. Regardless of whether d 2 • φ is a distance function or not, (2.5) implies that In the next section, we use the extended relative shape operator to derive a variational model for the mismatch of curvatures.

ENERGY FUNCTIONAL
Given two hypersurfaces M 1 and M 2 our ultimate goal is to describe best matching deformations φ, which map M 1 onto M 2 as the minimizer of a suitable energy. Thereby, different energy terms will reflect a set of matching conditions for a volumetric deformation φ : Ω → R n and without a hard constraint φ(M 1 ) = M 2 : • A membrane deformation energy E mem penalizes the tangential distortion measured through D tg φ. • A bending energy E bend penalizes bending as reflected by the relative shape operator.
• A matching penalty E match ensures a proper matching of the two hypersurfaces M 1 onto M 2 via a narrow band approach. • A volume energy E vol enforces a regular deformation on the whole computational domain Ω. Our approach is based on level sets. Hence, we replace the integration over a single hypersurface, i.e., M 1 , for the first three energies by a weighted integration over a narrow band of width σ with 0 < σ < dist(M 1 , sing d 1 ). To this end we will make use of a cutoff function η σ ∈ C ∞ 0 (R) with R η σ (t)dt = 1 and supp η σ = [−σ, σ]. Additionally, η σ is assumed to be even and strictly decreasing in [0, +∞).
In what follows we introduce the four energy contributions separately.
3.1. Tangential distortion energy. Picking up the insight gained in Section 2.1 we formulate the membrane energy in terms of the length and area change associated with the tangential distortion D tg φ: where W is a nonnegative polyconvex energy density vanishing at SO(n). The weight δ reflects the proper scaling of the tangential distortion energy in case of a thin shell model with shell thickness δ. The energy (3.1) vanishes only on deformations φ whose Jacobian matrix Dφ(x) maps T M In consequence, both tangential expansion and compression are penalized.
Let us remark, that the extension D tg φ+n 2 (φ(x))⊗n 1 (x) of the tangential derivative D tg φ defined in (2.1) with rank n−1 can degenerate or be orientation-reversing depending on the local configuration of M 1 and M 2 at x (cf. Figure 2 for examples). FIGURE 2. Configurations in which for the (obviously isometric) identity we have det D tg 1(x) + n 2 (1x) ⊗ n 1 (x) = 0 (left) and det D tg 1(x) + n 2 (1x) ⊗ n 1 (x) < 0 (right) and thus the extended tangential derivative degenerates or reverses orientation.
Furthermore, the energy density W should not satisfy W (B) → ∞ for det B → 0. A straightforward modification of the arguments of Ciarlet and Geymonat ( [18], [16] Theorem 4.10-2) leads to a smooth integrand W which has isometries as local minimizers, with the correct invariance properties, and with a Hessian for B = 1 which matches the quadratic energy integrand of the Lamé-Navier model of linearized elasticity. With given Lamé coefficients λ, µ > 0, we select the energy This density fits into the notation of (1.1), if we choose p = q = 2 and Γ(t) = ct 2 + de −(t−1) .

3.2.
Bending energy. Now, we discuss a variational formulation of the curvature matching condition , which is equivalent to a vanishing relative shape operator (cf. (2.4)), where S i = Dn i P i = P i D 2 d i P i for i = 1, 2 are the shape operators on the hypersurfaces M 1 and M 2 , respectively. At first sight, it appears natural to formulate a quadratic penalization and to define a bending energyẼ The weight δ 3 reflects the scaling of the bending energy for thin shells of thickness δ. However, this energy is in general not weakly lower semicontinuous. Indeed, consider a situation in which S 1 (x) = S 2 (φ(x)) = P (e n ). A short computation shows that the corresponding density is not convex in its matrix variable along the rank-one segment joining 1 − 3 4 e 1 ⊗ e 1 and 1 − 1 2 e 1 ⊗ e 1 , which precludes lower semicontinuity (cf. also Example 4.8 below on the lack of rank-one convexity). In fact, this kind of density is closely related to the Saint Venant-Kirchhoff energy, whose quasiconvex envelope is computed in [39].
Thus, we are asking for an alternative lower semicontinuous energy functional which gives preference to deformations φ for which Dφ T (S 2 • φ)Dφ is close to S 1 . We show that this can be achieved with the extended shape operators S ext i = P i D 2 d i P i + n i ⊗ n i for i = 1, 2 and factorization. For proving this we make use of the following lemma.
Moreover, assume that A ∈ R n×n satisfies then the following statements are equivalent: If we multiply this equation from left and right by P 1 M 1 2 P 1 and take into account that (P 2 N 1 2 P 2 ) 2 = P 2 N P 2 and (P 1 M 1 2 P 1 ) 2 = P 1 M P 1 we see that this is equivalent to Applying that AP 1 = P 2 A we finally achieve at the equivalent condition The proof of the converse follows the same steps in opposite direction.
If the assumptions of this lemma apply to M = S ext 1 (x), N = S ext 2 (y), and A = Dφ(x) with y = φ(x), then the curvature matching condition , Dφ(x)) ∈ O(n) and a lower semicontinuous energy functional penalizing deviations of Λ(S ext 1 (x), S ext 2 (φ(x)), Dφ(x)) from O(n) would be a proper choice for realizing curvature matching. Unfortunately, the positive definiteness assumption of Lemma 3.1 is not fulfilled if principal curvatures of M 1 or M 2 are negative. Hence, we are replacing the extended shape operator matrices S ext i by symmetric and positive definite curvature classification matrices C i = C(S ext i ), i = 1, 2, respectively. We have experimented with two different choices for C: where −µ is a strict lower bound of the principal curvatures. But in applications surfaces are frequently characterized by strong creases or rather sharp edges, leading to very large µ. As a consequence the relative difference of the eigenvalues is significantly reduced when dealing with the resulting curvature classification matrices. Thus, the variational approach is less sensitive to different principal curvatures of the input hypersurfaces.
• Another option is to use a truncation of the absolute value function for the eigenvalues of symmetric matrices. For a symmetric matrix B ∈ R n,n with eigenvalues λ 1 , . . . , λ n and a diagonalization B = Q T diag(λ 1 , . . . , λ n )Q we use the classification operator where |λ| τ = max{|λ|, τ } for some τ > 0. This approach properly represents the exact shape operator matching objective in case of principal curvatures of equal sign and absolute value larger than τ . A disadvantage of this construction is that it is not able to force the deformation to correctly match curvature directions on the hypersurface with the same absolute value of the principal curvatures but with different signs. That is, locally a saddle point of the hypersurface may be mistaken for an elliptical point. However, this effect is usually compensated globally, and in applications the ansatz performs well, in particular in matching regions of edges and creases (see Section 6).
Like for the membrane energy (3.1), if Dφ(x) is ensured to be orientation-preserving (det Dφ > 0) and n 1 · (n 2 • φ) > 0 (cf. Figure 2), the curvature matching condition is equivalent to . Based on these considerations, a suitable choice for the bending energy is where W can be chosen as the same polyconvex density already used for E mem .
3.3. Mismatch penalty and volumetric regularization energies. So far, we have defined tangential membrane and bending energies which quantify the appropriateness of deformations φ : Ω → R n in a narrow band around the hypersurface M 1 . In the derivation of these energies we assumed the constraint φ(M 1 ) = M 2 . However, such a constraint would be very hard to enforce numerically. Thus we use a weaker mismatch penalty instead: where 1/ν is a penalization parameter. Moreover, we aim for a regular deformation on the whole computational domain Ω which is globally injective. This, in particular, prevents from self-intersections of the deformed hypersurface φ(M 1 ). To achieve this we introduce the following volume regularization term based on a polyconvex densityŴ that enforces orientation preservation with p > n, q > n, s > (n − 1)q/(q − n), and with α p , β q , γ s > 0 ensuring that the densityŴ attains a local minimum when Dφ T Dφ = 1. As mentioned in the introduction, such an energy is weakly lower semicontinuous in W 1,p (Ω; R n ) when restricted to deformations whose Jacobian determinant is positive almost everywhere, and this condition is closed under weak convergence.
3.4. Total energy. Summing the above terms, our energy for shape-aware level set matching reads where the different terms depend on the fixed input geometries M 1 and M 2 through d 1 and d 2 .

EXISTENCE OF OPTIMAL MATCHING DEFORMATIONS
First we prove the following weak continuity lemma, which is a generalization of the classical result given in [50,Theorem 4.1]. Here the coefficients may depend on the deformed configuration. we have Proof. To prove (4.1) let ζ ∈ L p p−n (Ω). We show that Moreover, we denote Using the inequality (cf. [27,Theorem 4.7]) and Hölder's inequality it follows that Here, we have used that By the Rellich-Kondrakov embedding theorem ([1], Theorem 6.3 III) there exist subsequences of v k i , i = 1, 2, which for simplicity of notation are again denoted by v k i , i = 1, 2, that converge uniformly to v i , i = 1, 2, respectively. Taking into account the Lipschitz continuity estimate Next, we replace v i , i = 1, 2 inĪ k by a piecewise constant approximation on a grid superimposed to the computational domain Ω. Explicitly, we consider the finitely many non empty intersection ω z δ = δ(z + [0, 1] n ) ∩ Ω of cubical cells with Ω for z ∈ Z n and definē where z δ is any point inΩ ∩ ω z δ if this set is nonempty. Using analogous estimates as above we obtain where v 1,δ and v 2,δ are piecewise constant functions in Using the uniform continuity of v 2 and v 1 on Ω we obtain that I k δ − I k ≤ β(δ) for a monotonically increasing continuous function β : R + 0 → R with β(0) = 0. In particular the convergence is uniform with respect to k. The same argument applies for the difference of I and and we get Ī δ − I < Cβ(δ). Using (2.3) it follows that Thus det(P(v 2 (z δ ))AP(v 1 (z δ ))+v 2 (z δ )⊗v 1 (z δ )) = det(Ã) represents an (n−1)×(n−1) minor of the linear mapping corresponding to the matrix A with respect to different orthogonal basis in preimage space (associated with P(v 1 (z δ )) and v 1 (z δ )) and the image space (associated with P(v 2 (z δ )) and v 2 (z δ )). Indeed, denoting where we have used the orthogonal change of variables y = Q T 1 x and Cof nn denotes the minor obtained by erasing the last column and the last row. This change of orthogonal coordinates is fixed on each cell ω z δ . Since for each δ the domain Ω is covered by finitely many cells ω z δ , using the above computation and standard weak continuity results [20,Theorem 8.20] for determinants of minors of the Jacobian we obtain thatĪ k δ →Ī δ for k → ∞. Finally, for given we first choose δ small enough to ensure that Ī δ − I + Ī k δ −Ī k ≤ 2 . Then we choose k large enough to ensure that I k −Ī k + Ī k δ −Ī δ ≤ 2 . This proves that a subsequence of I k converges to I for k → ∞. Since the limit does not depend on the subsequence, we finally obtain weak convergence for the whole sequence.
To prove (4.2), consider the three sequences of matrix functions The determinant of the second expression above converges weakly as k → ∞ by the first part of the lemma, while the determinants of the first and third can be assumed to converge uniformly. Moreover, the matrices in (4.3) have the block structure shown in (2.3), so multiplying the three together and taking into account that P is a projection (depending on the argument) recovers the matrix appearing in the statement. Multiplicativity of the determinant and the fact that a product of strongly converging and one weakly converging sequence converges weakly then finishes the proof.
We are now in a position to prove existence of a minimizing deformation for the hypersurface matching energy E in a suitable set of admissible deformations. Of particular difficulty is that derivatives of d 2 are not defined in the whole of Ω and that in the functional these derivatives are evaluated at deformed positions. We handle this by ensuring that the involved deformations are such that terms involving these derivatives are not evaluated near the singularities. We obtain the following theorem: Assume further that where sing d i is the set of points where d i is not differentiable, and that C : R n×n → SPD(n) is continuous. Then there exists a constant 0 < ν 0 := ν 0 (Ω, M 1 , M 2 , σ, p, α p ) such that for 0 < ν ≤ ν 0 , the functional E ν has at least one minimizer φ among deformations in the space W 1,p 0 (Ω; R n ) + id. Moreover, φ is a homeomorphism of Ω into Ω, and φ −1 ∈ W 1,θ (Ω; R n ), where θ is given by θ = q(1 + s)/(q + s).
Proof. We proceed in several steps.
First, notice that with the assumption (4.4) we have Let φ k be a sequence of ρ-admissible deformations with E vol [φ k ] ≤ C. By (4.5) and using the Banach-Alaoglu and Rellich-Kondrakov theorems, a subsequence (again denoted by (φ k )) converges to a deformation φ, both in the W 1,p -weak and uniform topologies. Now, we have ([20, Theorem 8.20]) Additionally, since (4.7) holds and because E ν [φ k ] is bounded, Ω (det Dφ k ) −s dx is bounded by the definition ofŴ and det Dφ k ≥ 0 a.e. Together with (4.7), we have (4.8) det Dφ(x) > 0 a.e., so that φ is again ρ-admissible. Notice also that by a.e. positivity of the determinants, (4.7) and a standard lower semicontinuity result for convex integrands (see e.g. [20,Theorem 3.23]) implies and uniform convergence of φ k immediately leads to We claim that under the assumptions of this theorem, we also have that and (4.10) To see this, notice that φ k , φ being ρ-admissible ensures that the normal vectors satisfy Consequently, the first part of Lemma 4.1 (with V i = n i ) implies (4.11) with χ {|d 1 |≤σ} denoting the indicator function. Combining (4.11) with the polyconvexity of W , defining the function E mem , both introduced in (3.1) we find the assertion (4.9). Furthermore, by our assumptions on M i (see section 2), we have that Since C produces uniformly positive matrices, we have χ {|d 1 |≤σ} (C(S ext 1 )) −1 ∈ C 0 (Ω; R n×n ). We can then use a continuity result for square roots of nonnegative definite matrix-valued functions defined on Ω [15, Theorem 1.1] to see that The second part of Lemma 4.1 implies the weak convergence from which (4.10) follows by using the polyconvexity of W .
Step 3: Existence of minimizers restricted to admissible deformations. Since we have already seen that the set of ρ-admissible deformations is weakly closed and weakly compact, and that every term of E is weakly lower semicontinuous on this set, we just need to check that for all fixed ν > 0, the set of ρ-admissible deformations, with adequate ρ, is not empty. For some given σ satisfying dist(M 2 , sing d 2 ) − σ > 0 let ρ satisfy We construct a deformationφ, which is ρ-admissible and satisfies E ν [φ] < ∞. By assumption, there exists a diffeomorphism ϕ : M 1 → M 2 . Thus, we construct an extension of this diffeomorphism to {|d 1 | ≤ σ} along the normal directions using We can then extendφ to the inside and outside components Ω i , Ω o of Ω \ {|d 1 | ≤ σ} by solving the minimization problems for E vol with Dirichlet boundary conditions given by (4.13) on ∂Ω i and ∂Ω o \ ∂Ω, and byφ(x) = x on ∂Ω. For the resultingφ we have where the first two statements follow by construction, and the last two by virtue of ϕ being a diffeomorphism and the choice of σ. Moreover, we note that sinceφ has finite energy and the growth conditions assumed forŴ (see (3.7)), the condition det Dφ(x) > 0 for a.e. x is also satisfied [4].
Step 4: A priori estimate to remove the constraint. Next, we show that for any ρ satisfying (4.12) there exists a parameter ν 0 > 0 such that for all 0 < ν < ν 0 the constrained minimizers of E ν subject to (4.4) solves the unconstrained optimization problem, consisting in minimizing E ν on W 1,p 0 + id. To this end, we verify that every φ that satisfies (4.14) is ρ-admissible. It is immediate from (4.14) that E vol (φ) < +∞, and from the definition ofŴ in (3.7) it follows with the same arguments as in (4.8) that det φ > 0 a.e. We prove now that for all deformations φ satisfying (4.14) also satisfy This is sufficient because from (4.15) it follows for all x satisfying |d 1 (x)| ≤ σ by the triangle inequality that which is the third property of a ρ-admissible deformation φ.
To prove (4.15) we use the triangle inequality and estimate (4.16) By the monotonicity of η σ and the fact that the signed distance functions d i are Lipschitz continuous with constant 1 we have, for eachσ ∈ (0, σ) that Estimates (4.5) and (4.14) imply in turn Finally, combining (4.16), (4.17), and (4.18) we obtain Now we can apply Ehrling's lemma [54,Theorem 7.30] for the embeddings W 1,p (Ω) ⊂⊂ L ∞ (Ω) ⊂ L 2 (Ω) to control the last term in (4.19). Taking into account the Poincaré inequality and Dirichlet boundary conditions, we obtain for any > 0 a constant C( ) > 0 such that Now, for the first term in the right hand side of (4.20) we can estimate For the second term, denoting diam Ω = sup x,y∈Ω |x − y|, where we have applied the product rule, the definition of E match , η σ ∈ C ∞ 0 , η σ ≤ C, that |∇d i | = 1 a.e., i = 1, 2, the chain rule, and (4.14). The use of the chain rule is justified by [46,Theorem 2.2], since d 2 has Lipschitz constant 1.
Together, (4.20), (4.21), and (4.22) imply In light of (4.19) and (4.23), and since E ν [φ] is independent of ν, we can now choose firstσ, then and finally ν small enough to obtain Step 5: Injectivity. The injectivity and regularity of the inverse follow by the growth conditions satisfied by E vol and classical results of Ball [4,Theorems 2 and 3]. Note that Theorem 3 in [4] is stated in the mechanical application context in dimension n = 3, but it holds also in R n following the same proof and using the condition s > (n − 1)q/(q − n).
We have particularized the statement of Theorem 4.2 to the case of Dirichlet boundary conditions to ensure global invertibility. In fact, we also have existence of minimizing deformations for the case of Neumann boundary conditions. such that for 0 < ν ≤ ν N , the functional E ν possesses at least one minimizer among deformations in the space W 1,p (Ω; R n ).
Proof. The proof follows the same arguments used for Theorem 4.2, so we only point out the necessary modifications. We need a replacement for the coercivity estimate (4.5) and claim To verify this let us consider ω := {|d 1 | ≤ σ/2}. An adequate Poincaré inequality (see e.g. [41,Theorem 12.23]) implies that and we estimate the second term in the right hand side by where Hölder's inequality has been used to compare L 1 and L 2 norms. Therefore, (4.24) follows. The proof of the estimate for d 2 • φ L ∞ ({|d 1 |≤σ}) (to ensure that deformations stay away from the singularities of d 2 ) is still valid with minor modifications, since ν appears in (4.24) multiplicatively.
We conclude this section with the following proposition, which explores the penalization limit in which the parameter ν tends to zero. Proposition 4.4. Let {ν k } k∈N , be a sequence of penalty matching parameters such that ν k → 0 as k → ∞, and φ k be solutions of the Dirichlet minimization problem for E ν k . Then, up to a choice of subsequence, the φ k converge strongly in W 1,p to a minimizer of E mem + E bend + E vol in W 1,p 0 (Ω; R n ) + 1 under the constraint φ(M c 1 ) = M c 2 for all c ∈ (−σ, σ). Proof. First, notice that the energy E may be written as where H : R + × R n×n × R n×n × R × R n×n × R → R + is smooth and convex. Denote byφ the extension of a diffeomorphism between M 1 and M 2 used in the proof of Theorem . By the coercivity estimate (4.5) the φ k are then bounded in W 1,p and we may extract a (not relabelled) subsequence converging uniformly and weakly in W 1,p to some limit φ. Since {E ν k [φ k ]} is bounded and ν k → 0, the uniform convergence of φ k implies that is the uniform limit of the maps φ k M c 1 which are surjective onto M c 2 and M c 1 is compact, we conclude that φ(M c 1 ) = M c 2 for all c ∈ (−σ, σ). Therefore, φ is admissible for all ν k and E ν k [φ k ] ≤ E 1 [φ]. Combined with lower semicontinuity and (4.25), the above implies From this identity, the fact that H is convex and differentiable, and Dφ k Dφ in L p it follows that Together with the weak lower semicontinuity of the L p -norm, the above shows that Because L p (Ω) has the Radon-Riesz property ([47, 2.5.26]), weak convergence and convergence of the norm guarantee strong convergence. Since φ k was assumed to converge uniformly, we have also φ k → φ in L p , and this shows that φ k → φ in W 1,p (Ω; R n ).
That φ is a minimizer of the constrained problem follows directly ( [8], Theorem 1.21) from the fact that the E ν k are an equicoercive family of functionals, Γ-converging in the weak topology of W 1,p . Indeed, equicoercivity follows easily from the above, while Γ-convergence is implied by the fact that E ν k is an increasing sequence ( [8], Remark 1.40), because ν k → 0 appears as a denominator in E match .  Contrary to what might be expected, the limit problem we have obtained is not a surface problem, since all the level sets are still coupled through the volume energy E vol . The line of reasoning above depends heavily on the fact that the coefficients of the volume term are held fixed, since the equicoercivity and uniform strict quasiconvexity (in the language of [25]) both require the presence of Dφ p L p (Ω) in the functional. 4.1. Oscillations and lack of rank-one convexity for the naive approach. To model the tangential distortion energy we have considered a frame indifferent energy density with the argument D tg φ + (n 2 • φ) ⊗ n 1 . Let us now consider the case n = 2 and a simpler version of the membrane energy (3.6), where we use as an argument of the energy density directly the tangential Cauchy-Green strain tensor (cf (2.2)) (D tg φ(x)) T (D tg φ(x)) + n 1 (x) ⊗ n 1 (x), and define the membrane energy , and W : R 2×2 → R a frame indifferent energy density that has a strict minimum at SO (2). In fact, this energy is no longer lower semicontinuous and we will present counterexamples.
Example 4.7 (Oscillation patterns). We construct an explicit sequence for which lower semicontinuity of the membrane energyẼ mem fails. Fix 0 < R < 1 and M 1 = S 1 with the parametrization ξ → e iξ . Consider a sequence of deformations ϕ k : S 1 → R 2 defined in polar coordinates of (r, θ) by the condition (4.29) ∂ ξ ϕ k (ξ) = (R sin kξ) e r r(ϕ k (ξ)), θ(ϕ k (ξ)) + 1 − R 2 sin 2 kξ 1 2 e θ r(ϕ k (ξ)), θ(ϕ k (ξ)) , where e r = (cos θ, sin θ) T , e θ = (− sin θ, cos θ) T for given φ k (0). Note that for any k and θ that |∂ θ ϕ k (θ)| = 1, so that the transformations are tangentially isometric. We define ϕ k (0) via two integration constants r 0 and θ 0 for the initialization of r and θ at ξ = 0. We set θ 0 = 0 and choose r 0 such that the curve ϕ k is closed and simple, which imposes r 0 = r(ϕ k (0)) = r(ϕ k (2π)) since the first term in (4.29) has zero average. From the second term, taking into account that e θ (r, θ) is independent of r, we get the condition where we have applied the change of variables ζ = kξ. By periodicity the right hand side (an incomplete elliptic integral of the second kind with modulus R) is independent of k and thus determines r 0 . The resulting ϕ k for several values of k are depicted in Figure 3. We observe that ∂ θ ϕ k (θ) r 0 e θ in L p , for any 1 ≤ p < ∞ (and also weak-* in L ∞ ). Therefore, the weak W 1,p -limit ϕ of the ϕ k is the function defined by ϕ(θ) = r 0 e r and obviously not an isometry. Assuming 0 < σ < 1 and extending ϕ k , ϕ along the radial direction e r to the annulus {1 − σ ≤ r ≤ 1 + σ}, we obtain corresponding deformations given by where Q π 2 stands for clockwise rotation by π/2, so that Q π 2 ∂ θ ϕ k (θ) is the unit outward normal to ϕ k (S 1 ). Clearly also φ k φ in W 1,p on the annulus. We observe thatẼ mem [φ k ] = 0, butẼ mem [φ] > 0. Hence,Ẽ mem is not weakly lower semicontinuous. The celebrated Nash-Kuiper theorem [51,37] states that it is possible to uniformly approximate any short C ∞ immersion by C 1 isometric ones. Our explicit oscillations around r 0 S 1 is just one example of this phenomenon. Notice that a bending term of the type E bend introduced in our model only compares the curvatures of M  [Lack of rank-one convexity] We present an additional example of a configuration for which the integrand of an energy of the typeẼ mem is not rank-one convex. Rank-one convexity of the complete energy density, i.e., , convexity in t ∈ R when composed with the function A + tB for any matrix A and any rank one matrix B, is known to be a necessary condition for quasiconvexity ([20], Theorem 5.3). Quasiconvexity, in turn, is necessary for weak lower semicontinuity of integral functionals in Sobolev spaces ( [20], Theorem 8.1 and Remark 8.2).

FINITE ELEMENT DISCRETIZATION BASED ON ADAPTIVE OCTREES
We adopt a 'discretize, then optimize' approach and consider a finite element approximation and optimize for the coefficients of the solution. Since the energy E ν is highly nonlinear and nonconvex, we use a cascadic multilevel minimization scheme in which the solution for one grid level is used as the initial data for the minimization on the next finer grid. We use adaptive refinement of the underlying meshes around the surfaces M 1 , M 2 ⊂ (0, 1) n for n = 2, 3 (Algorithm 1).
One of the main characteristics of our functional is the pervasive presence of coefficients depending on the deformed position φ(x). Indeed, this is how the functional takes into account the geometry of target surface, through the projection P 2 and shape operator S 2 . From an implementation perspective, however, this means that frequently discrete functions have to evaluated at deformed positions. Therefore, the ability to efficiently search the index of an element containing a given position is of paramount importance, so a hierarchical data structure that allows for efficient searching is needed. The model only contains first derivatives of the unknown deformation. Hence, multilinear finite elements already allow a conforming discretization. For these reasons we use multilinear FEM on octree grids. The grids used are such that all of the elements are either squares or cubes of side length h = 2 − , for an integer to which we refer as grid level of the element. In what follows let us detail the different ingredients of the algorithm.

Algorithm 1 Cascadic minimization scheme.
Starting grid: Uniform of level min , h = 2 − min φ ← 1 for ← min to max do Regenerate d 1 , d 2 on grid by aFMM. Compute n 1 , n 2 , S 1 , Refine the grid (h ← 2 − +1 ). end for return φ Multilinear Finite Elements on Octrees. We assume n = 3 for the presentation here. Using an adaptive octree grid based on cubic cells leads to hanging nodes (see Figure 4), nodes which are on the facet of a cell without being one of its vertices. Enforcing continuity of the finite element functions leads to constraints for function values on hanging nodes and these hanging nodes are not degrees of freedom. Additionally, to minimize the complexity of the required interpolation rules, the subdivision is propagated in such a way that the grid level of neighboring elements sharing a cell facet differs at most by one. Octrees and the access to degrees of freedom via hashtables. Even though the tree structure gives a natural hierarchical structure to the elements of the mesh, maintaining consistent linear indices for degrees of freedom, hanging nodes, and elements can be delicate. Consistent rules could be devised to maintain consistency with the element octree for a given mesh, but these would not be easy to update when the grid is refined. In order to keep track of vertex indices in a simple manner without sacrificing efficiency, hash maps ( [19], Chapter 11) are maintained to keep track of the indices of degrees of freedom, hanging nodes, and cells. The keys used in the hashmap are a combination of a level value and point coordinates as integer multiples of h = 2 − . These keys uniquely identify nodes or elements, with the convention that an element is identified with its lower-left-back corner. Whenever a query for a node or cell is made, there are two possible outcomes. If it is already contained in the corresponding hash table, a linear index for it can be retrieved. Otherwise, a new entry of the hash table is created and the node or cell is given the next unused index. Since we do not require coarsening of the mesh, this scheme guarantees a consistent linear set of indices with a computational cost for insertions and queries that is, on average, independent of the mesh size. Computing distance functions on octrees. In our model we have assumed that the distance functions to our surfaces are given. In practice, especially when using adaptive grids, we need to compute signed distance functions on such grids. This has been accomplished by a straightforward adaptation of the Fast Marching Method on cartesian grids [57] exploiting the fact that our grids still are subgrids of a regular cartesian grid. In the implemented variant hanging nodes are not taken into account for the propagation, their values being linearly interpolated to accommodate the constraints needed for conformality. The initialization for the distance computation has been performed starting from triangular meshes of the surfaces (for n = 3; for n = 2 two-bit segmentation of interior and exterior of the curves has been used). The signs of the distance functions have to be computed separately, by detecting which points of the grid are inside (resp. outside) the initial surface data. In our case, they have been computed with the provably correct algorithm given in [2]. Computation of the coefficients. The discretization for the unknown deformation φ, as already mentioned, is done by multilinear finite elements. However, the coefficients of our model include first and second derivatives of the signed distance functions d i , for the normal vectors n i and shape operators S i (i = 1, 2), respectively. The approximations are required to be robust, since they appear in the highest order terms of the model. For the normal vectors n i , we compute the L 2 projection of the finite element derivative of d i to recover the nodal values of a piecewise multilinear function, followed by a orthogonal projection to the unit sphere to restore the constraint |n i | = 1. In the case of the shape operators, our approach is to approximate the distance functions d i by a quadratic polynomial supported on a neighborhood of each point. Given a fixed integer neighborhood size r, for each non hanging node x k (i.e., the neighborhood B r (x k ) contains the r closest other degrees of freedom of the adaptive grid) the local quadratic polynomial p k is defined as the one minimizing the least-squares error which can be easily computed by inverting a small matrix. The Hessian of d i at the node x k is then approximated by the Hessian of p k .
For the computation of matrix square roots and their inverses, we have used the method described in [26], taking appropriate care to truncate almost-singular matrices, since the resulting square roots also appear inverted. Minimization strategy. For the minimization at each level, we have opted for a Fletcher-Reeves nonlinear conjugate gradient method ( [52], Section 5.2). The L 2 gradient of E ν , whose computation is involved but elementary, was implemented directly. The parameter α is progressively reduced when a further feasible descent step is not found, according to an Armijo line search ( [52], Section 3.1).

NUMERICAL RESULTS
All of our results have been computed on the unit cube Ω = [0, 1] 3 for the matching of surfaces in 3D, and the unit square [0, 1] 2 for the matching of contour curves in 2D. In practice, we have used homogeneous Neumann boundary conditions, since this allows to have relatively large shapes M i in comparison to the size of the domain Ω without creating excessive volume energies near the boundary (for the justification we refer to Corollary 4.3). However, if the boundary is not fixed, the deformed domain φ(Ω) is not necessarily contained in Ω, so evaluation of coefficients on deformed positions has to be appropriately handled numerically. We use a projection of outside position onto the boundary of Ω for sufficient large dist(M 2 , ∂Ω).
For the membrane and the bending energy we use the material parameters λ = µ = 1, corresponding to the density In the bending term, the shape operators have been regularized through the truncated absolute value function with τ = 1. Since we work on the unit cube, this corresponds to a comparatively large curvature radius. For the volume term, given that enforcing orientation preservation in a finite element framework is a far from straightforward, it is advantageous to work with the simplified version c vol Ω W (Dφ) dx.
We have run the minimization scheme of Algorithm 1 beginning from a uniform grid of level min = 2 or min = 3 with 9 3 = 729 nodes, and refined up to max = 8 for 3D examples. For 2D cases a reasonable range turned out to be min = 4, max = 10. The finest grids used for two of the examples below are depicted in Figure 4. The width of the narrow band was chosen proportional the finest resolution of the mesh (σ = 2h) since a small value of σ clearly produces inaccurate results when η σ is evaluated on coarse grids. However, the constraint Ω η σ = 1 ensures that the overall strength of the surface terms E match , E mem and E bend is not affected. The value of the penalty constraint ν was divided by 10 for each grid refinement, which is justified by Proposition 4.4. Furthermore, the volume weight c vol was also halved per level to allow for simultaneously higher initial regularization and close final matches. Note that this reduction is much slower than that of the matching parameter.
In all examples, we have used the identity as the initial deformation. It should be noted that although the energy is geometric by design, we are using a first-order descent method for its minimization. In consequence, an adequate rigid pre-alignment can be beneficial for intricate shapes. Figure 8 shows results for the matching of two different dolphin shapes. Our variational approach is highly nonlinear and non-convex. Thus, the numerical approximation of the globally optimal deformation depends on the initialization of the deformation. Figure 9 shows that the identity deformation as the initial deformation is advisable only if the expected optimal deformation is not too large. This is demonstrated by applying different rigid body motions to M 1 .
All figures have been produced by deforming the input data (polygonal curve or triangulated surface) via the resulting deformation φ. This is in contrast to deforming the grid and plotting the resulting extracted level sets (which effectively visualizes the inverse deformation), as commonly done in the registration literature, and also in [32]. Test case. First we present a simple test case to underline the qualitative properties of our model. Figure 5 shows a configuration in which a high amount of compression, combined with rotation, is required. Our model finds the intuitively correct deformation, but oscillations typical for the lack of lower semicontinuity of the underlying energy are induced when P 2 is not used in the membrane and bending terms. The bending term assists in matching the curvatures even if the deformation is not rigid. Note, however, that for the optimal match the curvature energy E bend is not expected to vanish, as can easily be seen from (2.6), (3.4) and the related discussion in Section 3. Shape matching applications. We now turn our attention to high resolution examples with real data. Figure 6 demonstrates the effect of the multilevel descent scheme, in which details are added progressively to avoid spurious local minima. In Figure 11 a high-resolution 2D example is presented. Figures  7, 8 and 10 show 3D examples in which the influence of the curvature matching is indispensable to obtain shape sensitive matching deformations. For these examples, the shell parameter δ was chosen quite high, since the curvature matching term E bend is a major driving force to obtain correct matching of geometric features. Table 1 lists the parameter values used, and run times for our implementation. We have split the timings between the highest-resolution level and the combined previous ones, since in many applications a very high level of detail might not be necessary, thereby significantly reducing the required computational effort.   TABLE 1. Parameters and running times on a workstation with a single Intel Xeon E5-1650 CPU (6 cores, 3.2Ghz). Our implementation splits the computation of the different terms of the energy and the corresponding derivatives in different threads (obtaining a speedup factor ≈ 2), but no further parallelization is used.   . From left to right: Initial shape of Figure 8 after undergoing a rotation of π/6, deformed shape after level 8 in the minimization (correct matching), after a rotation of π, and corresponding result (incorrect matching). Moderate changes in the initial alignment are handled correctly, while drastic ones are not.