Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Int J Comput Vis. Author manuscript; available in PMC 2010 April 22.
Published in final edited form as:
Int J Comput Vis. 2008 December 1; 80(3): 317–336.
doi:  10.1007/s11263-008-0141-9
PMCID: PMC2858418

Large Deformation Diffeomorphic Metric Curve Mapping


We present a matching criterion for curves and integrate it into the large deformation diffeomorphic metric mapping (LDDMM) scheme for computing an optimal transformation between two curves embedded in Euclidean space Rd. Curves are first represented as vector-valued measures, which incorporate both location and the first order geometric structure of the curves. Then, a Hilbert space structure is imposed on the measures to build the norm for quantifying the closeness between two curves. We describe a discretized version of this, in which discrete sequences of points along the curve are represented by vector-valued functionals. This gives a convenient and practical way to define a matching functional for curves. We derive and implement the curve matching in the large deformation framework and demonstrate mapping results of curves in R2 and R3. Behaviors of the curve mapping are discussed using 2D curves. The applications to shape classification is shown and experiments with 3D curves extracted from brain cortical surfaces are presented.

Keywords: Large deformation, Diffeomorphisms, Vector-valued measure, Curve matching

1 Introduction

Matching curves embedded in R2 and R3 is important in the fields of medical imaging and computer vision (Zhang 1994; Besl and McKay 1992; Feldmar and Ayache 1996). In medical imaging, the most striking gross morphological features of the cerebral cortex are the sulcal fissures on the cortical surface (Welker 1990). They are selected as the basis for structural analysis on the internal surface anatomy of the brain because they separate functionally distinct regions and provide a natural topographic partition of the anatomy. Even though major sulci consistently appear in all anatomies, they exhibit variability in size and configurations (Welker 1990). A quantitative study of sulcal anatomical structures requires developing mathematical models for characterizing shape and geometry of curves, and for accommodating the variability present across populations (Durrleman et al. 2007; Fillard et al. 2007; Thompson et al. 1996; Rettmann et al. 2002). Moreover, image boundary, contours, and skeletons are often used to represent shapes of objects in many applications of computer vision, such as tracking human body motion in motion tracking, finding correspondence between maps and terrain images, object recognition and classification. All of these formulate two problems: representation of curves and mathematical models of transformations.

In non-rigid registration problems, desired transformations are often constrained to be diffeomorphic—one to one (invertible) and smooth with smooth inverse so that connected sets remain connected, disjoint sets remain disjoint, smoothness of features such as curves is preserved, and coordinates are transformed consistently. Consequently, many researchers have paid great attention to such diffeomorphic transformations particularly in the field of Computational Anatomy for studying the geometric variation of human anatomy (e.g. Grenander and Miller 1998; Gee and Bajcsy 1999; Miller et al. 2002; Twining et al. 2002; Beg et al. 2005; Joshi et al. 2004; Avants and Gee 2004; Thompson and Toga 1996; Thompson et al. 2004).

Over the past several years we have been using large deformation diffeomorphic metric mapping (LDDMM) (Joshi and Miller 2000; Camion and Younes 2001; Beg et al. 2005) to map 3D volume coordinates in the brain as well as in the heart (Helm et al. 2006). This method not only provides a diffeomorphic correspondence [var phi] between shapes C and S, but defines a metric distance in a shape space as well. The basic diffeomorphic metric mapping approach taken for understanding structures of object shapes is to place the set of objects into a metric space. The idea is to model the mapping of one from the other via a dynamic flow of diffeomorphisms (time dependent deformation) t [set membership] [0, 1] of the ambient space Rd. Instead of working directly in the space of diffeomorphisms, we work with time dependent vector fields, υt : RdRd for t [set membership] [0, 1], which model the infinitesimal efforts of the flow. The flow [var phi]tυ is then defined via the differential equation

[partial differential][var phi]tυ[partial differential]t=υt([var phi]tυ),

with [var phi]0υ(x)=x, and the resulting diffeomorphism is given by the flow at time t = 1, [var phi]1υ. In order to ensure that the resulting solution to this equation is a diffeomorphism, υt must belong to a space, V, of regular vector fields (see Trouvé 1995; Dupuis et al. 1998 for specific requirements), for t [set membership] [0, 1] with 01υtVdt<. V is a Hilbert space with an associated kernel function kV, and norm ‖·‖V. This norm models the infinitesimal cost of the flow.

Assume that for any pair objects of C and S, there exists some vector fields υt such that the diffeomorphism [var phi]1υ transforms shape C to shape S, which we write [var phi]1υ·C=S. Then one can look at the optimal vector fields [upsilon macr]t satisfying [var phi]1υ^·C=S and minimizing the integral of the infinitesimal costs. The metric distance between shapes, ρ(C, S), is defined as the length of the geodesic path [var phi]1υ^·C for t [set membership] [0, 1] generated from C to S in the shape space. Such a metric between C and S has the form

ρ(C,S)=infυt(01υtV2dt)12such that[var phi]1υ·C=S.

In practice, we often define an inexact matching problem: find a diffeomorphism [var phiv with circumflex]t between two objects C and S as a minimizer of

JC,S([var phi]t)=˙γρ(C,[var phi]1·C)2+E([var phi]1·C,S).

Equivalently, we have [var phi]^t=[var phi]tυ^, where the vector fields [upsilon macr]t minimize the energy

JC,S((υt)t[set membership][0,1])=˙γ01υtV2dt+E([var phi]1υ·C,S),

and γ is a trade-off parameter. The first term in (3) is a regularization term to guarantee the smoothness of deformation fields and the second term is a matching functional to quantify the mismatching between the objects [var phi]1υ·C and S.

In this LDDMM framework, two matching techniques are commonly used that result in different models for E([var phi]1υ·C,S). The first one is volume image matching, which works directly with the raw 3D volume images, multi-valued vectors, as well as tensor matrices arising from DTI studies (Beg 2003; Beg et al. 2005; Cao et al. 2005a, 2005b), but is computationally intensive. The second one is landmark matching which requires a manual or semi-automated selection of landmarks on the images (Joshi and Miller 2000). In the case of brain images, however, landmarks are usually not just isolated points, but lie along sulcal or gyral curves extracted from cortical surfaces which are themselves segmented from the 3D images. In such a situation, the landmark matching tool is commonly used (Joshi et al. 2007; Bakircioglu et al. 1999). It is imperfect, however, because it first requires to select the same number of points on both source and target curves, and second it specifies arbitrary correspondences between specific points on these curves, giving artificial constraints to the registration process. To overcome this issue, some new methods were proposed, e.g. using level-sets representations of the curves (Leow et al. 2004). In Glaunès et al. (2004), an extension of the LDDMM approach in which shapes are considered as unlabelled sets of landmarks has been proposed. However, as geometric objects, curves or surfaces should not be treated as a sequences of points because the higher order information (tangent vectors or curvature) is discarded when reducing 1D or 2D manifolds into 0-dimensional point sets. In the case of surfaces, a representation via currents, or vector-valued measures, has recently been developed (Vaillant and Glaunès 2005). It allows to put into the metric both location and tangential information of the surfaces. Tools presented in Glaunès et al. (2004) and Vaillant and Glaunès (2005) actually derive from the same general framework that defines computable metrics on any m-dimensional submanifolds of Rd with md. The case of curves (m = 1) has been introduced shortly in Glaunès et al. (2006), Qiu et al. (2007). It is the purpose of this paper to present a full mathematical explanation and derivation of this representation of curves in the LDDMM framework.

This model for curve matching is a special case of a general framework for submanifolds, which could be called measure or current matching, and uses ideas from geometric measure theory and reproducing kernels while setting up into Grenander’s theory of deformable templates (Grenander and Miller 1998). Other examples had already been presented: surfaces in Vaillant and Glaunès (2005) and unlabelled point sets (which correspond to the 0-dimensional case) in Glaunès et al. (2004). In particular, the study and use of the metric distance defined in the shape space for shape analysis has received great attention in computer vision and medical imaging as partly shown here in Sect. 5. The common approach consists in defining a Riemannian metric on the tangent space to a curve, i.e. the space of infinitesimal transformations of a curve, and extending it to the whole space. Sharon and Mumford (2006) proposed a very elegant mathematical model for closed planar curves, which can not be generalized. Younes (1998) proposed to use an elastic metric that measures the cost of deformation of the ambient space, which corresponds to the LDDMM setting described previously and considered in the present paper. Recent works (Klassen et al. 2003; Mio and Srivastava 2004; Schmidt et al. 2006) have emphasized the efficiency of the shape space model or developed its mathematical theory (Michor and Mumford 2007). However these constructions require that a smooth transformation actually exists between any two curves in the space, which usually restricts the study to the subclass of closed, non crossing curves. Our approach is different in the sense that we do not claim to find exact correspondences between curves, but rather find a diffeomorphic transformation that match two curves up to an acceptable error. What we propose here is a way to compute this error, i.e. a model for evaluating the differences between shapes that remain after matching. It allows to compare curves that may present noisy parts, or be composed of different numbers of connected components, or to match a closed curve to an open curve, etc. Such situations are common in practical applications, which indicates that the curve mapping presented in this paper will be more suitable to applications in computer vision and medical imaging analysis.

This paper is organized as follows. We start with the geometric representation of curves that incorporates the geometry of the curves due to the sensitivity to both location and their tangent vectors (the first order local geometric structure) in Sect. 2.1. In Sect. 2.2, we describe that such a representation of curves as mathematical objects provide elements in a vector space equipped with a computable normed distance. In Sect. 3, “closeness” between two curves is then quantified by the norm-square distance between their associated representations, which serves as matching criterion for finding an optimal transformation from one curve to the other. In Sect. 3.2, we integrate this matching functional with a variational optimization problem of the LDDMM and compute its gradient with respect to a linear transformation of the velocity, momentum. Appendix B, demonstrates the other use of the matching functional in the derivation of the similarity transformation (rotation, translation, and isotropic scaling). In Sect. 5, we finally demonstrate the application of this method to curves in R2 and R3 and compare matched results and deformation fields between this method and the landmark matching method (Joshi and Miller 2000) in the LDDMM framework.

2 Curves as Vector-Valued Measures

We define a norm that quantifies the similarity between two curves by first introducing a geometric representation of the curve and then inducing its norm in a Hilbert space via a differential operator.

2.1 Geometric Curve Representation

We follow the approach presented (Vaillant and Glaunès 2005) for surfaces. However, we will not use the term “currents” in this paper, nor use any notion of exterior calculus. Let C [subset or is implied by] Rd be an open or closed curve. We assume that C is both piecewise C0 and piecewise C1, which means that it may be composed of several disconnected components, and may present a finite number of angular points. Let γC : [0, 1] → Rd, s [mapsto] γC(s) = [x1(s), …, xd (s)] be a parametrization of this curve. At almost any point on the curve, the tangent vector is given by γC(s)=[dx1ds,,dxdds]. The existence of these tangent vectors along the curve gives us a natural action of the curve on vector fields w : RdRd by the rule:

left angle bracketµC|wright angle bracket=˙01γC(s)·w(γC(s))ds,

where · denotes the dot product in Rd. Intuitively, left angle bracketµC|wright angle bracket is large when the direction of the vector field w(x) is close to the tangent vectors γC(s) along the curve. It is close to zero when w(x) is almost perpendicular to the curve, or when w(x) take small values (in norm) along the curve. Figure 1 demonstrates the evaluation of left angle bracketµC|wright angle bracket in several cases when w(x) is parallel to the tangents of C (first panel), perpendicular to C (middle panel), and the norm of w(x) is close to zero along the curve (last panel).

Fig. 1
Evaluation of a vector field (blue arrows), w(x), along a curve (in red), C. Left: Directions of vectors match those of the tangents to the curve. Middle: Directions are perpendicular to the curve. Right: Directions match, but norm of vectors are close ...

Hence both location and tangential information of the curve is integrated in this expression. Another important point here is that it is invariant under re-parametrization of the curve, so that if we have a new parameter t = ψ(s), with ψ increasing, it is simple to show that left angle bracketµC|wright angle bracket is unaffected. However using t = ψ(s) with ψ decreasing will give an opposite value for left angle bracketµC|wright angle bracket, which means that the representation is orientation dependent.

Equation (4) associates to every oriented curve C a unique representant µC which is a linear functional acting on vector fields of Rd. Equivalently, µC can be understood as a vector-valued Borel measure. Linear combinations of measures are defined with the rule left angle bracketαµC + βµS|wright angle bracket = α left angle bracketµC|wright angle bracket + β left angle bracketµS|wright angle bracket for any α, β [set membership] R and curves C, S. In general, these combinations are not associated to any curve, unless |α| = |β| = 1: for example, µC − µS corresponds to the concatenation of curve C and curve S with opposite orientation.

The advantage here is to embed the curves in a vector space (of vector-valued measures) on which it is possible to use functional analysis tools to define norms and compare curves via these norms. This will be developed in Sect. 2.2.

2.2 Norm on Vector-Valued Measure

We would like to define a norm that computes a geometrical distance between two curves. For instance, Fig. 2 shows two parallel segments C and S with length 1 in R2, separated by a small distance ε. We want the norm of the difference µC − µS of the corresponding vector-valued measures to be small since the segments are close to each other.

Fig. 2
Two parallel segments C and S, separated by a small distance ε

This can be achieved if one puts constraints on the variations of the test vector fields w so that the evaluation of the curve representors µC and µS on w gives close results. To do so, we define a norm on vector fields in the form of


where L is a linear differential operator of sufficient order such that the following property on the norm is guaranteed

w1,=supx[set membership]Rd{|w(x)|+i=1d|[partial differential]w[partial differential]xi(x)|}cWwW

for some fixed constant cW > 0. In this setting, test vector fields w will belong to the Hilbert space W associated with the norm ‖·‖W.

Proposition 1

If property in (6) is satisfied, then the linear functional µC defined in (4) is a continuous linear functional on W.

Proof For any w [set membership] W,

|left angle bracketµC|wright angle bracket|01|γC(s)·w(γC(s))|ds01|γC(s)||w(γC(s))|dsw1,01|γC(s)|ds=w1,lCcWwWlC,

where lC is the length of C.

We can therefor compare vector-valued measures via the dual norm, defined by:

µCW*=˙sup wW1left angle bracketµC|wright angle bracket.

Going back to example of Fig. 2, the evaluation of the difference µC − µS is

left angle bracketµSµC|wright angle bracket=01γS(s)·w(γS(s))γC(s)·w(γC(s)) ds.

With γS(s) = (s, 0), γC(s) = (s, ε), γC(s)=γS(s)=(1,0), we get

left angle bracketµSµC|wright angle bracket=01(w1(s,0)w1(s,ε)) ds,

where w1(s, 0) denotes the first coordinate of vector w(s, 0). The integrand w1(s, 0) − w1(s, ε) is bounded by ε ‖w1,∞, which in turn is bounded by εcWwW from property (6). Hence the dual norm ‖µS − µCW* is bounded by εcW. When ε decreases, the norm ‖µC − µSW* goes to zero, as expected.

2.3 Explicit Expression of the Dual Norm with Kernels

The explicit expression of the dual norm involves the Green’s function kW(x, y) of the operator L*L, where L* is the adjoint of L. kW(x, y) is also the reproducing kernel of the space W. It is defined by the property

RdLkW(x,y)ξ·Lw(y) dy=left angle bracketkW(x,·)ξ,wright angle bracketW=w(x)·ξ,

where left angle bracket·, · right angle bracketW represents the inner product of two vector fields in Hilbert space W. Note that since we are dealing with vector fields, kW(x, y) is a matrix operator, we need to apply it to a vector ξ [set membership] Rd to get a proper definition.

To fix ideas at this point, one may choose L such that L*L = (Id − α2[nabla]2)k for some constants α > 0 (scale parameter) and integer k, which gives the Sobolev Hk Hilbert space. The corresponding kernel operator has explicit expression in terms of modified Bessel functions of the second kind (McLachlan and Marsland 2007).

Lemma 1

The norm of the vector-valued measure µC in the dual space W* is equal to

µCW*2=0101[kW(γC(s),γC(r))γC(s)]·γC(r) ds dr.


left angle bracketµC|wright angle bracket=01γ(s)·w(γ(s))ds=01left angle bracketkW(r(s),·)γ(s),wright angle bracketWds=left angle bracket01kW(r(s),·)γ(s)ds,wright angle bracketW.

By Cauchy-Schwartz inequality, this achieves the maximum for ‖wW = 1 if


Therefore, we have

µCW*2=01kW(γ(s),·)γ(s)dsW2=0101left angle bracketkW(γC(s),·)γC(s),kW(γC(r),·)γC(r)right angle bracketW ds dr=(a)0101[kW(γC(s),γC(r))γC(s)] ·γC(r) ds dr,

where (a) follows from the reproducing property of the kernel kW stated in (7).

The expression of the dual norm indicates that the central object in computations is the kernel kW. Mathematically it is possible to build the Hilbert space W and its dual starting from a given positive kernel (without any knowledge of operator L), with enough regularity to satisfy the property in (6). Links between properties of the kernel and properties of the space in this setting can be found in Glaunès (2005). In short, kW should be twice differentiable and vanish at infinity up to order 1. If one additionally requires invariance with respect to rigid motions, one can prove that kW is defined by two scalar functions h, h[perpendicular] : R+R such that for every x, y [set membership] Rd, kW(x, y)α = h(|xy|2)α if α is colinear to xy, and kW(x, y)α = h[perpendicular](|xy|2)α if α · (xy) = 0. For our applications we only consider “scalar” kernels kW(x, y) = h(|xy|2)Id i.e. such that h = h[perpendicular], and in practice we tried Gaussian h(r2)=er2/σW2 and Cauchy h(r2)=11+r2/σW2 kernels. However it would be interesting to study the effect of using non scalar kernels.

3 Curve Matching via Deformation Maps

3.1 A General Variational Formulation

In the previous section, we introduced representations of curves by vector-valued measures and defined a computable norm in the space of measures to compare them. These tools give us a practical way to define a matching functional between two curves C and S in Rd. Let [var phi]: RdRd be a deformation map which is supposed to give an approximate correspondence between C and S. The accuracy of this matching can be measured by:

E([var phi](C),S)=˙µ[var phi](C)µSW*2

which can be rewritten in terms of kernel kW, according to (8),

E([var phi](C),S)=0101[kW(γ[var phi](C)(s),γ[var phi](C)(r))γ[var phi](C)(s)]·γ[var phi](C)(r) ds dr20101[kW(γ[var phi](C)(s),γS(r))γ[var phi](C)(s)]·γS(r) ds dr+0101[kW(γS(s),γS(r))γS(s)]·γS(r) ds dr.

The first and last terms are intrinsic energies of the two curves. Roughly, the more curved they are, the higher energy they possess, since flat shaped parts of the curves tend to vanish in the space of vector-valued measures. The middle term penalizes to mismatching between tangent vectors of S and those of [var phi] (C).

Optimal matching between curves will be defined via minimization of a functional composed of the matching term E and optionally a regularization term to guarantee smoothness of the transformations. We model deformation maps [var phi] in the LDDMM framework. Certainly any other model of transformations could be directly applied at this point. For instance, the similarity transformation (rotation, translation, and uniform scale) is derived in Appendix B.

3.2 Large Deformation Diffeomorphic Metric Curve Mapping

In the LDDMM setting, we define the optimal matching, [var phiv with circumflex] between two curves C and S as [var phi]^=[var phi]1υ^t, [upsilon macr]t being the minimizer of the energy

JC,S((υt)t[set membership][0,1])=˙γ01υtV2 dt+µ[var phi]1υ(C)µSW*2

and [var phi]1υ is defined in (1).

The first term is the regularization term which controls the smoothness of the diffeomorphism, and the second term is the data attachment to quantify mismatching between deformed C and S. A general result (Trouvé 1995; Dupuis et al. 1998) guarantees existence of the solution to this minimization problem.

In Glaunès et al. (2006), such variational problems for 2D curve matching in the large deformation framework were studied, and precise forms of the optimal solutions were given. In this paper we will focus on the discrete formulation, leading to a practical matching algorithm.

4 The Discrete Model

4.1 Representation with Vector-Valued Dirac Functionals

First suppose that the curve C is continuous and discretized in a sequence of points x = (x1,…, xn) such that each xi belongs to the curve: for any parametrization γC we have xi = xsi for si [set membership] [0, 1]. We do not make any assumption on the way points xi are distributed along the curve. If C is a closed curve then we assume xn = x1. We can associate to this sequence of points a specific measure given by vector-valued Diracs:


where cx,i=12(xi+xi+1) is the center between points of xi and xi+1, and τx,i = xi+1xi gives an approximate tangent vector. This vector-valued measure is also defined by its pairing with smooth vector fields w:

left angle bracketµx|wright angle bracket=i=1n1τx,i·w(cx,i).

The important point here is that both the curves and their discrete representants (possibly with different sampling rates) are embedded into the same space of vector-valued measures. Hence it will be possible to compare them via the dual norm introduced earlier.

If C is piecewise continuous, we build a discrete representation by summing the discrete measures µx corresponding to each connected component.

4.2 Expression of the Dual Norm

The dual norm of the discrete representant µx expresses conveniently in terms of the kernel kW:


Note that this is not an approximation of the dual norm, but an exact expression of the dual norm of the discrete approximation of the curve. The following theorem shows that this discrete representant is a good approximation of the curve in space W.

Theorem 1

If the length of the piece of curve connecting two consecutive points dC(xi,xi+1)=sisi+1|γC(s)|ds is bounded by δ, then ‖µx − µCW* ≤ cWδlC.

Proof For any w [set membership] W,

left angle bracketµxµC|wright angle bracket=left angle bracketµx|wright angle bracketleft angle bracketµC|wright angle bracket=i=1nτi·w(ci)01γC(s)·w(γC(s))ds=i=1n(τi·w(ci)    sisi+1γC(s)·w(γC(s))ds).

We have τi=γC(si+1)γC(si)=sisi+1γC(s)ds, so we can write:

left angle bracketµxµC|wright angle bracket=i=1nsisi+1γC(s)·(w(ci)w(γC(s))) dsi=1nsisi+1|γC(s)||w(ci)w(γC(s))|ds.

Now |w(ci)w(γC(s))||w1,|ciγC(s)|cW·wWδ, so we get

left angle bracketµxµC|wright angle bracketcWwWδi=1nsisi+1|γC(s)|dscWwWδlC.

4.3 Matching Functional, Discrete Case

When curves C and S are discretized by point sequences x=(xi)i=1n and y=(yj)j=1m, the sequence of transformed points z=([var phi](xi))i=1n gives a discretization of the deformed curve [var phi](C). Here we will derive algorithms based on the following matching functional:


which is explicitly


Note that curves C and S are not necessarily discretized into the same number of points from the above equation. This expression was also given in Qiu et al. (2007). We show the derivation of gradient of E(z, y) with respect to [var phi] in Appendix A.

4.4 LDDMM Variational Formulation in the Discrete Setting

Let x=(xi)i=1n and y=(yi)i=1m be discrete samples of curves C and S respectively, and define the trajectories xi(t) [equals, single dot above] [var phi]t (xi) for i = 1, …, n. The optimal matching between samples x and y will be defined as a minimizer of

J=γ01υtV2 dt+i=1nτz,iδcz,ij=1mτy,jδcy,jW*2,

where z=([var phi]1(xi))i=1n=(xi(1))i=1n. The second term of J depends only on the positions of the finite number of points ([var phi]1(xi))i=1n. In such a situation, like all point-based matching problems in the large deformation setting (Joshi and Miller 2000; Camion and Younes 2001), the following reasoning apply: starting from any arbitrary state [upsilon macr]t of the time-dependent vector fields, let us find the optimal vector fields υt that keep trajectories xi(t)=[var phi]tυ^t(xi) unchanged and minimize the functional. Since trajectories are fixed, in particular the end-points (xi (1)) are also fixed; hence the second term of J is constant for this sub-problem and we are left to find vector fields which minimize the V-norm and satisfy for each time t and each 1 ≤ in the constraint υt (xi (t))= xi (t) where xi (t) are the fixed trajectories. This is a classical interpolation problem in a Hilbert space setting, for which the solution is expressed as a linear combination of spline vector fields involving the kernel operator kV:


where αk(t) are referred to as momentum vectors. These vectors can be computed from trajectories xi(t) by solving the system of linear equations

dxi(t)dt=j=1nkV(xj(t),xi(t))αj(t),   i=1,,n.

Now coming back to the main variational problem, we can restrict the domain of minimization to the subspace of vector fields of the form given by expression in (14), and choose either trajectories xi(t) or momentum αi(t) as minimization variables. The choice of αi(t) avoids the inversion of linear systems, since in this case (15) is seen as an ODE system in the unknown xi(t).

The regularization term can be rewritten as a function of variables αi (t) and xi (t). From (14) and the reproducing property of kernel kV, we get

01υtV2 dt=01i=1nj=1nleft angle bracketkV(xi(t),·)αi(t),kV(xj(t),·)αj(t)right angle bracketV dt


Hence the explicit expression of functional in (13) is the following:


The derivation of the gradient of J with respect to αi(t) is given in Appendix D.

4.5 Algorithm

We used a conjugate gradient routine to perform the minimization of functional J in (18) with respect to variables αi (t). Of course any other optimization scheme could be considered at this point. The different steps required to compute the functional and its gradient for each iteration are the following:

  1. from momentum vectors αi (t), compute trajectories xi (t) by integrating the system of ordinary differential equations (ODE) using (15)
  2. evaluate J from (18)
  3. compute vectors ηi (t) by integrating the system of ODE in (28) in Appendix D
  4. compute gradient ([nabla]J)i (t) = 2γαi (t)+ ηi (t).

All time-dependent variables were evaluated on a uniform grid t1 = 0, …, tT = 1 and a predictor/corrector centered Euler scheme was used to solve the systems of ODE in (15) and (28). The complexity of each iteration is of order dTN2, where N = max(n, m). To speed up computations when N is large, all convolutions by kernels kV and kW are accelerated with a multipole method (Yang et al. 2003), which reduces the complexity to dTN log(N).

5 Results

In our experiments, we first used the similarity transformation to transform curves into the same orientation and scale, as described in Appendix B. Then, we applied the curve matching under the LDDMM setting to deform one curve to the other. Kernels of kV (x, y) and kW(x, y) were radial with the form e|yx|2σ2Id, where Id is the d × d identity matrix. σV and σW represent the kernel sizes of kV and kW, respectively. They were experimentally adjusted.

5.1 2D Curve Matching


Shown in Fig. 3 are two-dimensional examples where the shape of curves represent bones (the first row), or birds (the second row), or hands (the third row). In this experiment, we chose σV = 0.25 and σW = 0.1 and deformed source shapes to target shapes. The first column shows the target shape, while the rest columns show source shapes. In each panel of the second and third rows, the blue curve represents the source shape before the curve matching and the green curve represents the deformed source shape after the curve matching. The background grid shows the deformation field. Figure 4 shows the time sequence of deformation of the source configuration. The first panel shows the source hand in green and the target hand in red. The last panel illustrates how close the deformed source curve (green) are to the target curve (red) after the curve matching. The middle panels demonstrate the smooth deformation as the source hand moves to the target hand along the time.

Fig. 3
Examples for plane curve matching. Bone, bird, and hand examples are respectively shown in rows. The first column shows target shapes. The second and third columns show source shapes. Blue curves are source shapes, while green curves are deformed source ...
Fig. 4
Panels from the left to the right depict the sequence of geodesic mappings connecting the source hand to the target hand at time t = 0, 0.2, 0.4, 0.6, 0.8, 1. The source shape and target shape are respectively represented in green and red

Robustness Against Noise

One application of studying transformations is to classify objects into different shape groups. In our LDDMM framework, the velocity vectors transform one object to the other object and give the shortest geodesic path connecting these two objects. The length of the geodesic path defines a metric distance given in (2) and feature in classification. As the purpose of classification, the metric distance should be robust against noise to achieve accurate classification. In this section, we claim that the curve matching provides a more accurate metric distance measurement and overcomes noise effect on matching results compared with the landmark matching (Joshi and Miller 2000; Allassonnière et al. 2005) when outliers are chosen as landmark points.

Example A in Fig. 5 shows boats without and with a mast as the source and target curves in panels (a, b), respectively. The curve and landmark matchings are applied to transforming the source boat to the target boat and results are shown in panels (c, d). The matching term in (18) for the template is near to zero in the region of the mast so that there is no force to drive the source boat to match the mast. The peak point of the mast is chosen as landmark in the landmark matching procedure so the source boat is deformed to the target boat at the mast. However, the metric distance (1.63) in the landmark matching is significantly larger than the one (0.13) in the curve matching. If a threshold is selected to classify these boats, the landmark matching probably considers them in different shape categories while the curve matching recognizes them as boats. When classifying shapes this problem typically arises as the result of imperfect extraction. The curve matching has less noise effect on the metric distance compared with the landmark matching. A similar example is given in Fig. 5(B), one curve without a gap (source shape) and the other with a small gap (target shape). The metric distance (1.37) in the landmark matching is larger than the one (0.72) in the curve matching.

Fig. 5
Comparison between the landmark and curve matchings. In each example (A) or (B), panels (a, b) show the source and target curves. Panels (c, d) show results from the curve and landmark matchings, respectively. The target and the deformed source shapes ...

Box Bump Experiment

To study the behaviour of our curve matching algorithm, we computed matchings between two shapes (cyan and red curves in Fig. 6) composed of a box with a semi-circle. Each row in Fig. 6 shows one example representing the time evolution of the optimal diffeomorphism by green curve in each panel. The desired behaviour, in which the bump slides along the edge, occurs only when the bumps locations on the two shapes are close (see the top row of Fig. 6). Otherwise the deformation that flattens one bump and reforms the other one is prefered, because the deformation cost is lower (see the bottom row of Fig. 6). Note again that our algorithm matches shapes up to a small error, measured by the metric W*. Avoiding such solutions would require to put higher-order (typically curvature) information into the metric W*, while our metric only relies on tangential information.

Fig. 6
Box bump experiments. Each row shows one mapping from a template (cyan curve) to a target (red curve) with evolution of [var phi]tυ^ at different times t [set membership] [0, 1] denoted by green curves. Top row: first experiment when two bumps on the ...

Computing all pairwise matchings between 17 shapes composed of a box with a semi-circle when the locations of the bump on each curve subsequently moves further away, we mapped the resulting distance matrix into a 2-dimensional coordinate system using classical multidimensional scaling (MDS) (Cox and Cox 2001). Figure 7 shows the locations of each shape in the first two scaling dimensions. This confirms the previous analysis: while shapes are correctly ordered into a sequence, they seem to follow a highly curved path, because the “flatten and reform” solutions bound the metric distances between shapes.


We demonstrate pattern classification of metric distances among curves shown in Fig. 8. The shapes in Fig. 8 were chosen from hand (10), dude (10), turtle (7), and fish (5). Special cases (hands and dudes with noise effect) are shown in the last two columns of rows (a, b), which are similar cases as in Fig. 5. The reason for choosing fish and turtles is that some of fish are outliers which are similar to turtles.

Fig. 8
Rows (ad) respectively show hand, dude, turtle, and fish

Our classification approach has three steps: map any two curves to each other and compute metric distances between them using (2), map metric distances into Euclidean space via MDS (Cox and Cox 2001), classify shapes via K-means clustering. Figure 9 shows the scatter plot of the curves in the features spaces generated from MDS. The numbers (1, 2, 3, 4) give the classification results and respectively denote the categories of hand, turtle, fish, and dude. Only two curves were wrongly classified as turtles. One is the hand shown in the second column from the right in row (a), the other is the fish shown in the last column of row (d). As we mentioned in the previous section, the small noise (e.g. the hand and dude in the last column of rows (a, b)) does not inference the classification using the curve matching approach. The successful performance rate of our classification approach is 93.75%.

Fig. 9
A scatter plot of feature dimensions from multidimensional scaling. 1: hand, 2: turtle, 3: fish, 4: dude

5.2 3D Curve Matching

In this section, we illustrate matching results of 3D curves. All curves were extracted from brain cortical surfaces and matched to the template curve. For the purpose of visualization, we show deformed brain surfaces where the deformation fields obtained from the curve matching were interpolated.


One major application of the curve matching in medical imaging is to transform one anatomical structure to the other using anatomical curve information. For instance, curves in the brain anatomy can be sulci or gyri that are conventional boundaries of cortical regions. In this experiment, we demonstrate how to warp one cortical surface to the other using the curve matching method. Cortical surfaces shown in the following were generated from MRI image volumes. A 3D subvolume encompassing the region of interest (ROI) was defined manually in each volume. A Bayesian segmentation using the expectation-maximization algorithm to fit the compartmental statistics was used to label voxels in the subvolume as gray matter (GM), white matter (WM), or cerebrospinal fluid (CSF) (Joshi et al. 1999; Miller et al. 2000). Surfaces were then generated at the GM/WM interface using a topology-correction method and a connectivity-consistent isosurface algorithm (Han et al. 2001, 2002).

Figure 10(a) shows the target of planum temporale (PT) surface located on the superior temporal gyrus (STG) posterior to the Heschl’s sulcus (HS) and extending to Sylvian fissure (SF). The structure of PT can be determined by three principal curves (STG, HS, SF) tracked by dynamic programming (Ratnanather et al. 2003). In this experiment, we fixed one surface as the target PT surface shown in Fig. 10(a). We extracted the three curves from twenty PT surfaces as curve objects and registered them with the curves on the target curves using our curve matching algorithm. Finally, the deformation field was interpolated to the interior region of PTs. The columns from left to right on rows (a)–(d) in Fig. 11 respectively illustrate the source PT surface, deformed source PT, source PT superimposed with the target PT, and deformed source PT with the target PT. As shown in Fig. 11, the performance for this dataset was quite good. Panel (a) in Fig. 13 shows cumulative distance distributions before and after the matching process. The cumulative distance distribution plots the percentage of vertices whose distance to the target surface is less than d as a function of distance d. Gray curves correspond to the source surfaces, while black curves are for the deformed source surfaces after the matching. Red curves are the mean surface distance graphs before and after the matching process from the bottom to the top, respectively. It is obvious that more vertices move close to the target surface after the matching.

Fig. 10
Panel (a) shows the curve defined as boundary curve for planum temporale cortical surface. Panel (b) depicts six curves describing the main shape of cingulate gyrus. Each curve is colored differently with indices in the same color scheme. The surfaces ...
Fig. 11
Three-dimensional deformation of planum temporale cortical surface. Panels (a)–(d) depict the deformation applied to planum temporale cortical surfaces. From the left to the right, panels respectively show the source and deformed source surfaces ...
Fig. 13
Panels (a) and (b) depict the cumulative distance distributions for planum temporale and cingulate surfaces. The cumulative distance distributions for the source surfaces are in gray, while the distributions for the deformed source surface are in black ...

We applied the same strategy to match ten cingulate gyri (CG) to the target CG shown in panel (d) of Fig. 12. CG is a part of the limbic cortex and has functions related to symptoms of Schizophrenia and Alzheimer’s disease. Panel (b) in Fig. 10 shows the six curves defined on the target CG. Each curve is color coded with the index in the same color scheme. Panels (a)–(c) in Fig. 12 give examples of the source CG surfaces and their deformed versions. To validate the matched results, we also computed the cumulative distance distributions shown in panel (b) of Fig. 13. This set of the CG curve mapping results have been used to detect the group differences in the CG thickness between healthy comparison controls and patients with schizophrenia, which illustrated a successful application of the curve mapping in neuroimaging studies (Qiu et al. 2007).

Fig. 12
Three-dimensional deformation of cingulate cortical surface. Panels (a)–(c) depict the deformation applied to cingulate cortical surfaces. In each panel, the top row respectively shows the source and deformed source surfaces from the left to the ...

Comparison with Landmark Matching

The first purpose of this section is to compare the matched results of PT surfaces via the curve matching algorithm with those from the landmark matching algorithm. Twenty PTs were used in this experiment. As for the PT structure, three corner points and three boundary curves of the PT are defined as point and curve objects across the population. In the landmark matching process, five points equally spaced on each boundary and the three corner points were chosen as point landmarks on each PT surface. Then, the landmark matching algorithm (Allassonnière et al. 2005; Joshi and Miller 2000) was applied to PTs to obtain the deformation field that was used to deform PT surfaces to the template. Since the PT structure is relatively flat and simple compared with other structures of the brain, the cumulative distance distribution between deformed source and target surfaces, as described in the previous section, can be used as quantitative validation measurement. Figure 14 shows the averages for the original source surfaces (solid line), deformed ones via the landmark matching (dashdot line), and deformed ones via the curve matching (dash line). This figure suggests that the cumulative distance distribution from the curve matching algorithm is above that from the landmark matching algorithm, which implies that the curve matching algorithm introduced here significantly improves matching results in terms of the surface distance measurement compared to those from the landmark matching method.

Fig. 14
Comparison of the cumulative distance distributions between the landmark and curve matchings. Solid, dashdot, and dash lines represent average cumulative distance distributions over the original source, deformed source surfaces via the landmark and curve ...

The second aim of this section is to compare deformation fields of CG surfaces from the landmark and curve matching algorithms. Landmark matching methods have been used to map curves and surfaces by extracting a pair labeled point sets from the curves as landmarks. However, manually labeling landmarks is labor-intensive as well as highly variable, and landmark matching methods are sensitive to the labeling. Consequently, it may result in giving incorrect mapping. In this experiment, we give such an example to show matching results and deformation fields via both landmark (Joshi and Miller 2000; Allassonnière et al. 2005) and curve matching methods under the LDDMM framework. In the landmark matching process, seventy two landmarks were selected from the six curves shown in panel (b) of Fig. 10, twelve landmarks per curve (Bakircioglu et al. 1998). After the landmark matching, the deformation field was interpolated for the other vertices on the CG. Figure 15 and Figure 16 give examples what the deformation fields are from the landmark and curve matchings. Panels (a)–(d) in both figures respectively show the target, source, deformed source CG surface after the landmark matching, and deformed source CG surface after the curve matching. Two red points are marked on the target surface in Figs. 15(a) and 16(a). Compared to the positions of these two points on the source surface, the relative positions are changed after the landmark matching process, which can be observed in panel (f) on both figures (zoom in version of the column in the triangulated mesh format). However, we do not observe such a change in the curve matching shown in panels (g) of Figs. 15 and and16.16. Thus, we conclude that the curve matching gives more proper deformation field due to incorporating the geometry information of the curves.

Fig. 15
Comparison of deformation fields between the landmark and curve matchings. Panels (a–d) show the target, source, deformed source surface after the landmark matching, deformed source surface after the curve matching. Red dots represents paired ...
Fig. 16
Comparison of deformation fields between the landmark and curve matchings. Panels (a–d) show the target, source, deformed source surface after the landmark matching, deformed source surface after the curve matching. Red dots represents paired ...

6 Conclusion

In this work we have proposed a mathematical framework in which both curves and their discrete samples are represented as vector-valued measures and embedded in a Hilbert space structure. Coupled with a model of deformations of the ambient space, we designed variational problems for curve matching takes convenient explicit formulations due to the use of the kernel of the Hilbert norm. We derived the case of the large deformations, given by integration of time-dependent vector fields. We also derived the case of the similitude transformation in Appendix B. We presented experiments for 2D and 3D curve mappings, and validated the accuracy of matching as well as demonstrated the behaviors of the curve mapping. Results indicate that the curve matching gives good matching and proper deformation fields compared to the landmark matching algorithm (Joshi and Miller 2000). We also proposed to apply this technique for shape classification using the metric distances provided by the large deformation framework.

Fig. 7
2-dimensional scaling plot of the distance matrix computed from pairwise matchings of 17 box bumps shapes


The work reported here was supported by grants: NIH P41 RR15241, NSF DMS 0456253, and National University of Singapore start-up grant R-397-000-058-133.

Appendix A: Gradient of E

As the first step, we compute the derivative of the dissimilarity measure between two sequences of points (zi)i=1n and (yj)j=1m with respect to z, given in (19):


Let 1 ≤ qn and zq + εzq be a small perturbation of point zq [set membership] Rd, with other points fixed. Since cz,i=12(zi+zi+1) and τz,i = zi+1zi, we have

[partial differential]cz,i[partial differential]ε={12z˜qif i=q or i=q1,0otherwise,[partial differential]τz,i[partial differential]ε={z˜qif i=q1,z˜qif i=q,0otherwise.

We consider kernels of the form kW(x, y) = h(|xy|2) and we note hx,y,i,j = h(|cx,icy,j|2) [set membership] R and hx,y,i,j=[partial differential]h(|cx,icy,j|2)[partial differential]cx,i=2h(|cx,icy,j|2)(cx,icy,j)[set membership]Rd for any sequences of points x, y and indices i, j. We get

[partial differential]E[partial differential]ε=j=1n1(hz,z,q,j·z˜q)(τz,q·τz,j)+j=1n1(hz,z,q1,j·z˜q)(τz,q1·τz,j)+2j=1n1(hz,z,q1,jhz,z,q,j)(z˜q·τz,j)2j=1n1(hz,y,q1,jhz,y,q,j)(z˜q·τy,j)j=1m1(hz,y,q,j·z˜q)(τz,q·τy,j)j=1m1(hz,y,q1,j·z˜q)(τz,q1·τy,j).

Hence the gradient of E with respect to zq is given by


Appendix B: Similarity Transformation

B.1 Quaternion Representation of Similitude Transformation

Here we set d = 3 and look at transformations of the form [var phi](x) = Rx + b, where R is a similitude matrix defined by the direct product of uniform scale with orthogonal spatial rotation and b is a translation vector. The algebra of quaternions is a useful mathematical tool for formulating the composition of arbitrary spatial rotations, and establishing the correctness of algorithms founded upon such compositions. Here, we use the quaternion q = [a, s] to represent a similitude transformation, where a is a vector with three entries a1, a2, a3 to describe spatial rotation and s is a scalar. We do not restrict the length of q to 1 so that the scaling factor is also considered.

The group (Q, *) is defined as the set of quaternions Q with the law of composition *: Q × QQ, given by


The inverse of q is q1=1s2+a·aq*, where q* is the conjugate of q defined as q* = [−a, s]. According to the law of composition in (20), the similitude operation R on x [set membership] R3 is described as [Rx, 0] = q * [x, 0] * q*, where


B.2 Similarity Alignment of Curves

In the continuous setting, the alignment of two curves C and S is performed through minimization of the following functional:

J(s,a,b)=E([var phi](C),S)=µ[var phi](C)µSW*2,

where [var phi](x) = Rx + b, with R parametrized by variables s and a via formula (21). When curves C and S are discretized via sequences x=(xi)i=1n and y=(yj)j=1m, we consider the functional


where zi = Rxi+ b. The derivative of this functional is computed in Appendix C.

Appendix C: Gradient of a Matching Functional with Respect to Parameters of Similarity Transformations

We need to compute the gradient of a dissimilarity functional


where z=(zi)i=1n=([var phi](xi))i=1n and [var phi](x) = Rx+ b, R computed from parameters s and a via formula (21): Rx = 2sa × x + 2(a · x)a + (s2 − |a|2)x. Let first compute the gradient of J with respect to parameter s. We have

[nabla]SJ=i=1n[partial differential][var phi](xi)[partial differential]s·[nabla][var phi](xi)E=i=1n(2sxi+2a×xi)·[nabla][var phi](xi)E.

The gradient of J with respect to parameters a can be computed as follows. Let aε = a + εã and take derivatives at ε = 0. We have

[partial differential]εJ=i=1n[partial differential]ε([var phi](xi))·[nabla][var phi](xi)E=i=1n2(a·a˜)(xi·[nabla][var phi](xi)E)+2s(a˜×xi)·[nabla][var phi](xi)E   +2(a˜·xi)(a·[nabla][var phi](xi)E)+2(a·xi)(a˜·[nabla][var phi](xi)E).

Since (ã × xi) · [nabla][var phi](xi)E = (xi × [nabla][var phi](xi)E) · ã, we obtain

[partial differential]εJ=i=1na˜·[2(xi·[nabla][var phi](xi)E)a+2s(xi×[nabla][var phi](xi)E)  +2(a·[nabla][var phi](xi)E)xi+2(a·xi)[nabla][var phi](xi)E)].

The first and third terms form a double wedge product: −(xi · [nabla][var phi](xi)E)a + (a · [nabla][var phi](xi)E)xi = (a × xi) × [nabla][var phi](xi)E so finally the gradient is

[nabla]aJ=i=1n2(a×xi)×[nabla][var phi](xi)E+2s(xi×[nabla][var phi](xi)E)+2(a·xi)[nabla][var phi](xi)E.

We can compute the derivative of the gradient of J with respect to translation b in the same way as for s and get

[nabla]bJ=i=1n[nabla][var phi](xi)E.

The gradients with respect to s, a, and b in (22, 23, 24) all involve the term [nabla][var phi](xi)E given in (19) (with zi = [var phi](xi)).

Appendix D: Gradient of a Point-based Matching Functional in the Large Deformation Setting

We rewrite J into a matrix form


where αt and xt are vectors of momentum and coordinates of vertices at time t. E denotes the matching functional.

Lemma 2

The gradient of J with respect to momentum αt with metric given by matrix kV (xt, xt) at time t is of the form ([nabla]J)t = 2γ αt + ηt, where vector ηt=[nabla]x1E+t1([partial differential]xs[kV(xs,xs)αs])*(ηs+γαs)ds and [nabla]x1E is given in (19).

Proof Our first goal is to compute variation of velocity fields υt and trajectories xt given a variation αϵ = α + ϵ[alpha].

We have




For a variation αϵ = α + ϵ[alpha], computing the derivative with respect to ϵ in (26) yields the corresponding variation of velocity field [upsilon with tilde] = [partial differential]ϵυϵ as the form

υ˜t(x)=[partial differential]xt[kV(xt,x)αt]x˜t+kV(xt,x)α˜t,

where [partial differential]i is the derivative with respect to the ith variable in kV. Similarly, the corresponding variation of trajectory x˜t=[partial differential]ϵ[var phi]tυϵ(x) from (15) is of the form

x˜t=0t[partial differential]xs[kV(xs,xs)αs]x˜s+kV(xs,xs)α˜sds.

To obtain xt a manipulation is needed. The above equation can be rewritten as

dx˜tdt=[partial differential]xt[kV(xt,xt)αt]x˜t+kV(xt,xt)α˜t,

which is a linear differential equation in variable xt. By the variation of constants method, we are able to write the solution of xt as


where operator Rst gives the solution of the homogeneous equation

dRstdt=[partial differential]xt[kV(xt,xt)αt]Rst.

It remains to express the variation of the energy [partial differential]ϵJ in the form

[partial differential]ϵJ=2γ01α˜t·kV(xt,xt)αtdt+γ01αt·[partial differential]ϵ(kV(xtϵ,xtϵ)αt)dt+dx1E[partial differential]ϵ(x1ϵ).

Substituting [upsilon with tilde]t (x) and xt in [partial differential]ϵJ, we have the second term as

01αt·[partial differential]ϵ(kV(xtϵ,xtϵ)αt)dt  =01(t1Rts*[partial differential]xs[kV(xs,xs)αs]*ds)·kV(xt,xt)α˜tdt.

The third term can be rewritten as

dx1E[partial differential]ϵ(x1ϵ)=01dx1ERt1kV(xt,xt)α˜tdt.

We finally obtain

[partial differential]ϵJ=01kV(xt,xt)(2γαt+ηt)·α˜tdt.

The gradient [nabla]J [set membership] L2([0, 1], (Rd)n) at time t is given


Vector ηt=γ01Rts*[partial differential]xs[kV(xs,xs)αs]*ds+Rt1*[nabla]x1E and it can be written in the backward integral equation

ηt=[nabla]x1E+t1([partial differential]xs[kV(xs,xs)αs])*(ηs+γαs)ds.

The choice of the metric over momentum αt in the space of L2 is arbitrary. A reasonable choice of metric is given by kV (xt, xt) with two-fold advantages of being closer to the metric inducing the space of velocity field υ and simplifying the formula for the gradient. Therefore, the gradient ([nabla]J)t can be simplified as


Contributor Information

Joan Glaunès, MAP5, CNRS UMR 8145, Université Paris Descartes, 75006 Paris, France.

Anqi Qiu, Division of Bioengineering, National University of Singapore, 7 Engineering Drive 1, Block E3A 04-15, Singapore 117574, Singapore ; gs.ude.sun@aqeib.

Michael I. Miller, Center for Imaging Science, Johns Hopkins University, Baltimore, MD 21218, USA.

Laurent Younes, Department of Applied Mathematics and Statistics, Johns Hopkins University, Baltimore, MD 21218, USA.


  • Allassonnière S, Trouvé A, Younes L. Geodesic shooting and diffeomorphic matching via textured meshes. EMMCVPR. 2005:365–381.
  • Avants B, Gee JC. Geodesic estimation for large deformation anatomical shape and intensity averaging. NeuroImage. 2004;23:139–150. [PubMed]
  • Bakircioglu M, Grenander U, Khaneja N, Miller MI. Curve matching on brain surfaces using frenet distances. Human Brain Mapping. 1998;6(5–6):329–333. [PubMed]
  • Bakircioglu M, Joshi S, Miller M. Image processing: Vol. 3661. Proc. SPIE medical imaging 1999. Bellingham: SPIE; 1999. Landmark matching on brain surfaces via large deformation diffeomorphisms on the sphere; pp. 710–715.
  • Beg MF. Variational and computational methods for flows of diffeomorphisms in image matching and growth in computational anatomy. Ph.D. dissertation. Johns Hopkins University; 2003.
  • Beg MF, Miller MI, Trouvé A, Younes L. Computing large deformation metric mappings via geodesic flows of diffeomorphisms. International Journal of Computer Vision. 2005;61(2):139–157.
  • Besl P, McKay N. A method for registration of 3-d shapes. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1992;14(2):239–256.
  • Camion V, Younes L. Geodesic interpolating splines. In: Figueiredo M, Zerubia J, Jain K, editors. Lecture notes in computer sciences: Vol. 2134. EMMCVPR 2001. Berlin: Springer; 2001.
  • Cao Y, Miller M, Winslow R, Younes L. Large deformation diffeomorphic metric mapping of vector fields. IEEE Transactions on Medical Imaging. 2005a;24:1216–1230. [PubMed]
  • Cao Y, Miller MI, Winslow RL, Younes L. ICCV. Los Alamitos: IEEE Comput. Soc.; 2005b. Large deformation diffeomorphic metric mapping of fiber orientations; pp. 1379–1386.
  • Cox MF, Cox MAA. Multidimensional scaling. Boca Raton: Chapman and Hall; 2001.
  • Dupuis P, Grenander U, Miller MI. Variational problems on flows of diffeomorphisms for image matching. Quaterly of Applied Mathematics. 1998;56:587–600.
  • Durrleman S, Pennec X, Trouve A, Ayache N. Measuring brain variability via sulcal lines registration: a diffeomorphic approach. Int. conf. med. image comput. comput. assist. interv. 2007:675–682. [PubMed]
  • Feldmar J, Ayache N. Rigid, affine and locally affine registration of free-form surfaces. International Journal of Computer Vision. 1996;18(2):99–119.
  • Fillard P, Arsigny V, Pennec X, Hayashi K, Thompson P, Ayache N. Measuring brain variability by extrapolating sparse tensor fields measured on sulcal lines. Neuroimage. 2007;34:639–650. [PubMed]
  • Gee JC, Bajcsy RK. Elastic matching: Continuum mechanical and probabilistic analysis. In: Toga AW, editor. Brain warping. San Diego: Academic Press; 1999. pp. 183–196.
  • Glaunès J. Transport par difféomorphismes de points, de mesures et de courants pour la comparaison de formes etl l’anatomie numérique. Ph.D. dissertation. Université Paris 13; 2005.
  • Glaunès J, Trouvé A, Younes L. CVPR. Los Alamitos: IEEE Comput. Soc.; 2004. Diffeomorphic matching of distributions: A new approach for unlabelled point-sets and sub-manifolds matching; pp. 712–718.
  • Glaunès J, Trouvé A, Younes L. Modeling planar shape variation via hamiltonian flows of curves. In: Krim H, Yezzi A, editors. Statistics and analysis of shapes. Boston: Birkhauser; 2006.
  • Grenander U, Miller MI. Computational anatomy: An emerging discipline. Quarterly of Applied Mathematics. 1998;56(4):617–694.
  • Han X, Xu C, Prince JL. CVPR’2001 (Kauai, HI) Vol. 2. Los Alamitos: IEEE Comput. Soc.; 2001. A topology preserving deformable model using level set; pp. 765–770.
  • Han X, Xu C, Braga-Neto U, Prince J. Topology correction in brain cortex segmentation using a multiscale, graph-based algorithm. IEEE Transactions on Medical Imaging. 2002;21:109–121. [PubMed]
  • Helm PA, Younes L, Beg MF, Ennis DB, Leclercq C, Faris OP, McVeigh E, Kass D, Miller MI, Winslow RL. Evidence of structural remodeling in the dyssynchronous failing heart. Circulation Research. 2006;98:125–132. [PubMed]
  • Joshi SC, Miller MI. Landmark matching via large deformation diffeomorphisms. IEEE Transactions on Image Processing. 2000;9(8):1357–1370. [PubMed]
  • Joshi M, Cui J, Doolittle K, Joshi S, Van Essen D, Wang L, Miller MI. Brain segmentation and the generation of cortical surfaces. NeuroImage. 1999;9:461–476. [PubMed]
  • Joshi SC, Davis B, Jomier M, Gerig G. Unbiased diffeomorphic atlas construction for computational anatomy. NeuroImage. 2004;23:151–160. [PubMed]
  • Joshi AA, Shattuck DW, Thompson PM, Leahy RM. Registration of cortical surfaces using sulcal landmarks for group analysis of meg data. International congress series: Vol. 1300. New frontiers in biomagnetism; Proceedings of the 15th international conference on biomagnetism; 21–25 August 2006; Vancouver, BC, Canada. 2007. pp. 229–232. [PMC free article] [PubMed]
  • Klassen E, Srivastava A, Mio W, Joshi SH. Analysis of planar shapes using geodesic paths on shape spaces. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2003;26(3):372–383. [PubMed]
  • Leow A, Thompson PM, Protas H, Huang S-C. ISBI. Los Alamitos: IEEE Comput. Soc.; 2004. Brain warping with implicit representations; pp. 603–606.
  • McLachlan RI, Marsland S. N-particle dynamics of the Euler equations for planar diffeomorphisms. Dynamical Systems. 2007;22(3):269–290.
  • Michor PW, Mumford D. An overview of the Riemannian metrics on spaces of curves using the Hamiltonian approach. Applied Computational Harmonic Analysis. 2007;23(1):74–113.
  • Miller MI, Massie AB, Ratnanather JT, Botteron KN, Csernansky JG. Bayesian construction of geometrically based cortical thickness metrics. NeuroImage. 2000;12:676–687. [PubMed]
  • Miller MI, Trouvé A, Younes L. On the metrics and Euler-Lagrange equations of computational anatomy. Annual Review of Biomedical Engineering. 2002;4:375–405. [PubMed]
  • Mio W, Srivastava A. Elastic-string models for representation and analysis of planar shapes. CVPR. 2004;(2):10–15.
  • Qiu A, Younes L, Wang L, Ratnanather JT, Gillepsie SK, Kaplan G, Csernansky JG, Miller MI. Combining anatomical manifold information via diffeomorphic metric mappings for studying cortical thinning of the cingulate gyrus in schizophrenia. NeuroImage. 2007;37:821–833. [PubMed]
  • Ratnanather JT, Barta PE, Honeycutt NA, Lee N, Morris NG, Dziorny AC, Hurdal MK, Pearlson GD, Miller MI. Dynamic programming generation of boundaries of local coordinatized submanifolds in the neocortex: application to the planum temporale. NeuroImage. 2003;20(1):359–377. [PubMed]
  • Rettmann ME, Han X, Xu C, Prince JL. Automated sulcal segmentation using watersheds on the cortical surface. NeuroImage. 2002;15(2):329–344. [PubMed]
  • Schmidt FR, Clausen M, Cremers D. Shape matching by variational computation of geodesics on a manifold. In: Franke K, Müller K-R, Nickolay B, editors. Lecture notes in computer science: Vol. 4174. DAGM-symposium. Berlin: Springer; 2006. pp. 142–151.
  • Sharon E, Mumford D. 2d-shape analysis using conformal mapping. International Journal of Computer Vision. 2006;70(1):55–75.
  • Thompson P, Toga A. A surface-based technique for warping three-dimensional image of the brain. IEEE Transactions on Medical Imaging. 1996;15(4):402–417. [PubMed]
  • Thompson PM, Schwartz C, Lin RT, Khan AA, Toga AW. Three–dimensional statistical analysis of sulcal variability in the human brain. Journal of Neuroscience. 1996;16(13):4261–4274. [PubMed]
  • Thompson PM, Hayashi KM, Sowell ER, Gogtay N, Giedd JN, Rapoport JL, de Zubicaray GI, Janke AL, Rose SE, Semple J, Doddrell DM, Wang Y, van Erp TG, Cannon TD, Toga AW. Mapping cortical change in alzheimer’s disease, brain development, and schizophrenia. NeuroImage. 2004;23:S2–S18. [PubMed]
  • Trouvé A. An infinite dimensional group approach for physics based models. 1995. (Technical report). Electronically available at
  • Twining C, Marsland S, Taylor C. Measuring geodesic distances on the space of bounded diffeomorphisms; Proceedings of the British machine vision conference (BMVC); September 2002; Cardiff. 2002. pp. 847–856.
  • Vaillant M, Glaunès J. Inform. proc. in med. imaging: Vol. 3565. Lecture notes in comput. sci. Berlin: Springer; 2005. Surface matching via currents; pp. 381–392. [PubMed]
  • Welker W. Why does cerebral cortex fissure and fold? Cerebral Cortex. 1990;83:3–136.
  • Yang C, Duraiswami R, Gumerov N, Davis L. Improved fast gauss transform and efficient kernel density estimation. IEEE international conference on computer vision; 2003. pp. 464–471.
  • Younes L. Computable elastic distances between shapes. SIAM Journal on Applied Mathematics. 1998;58:565–586.
  • Zhang Z. Iterative point matching for registration of free-form curves and surfaces. International Journal of Computer Vision. 1994;13(2):119–152.