PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Math Imaging Vis. Author manuscript; available in PMC 2010 July 6.
Published in final edited form as:
J Math Imaging Vis. 2006 January 31; 24(2): 209–228.
doi:  10.1007/s10851-005-3624-0
PMCID: PMC2897162
NIHMSID: NIHMS189824

Geodesic Shooting for Computational Anatomy

Abstract

Studying large deformations with a Riemannian approach has been an efficient point of view to generate metrics between deformable objects, and to provide accurate, non ambiguous and smooth matchings between images. In this paper, we study the geodesics of such large deformation diffeomorphisms, and more precisely, introduce a fundamental property that they satisfy, namely the conservation of momentum. This property allows us to generate and store complex deformations with the help of one initial “momentum” which serves as the initial state of a differential equation in the group of diffeomorphisms. Moreover, it is shown that this momentum can be also used for describing a deformation of given visual structures, like points, contours or images, and that, it has the same dimension as the described object, as a consequence of the normal momentum constraint we introduce.

1. Introduction

Over the past several years we have been studying natural shapes using homogeneous orbits of imagery constructed via the action of transformation groups on exemplars or templates. The mathematical structure of group action as a model in image analysis has been pioneered by Grenander [13], the idea being to introduce the group actions in the very nature of the objects themselves, through the notion of deformable templates. Roughly speaking, a deformable template simply is an “object or exemplar” Itemp on which a group G acts and generates, through the orbit An external file that holds a picture, illustration, etc.
Object name is nihms189824ig1.jpg = G · Itemp, a whole family of new objects. The interest of this approach is to concentrate the modeling effort on the group G, and not on the family of objects An external file that holds a picture, illustration, etc.
Object name is nihms189824ig1.jpg.

Since the earliest introduction by Silicon Graphics Incorporated of special purpose graphics hardware for object rendering, group action as a model in image analysis has been the subject of a wide range of research in computer vision. Naturally, the analytical and computational properties of the low-dimensional matrix Lie groups form the core dogma of modern Computer graphics. In sharp contrast, however, for the study of imagery generated from natural or biological shapes, the finite dimensional matrix groups are replaced by their infinite dimensional analogue, the more general diffeomorphisms [7, 11, 16, 22, 23, 36].

The anatomical orbit or deformable template is made into a metric space with a metric distance between elements by constructing curves through the space of diffeomorphisms connecting them; the length of the curve becomes the basis for the construction, the metric distance corresponding to the geodesic shortest length curves. This gives rise to a natural variational problem describing the geodesic flows between elements in the orbit, with the solution of the associated Euler-Lagrange equations giving the optimal flow of diffeomorphisms and thus the metric between the shapes. The obtained setting shares several similarities with the mechanics of perfect fluids, for which the Euler-Lagrange equation has been derived by Arnold (Eq. (1) of [2]) for the group of divergence-free volume-preserving diffeomorphisms. As well these results become another example of the general Euler-Poincaré principle of [19] applied to an infinite dimensional setting.

Interestingly, as emphasized by Arnold [1] in his study, one of the most beautiful aspects of studying diffeomorphisms with a Lie group point of view is that many fundamental aspects which can be proved in the finite dimensional case can be formally extended to retrieve well-known equations of mechanics. One of the purposes of this paper is to develop infinite dimensional analogues, for the study of high dimensional shapes via diffeomorphisms, of several of the well known properties of Lie groups in rigid body mechanics. In particular we shall focus on the interpretation of the Euler equation as an expression of the evolution of the generalized momentum of diffeomorphoc flow of least energy in both Eulerian and Lagrangian coordinates.

Such a point of view will link our geodesic formulation to a conservation of momentum law in Lagrangian coordinates providing a powerful method for studying and modeling diffeomorphic evolution of shape. It will imply that the momentum of the diffeomorphic flow at any place along the geodesic can be generated from the momentum at the origin, thus providing the vehicle for geodesic generation via shooting.

This same conservation of momentum of the diffeomorphic flow, allows us to derive equations for the geodesic evolution of the elements in the orbit It=Iφt1, t [set membership] [0, 1], I [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms189824ig1.jpg. This unifies various geodesic evolutions associated with orbits of sparse finite dimensional landmarked shapes as well as the evolutions of dense images. Of special interest is the fact that for the special case of image matching, geodesic evolution of elements in the orbit links us to the notion of normal motion familiar to the rapidly growing community working in level set methods. Interestingly, as we show, the momentum of the diffeomorphic flow is normal to the level sets associated with geodesic motion. By solving the partial differential equations which are associated with the conservation of momentum, we will be able to control by specifying the initial conditions (within a specific class of momentum which depends on the considered imaging problem) a wide range of arbitrarily large deformations; this provides new possibilities for learning shape models of deformable templates, or for designing new numerical matching procedures.

This second point of view in terms of the conservation of momentum law also sheds new light on a large number of high dimensional evolution based Active Model Methods in Computer Vision, including active snakes and contours [6, 12, 18, 20, 29, 31, 37, 40, 41, 43], active surfaces and deformable models [8, 9, 21, 24, 25, 28, 30, 33, 39, 40, 42]. In such approaches vector fields are defined which give the boundary manifold of the shape some velocity of motion, usually following the gradient of an energy to form an attractive force to pull the boundary. The power of such methods is that they parameterize motion only associated with a submanifold of the imagery, not the entire extrinsic background space. For example, to deform a planar simply connected shape via an active contour method, the dimension of the motion is determined by the dimension of the boundary of the region, which is substantially less than the dimension of the plane. Historically such approaches have not been studied globally as diffeomorphic action. In fact it is well known that such methods cannot prevent self intersection nor ensure topological consistency, for which the addition of other constraints become necessary [14, 15]. From the conservation law in Lagrangian coordinates describing geodesic motion in the metric space of diffeomorphic action, we introduce the normal momentum motion which constrains the momentum to the bounding manifold, and extends the velocity of motion of the shape to the entire background space, thereby giving the global property that the resulting integrated vector field generates a diffeomorphism on the entire extrinsic space. This in turn carries the smooth submanifold diffeomorphically. As the analysis shows, this global property seems to be required to generate geodesic motions.

2. The Basic Set Up

2.1. Right Invariant Metric on Group of Diffeomorphisms

The basic component of our models is the group of one-to-one, smooth, transformations (diffeomorphisms) of a bounded subset Ω [subset or is implied by] Rd. In this paper, we consider diffeomorphisms emerging as flows of non-autonomous differential equations. A time dependent vector field on Ω is a function:

v:[0,1]×ΩRd(t,x)v(t,x)

(v(t, x) will also be denoted vt(x)). The associated ordinary differential equation is

dydt=vt(y).

The flow of this ODE is a function [var phi]v which depends on time and space, such that

φvt(t,x)=vt(φtv(x))

and [var phi](0, x) = x for all x [set membership] Ω. We will also use the notation φtv(x) for [var phi](t, x) and

φsv1=φtvφsv1.
(1)

It is well-known that, under some smoothness conditions on v, such φtv is at all times a diffeomorphism of Ω.

The groups that we consider are precisely composed with such flows φ1v for v belonging to a specified functional class. More precisely, we assume that a Hilbert space g is given, the element of which being smooth enough vector fields on Ω, and denote the norm and inner product on this space by || ||g, left angle bracket,right angle bracketg. We now define (following [35]) the group G as the set of functions φ1v for time-dependent vector fields v satisfying

01||vt||gdt<,

i.e. for v belonging to L1([0, 1], g). We will always assume that g can be embedded in the space of ( C01(Ω,Rd), || ||1, ∞), containing vector fields on Ω, which vanish on [partial differential]Ω, where

||w||1,=||w||+||dw||.

From this definition, it appears that the main ingredient in the construction of G is the Hilbert space g.

Fixing v [set membership] g, one can define the linear form w [mapsto] left angle bracketv, wright angle bracketg, which will be denoted Lv. We therefore have the identity

(Lv,w)v,wg

(we use the standard notation (M, w) for the linear form M applied to w). By definition, Lv belongs to the dual, g* of g, and L can be seen as an operator L: gg* (this is the canonical duality operator of g on its dual). As we shall see, this operator turns out to be a key feature in our analysis. For the moment, we point out the fact that Lv is a linear form on g which is a space of smooth vector fields. Therefore, Lv itself can be a singular object (a generalized function, or a distribution). Here are a few examples of distributions M which qualify as elements of g*, under our running assumption that g is embedded in the space of C1 functions:

  1. L1 vector fields of Ω: if ψ: Ω [mapsto] Rd is integrable, define
    (M,v)=Ωψ(x),v(x)Rddx
  2. Let now μ be any measure on Ω, and ψ be μ integrable. Define
    (M,v)=Ωψ(x),v(x)Rddμ(x)
  3. Dirac measures: as a particular case of the previous, define, for x [set membership] Ω and a [set membership] Rd,
    (M,v)=a,v(x)Rd
    This will be denoted M=δx(a)g.
  4. Differential operators: if (fi, j, 1 ≤ i, jd) are integrable functions, define
    (M,v)=i,j=1dΩfijvixjdx

It is important to notice that, although L is defined in a rather abstract way in the previous lines, numerical procedures to compute geodesics can be derived most of the time from the knowledge of its inverse (of Green kernel) K = L−1. This K is a smoothing kernel, the choice of which, within a specific range of available kernels, being the starting point of any practical procedure. We do not detail numerical algorithms in this paper, but the reader can refer to [4, 5, 17, 23] for examples of choices of K.

2.2. Energy and Momenta

Consider a time-dependent diffeomorphism v [set membership] L1([0, 1], g), and let ( φ0tv, t [set membership] [0, 1]) be the associated flow, defined in the previous section. Along time, each point x [set membership] Ω, considered as a particle, evolves on the trajectory tφ0tv(x), its velocity at time t being by definition vt(φ0tv(x)). In other terms, for y [set membership] Ω, vt(y) is the instantaneous velocity of the particle which is at y at time t. It is called the Eulerian velocity at y at time t.

So, at each time, we have an Eulerian velocity field, y [mapsto] vt(y), and we define the kinetic energy of the system to be E(vt)=12||vt||g2. The total energy spent during the deformation path now is

Etotal(v)=1201||vt||g2dt.

Note that, in classical fluid mechanics, the kinetic energy is the sum of particle kinetic energies, which, for a homogeneous fluid with mass density given by ρ yields

||vt||2=ρ2Ωvt(y)2dy.

This is the L2 norm of v, which cannot be used in our context, since we require that g is embedded in C01 (we need some kind of Sobolev norms). However, in analogy with standard mechanical systems, we may define the global momentum of the system at time t to be the linear form Mt [set membership] g* such that E(vt) = (Mt, vt)/2, which, with the notation of the previous section, yields.

Mt=Lvt.

So, if vt is the Eulerian velocity field at time t, the momentum at time t is given by Lvt. It will be called the momentum in Eulerian coordinates.

2.3. Lagrangian and Eulerian Frames

The Eulerian frame, as introduced above, describes mechanical quantities as they are observed in the current configuration at each time. The Lagrangian frame, on the contrary, describes quantities as seen from the initial configuration. For example, the diffeomorphism φ0tv(x) provides the position at time t of the particle which was at x at time 0, which is a Lagrangian notion. For the velocity, we create a Lagrangian velocity field by pulling back the previously defined velocity vt, setting

vtl(x)=dds(φt1(φt+s(x)))s=0,i.e.vtl=(dφt)1(vtφt).

The operation

v(dφ)vφ1

defines a fundamental Lie group operation, and is called the adjoint action of G on its Lie algebra (which here is g), denoted Ad[var phi]v. We have obtained the relation

vt=Adφtvtl.

To interpret the adjoint action pictorially, the new vector field under the adjoint action v → (d[var phi])v * [var phi]−1 has to be interpreted as the transformation of v under the deformation generated by [var phi]. Figure 1 shows how the field vl at location x is transported by the flow to the value v(y) at location y = [var phi](x) by pushing forward (using [var phi]) the Lagrangian frame on which vl is drawn. Note that the orientation of the vector v(x) drawn on the deformed sheet is also changed (through the action of (d[var phi])).

Figure 1
Here is represented the deformation obtained by pulling back the Eulerian frame associated with ve and represents pictorially the adjoint action.

2.4. Momentum in Eulerian and Lagrangian Coordinates

The momentum Mt = Lvt, which has been defined in Eulerian coordinates, also admits a Lagrangian version. It can be computed by expressing the kinetic energy at time t, which is (Lvt, vt)/2, under the form (Mtl,vtl)/2,Mtl, being then the Lagrangian momentum. This is straightforward, since, by definition of an adjoint operator:

(Lvt,vt)=(Lvt,Adφtvtl)=(AdφtLvt,vtl).

This leads to the definition Mtl=AdφtLvt for the momentum in Lagrangian coordinates. The Lagrangian frame takes here the role of a Galilean, or reference frame, and we will retrieve below the fundamental principle of mechanics, which states that the Lagrangian momentum is constant over time along any energy minimizing path. Before this, we make a brief interuption in our discussion to describe the relation between the classical mechanics of a rigid body, and geodesic equation in matrix Lie groups. This simple description will help to understand the formalism in our infinite dimensional group of diffeomorphisms.

3. Euler Equation and Conservation of the Momentum for Lie Groups of Matrices

In this part, we derive the Euler equation for extremal paths of the kinetic energy in the case of Lie groups of matrices. This derivation is well-known in the context of classical solid mechanics [1], but this simpler case, which can be derived completly without too much technicalities, may be helpful to understand the case of diffeomorphisms.

Let G [subset or is implied by] x2133d(R) be a Lie group of d × d matrices with Lie algebra g. The case of interest is when G is a group of 3D rotations, which models the position of a rigid body with fixed center of mass. In this case, g is the vector space of antisymmetric 3 × 3 matrices. Let t [mapsto] gt be a trajectory in this group. Then, the angular velocity at [set membership] g is given by the equation dgtdt=gtat. This is to be related to our previous definition of Eulerian velocity, which was

dφtdt=vtφt

in which the (left) product of matrices is replaced by the (right) composition of functions.

Returning to the matrix case, we define the kinetic energy at time t to be (Jat, at)/2, for some symmetric positive definite operator J: gg*. In the case of the rigid body, the angular velocity can be identified to a 3-vector ωt, and J can be seen as a 3 × 3 matrix which only depends on the geometry of the object, called the inertia operator and, with some abuse of notation,

(Jat,at)=(Jωt,ωt).

(note that here the notation (,) refers to the sum of products of coordinates, i.e. the usual inner product on Euclidean spaces, after identification between g and g*). The total energy is

E(g)=1201(Jat,at)dt.

We retrieve again the analogy with the diffeomorphisms by letting

a,bg=(Ja,b)

so that J takes the role of L in the previous section.

We now compute the Euler equation for least energy paths between two fixed endpoints g0 and g1. We recall that the Lie bracket on g is [a, b] = abba.

Theorem 1

The Euler-Lagrange equation for the kinetic energy is given by

Jatada(Ja)=0.
(2)

where ada:gg is defined by duality through the equalities (adaf,b)=(f,adab)=(f,[a,b]).

Proof

Let (t [mapsto] g0(t)) be an extremal curve for the kinetic energy and ((t, h) [mapsto] g(t, h)) be a smooth deformation around h = 0 (g(t, 0) = g0(t)): Let a(t, h) and A(t, h) be such that

gt=gaandgh=gA.
(3)

Writing 2gth=2ght, we get gAa+gah=gaA+gAt i.e.

ah=At+[a,A]=At+adaA.
(4)

The curve A(t, h) can vary freely in g, with boundary conditions A(0, h) = A(1, h) = 0. From

ddh(a,agdt)h=0,

we get

a,At+adaAgdt=(Ja,At+adaA)dt=0
(5)

Using the duality relation, we get (Ja,adaA)=(ada(Ja),A) so that by integration by part, we finally obtain the Euler equation

Jatada(Ja)=0.
(6)

We know from Lagrangian mechanics that the motion of a body with inertial operator J without external forces are extremal paths of the kinetic energy. Hence, Eq. (6) is the evolution equation of a body. We recognize in this equation the momentum to the body MtbJat and the Euler equation is then:

Mbtada(Mb)=0.
(7)

The momentum in the body here is to relate to the momentum in Eulerian coordinates for diffeomorphisms. However, if we study the motion of the body in a fixed static reference frame, the momentum to the space denoted here Ms should remain constant in the absence of external forces. The momentum to the space is defined from Mbt by a change of reference frame:

MtbAdgt1(Mtb)
(8)

where Adg is the co-adjoint representation which is defined by duality through the equalities: (Adgf,b)=(f,Adgb)=(f,gbg1). We derive from the evolution equation for Mb, given by the Euler equation (7), the conservation of the momentum to the space Ms:

Theorem 2

Along extremal curves for the kinetic energy, Ms is constant:

dMsdt=0.
(9)

Proof

Indeed, we have

(dMtsdt,b)=ddt(Mts,b)=ddt(Jat,Adgt1b)

Since, ddt(Adgt1)=adatAdgt1, we get finally using Euler Eq. (6),

(dMtsdt,b)=(Jatada(Ja),Adgt1b)=0.
(10)

Thus, from the conservation of the momentum to the space, MtsM0s, we deduce that

Jat=Adgt(La0),
(11)

or equivalenty, for any b [set membership] g:

(Jat,b)=(Ja0,Adgtb)
(12)

These results are in fact true for any Lie-group with a left-invariant metric. As we now investigate, they can be formally extrapolated also for infinite dimensional groups of diffeomorphisms.

4. Geodesic Evolution of the Diffeomorphism and Conservation of Momentum

4.1. Euler Equation as Evolution Equation for the Momentum in Eulerian Coordinates

The derivation of the Euler equation for extremal paths of the kinetic energy in the case of finite-dimensional Lie groups can be carried out in full generality within the Lie theory framework, to lead to the law of conservation of momentum. A general computation can be found in [1]. In our infinite dimensional case, a rigorous derivation of this law is much harder, and must most of the time be obtained directly from variational and functional analysis arguments rather than with purely algebraic Lie group derivations. However, it is interesting, and quite informative, to use these derivations to obtain a formal proof of the conservation of momentum, without wondering too much about the well-posedness of the expressions. This will be done in the next paragraphs.

The first Euler equation provides the variations of the momentum in Eulerian coordinates. Before stating it, we need some definitions:

Definition 1

The adjoint action Ad of G on g and the associated adjoint action ad of g on itself are given with their dual operators Ad*, ad* by

Adφw=(dφ)wφ1,(Adφf,w)=(f,Adφw)
(13)

advw=[v,w]=(dv)w(dw)v,(advf,w)=(f,advw).
(14)

with [var phi] [set membership] G, w [set membership] g, f [set membership] g*.

Already at this point, one can point out the difficulty of the infinite dimensional problem: at the difference with the matrix case, if v, w belong to g, it cannot be guaranteed that it is still so for [w, v] = (dw)v − (dv)w: in situations of interest, g is, in fact, not a Lie algebra: Ad[var phi]w and advw do not necessarily belong to g. As a consequence, the definition of advf which has been given cannot hold without some restriction on f, in order to be able to extend it to vector fields which are brackets of elements of g. We however proceed with such formal computation without addressing these issues.

The geodesics are extremal curves for the kinetic energy. They satisfy an Euler equation giving the variation of the momentum in terms of the co-adjoint action operator on the momentum.

Statement 1

The Euler equation for the kinetic energy is given by

dLvdt+adv(Lv)=0.
(15)

When Lv [set membership] H (i.e. it is a function), one has

advLv=div(Lvv)+dvLv.
(16)

where div(u [multiply sign in circle] v) = duv + div(v)u.

These equations, which are derived below, are special cases of the Euler-Poincaré principle, described, for example in [19]. Equation (15) is formally identical to Eq. (2) in the matrix case, excepted for a sign difference arising from the switch from a left-invariance in the matrix case to a right-invariance in the diffeomorphism case.

Formal Justification

This is exactly as in the matrix case. Here again, let (t [mapsto] [var phi]t) be extremal and ((t, ε) [mapsto] [var phi]t,ε) be a smooth deformation around ε = 0, with the abuse of notation [var phi]t,0 = [var phi]t. Denote φt,εt=vt,εφt,ε,φt,εε=ηt,εφt,ε,vt,εε=ht,ε, still denoting vt,0 = vt, ηt,0 = ηt and ht,0 = ht. Our first step is to express ht in function of the other variables. For this, write

2φεt=2φtε

which yields

htφt+dφtvtηtφt=ηttφt+dφtηtvtφt

or (applying φt1 on the right to both terms) gives

ht=ηtt+dηtvtdvtηt=ηtt+[ηt,vt].

The1 first variation of the energy is given by

ddε01||vt,ε||g2dt=201vt,htgdt=01vt,dηtdt+[ηt,vt]gdt=01(Lvt,dηtdt)dt01(Lvt,advtηt)dt.

Since [var phi]t is extremal, this expression vanishes for all η (with η0 = η1 = 0), and a last integration by parts yields

dLvtdt+advtLvt=0,

which is Eq. (15).

We now prove Eq. (16) under the assumption that Lv is a function. By definition

(advLv,w)=(Lv,dvwdwv)=(dvLv,w)(Lv,dwv)

and the conclusion comes from Stokes’ theorem which states that (since v and w vanish on [partial differential]Ω)

(div(Lvv),w)=(Lv,dw.v).

It appears that the Euler equation (15) with advLv=divvLv+dvLv has been derived in [26] and subsequently [22] directly as the Euler-Lagrange equation for the kinetic energy by analytical means. This has been originally proved by Arnold in [3] for the motion of impressible fluid which corresponds to the case L = Id with the constraint div v = 0.

4.2. Conservation of Momentum in Lagrangian Coordinates

The Euler equation (15) is the evolution of the momentum in Eulerian coordinates. We recognize in this equation the momentum Mt [equals, single dot above] Lvt; the momentum in Eulerian coordinates evolves in time so as to balance the co-adjoint of the momentum thereby satisfying the associated Euler equation dMtdt+advt(Mt)=0 for extremal paths. However, the momentum in Lagrangian coordinates, identified in the introduction as Mtl=Adφt(Mt), remains constant in the absence of external forces, ddtMtl=0.

Statement 2

Along extremal curves for the kinetic energy, Mtl is constant:

dMtldt=0.
(17)

In particular, we have for all w [set membership] g,

(Lvt,w)=(Lv0,(dφt)1wφt).
(18)

Formal Derivation

Indeed, fix w [set membership] g and let f(ε)=(Mt+εl,w). We have, on the first hand f(0)=(dMtldt,w), and on the second hand (derivatives being evaluated at time ε = 0)

f(0)=ddε(Lvt+ε,Adφt+εw)=(dLvtdt,Adφtw)+(Lvt,ddεAdφt+εw)

Note here that Ad[var phi]*[var phi]= Ad[var phi] Ad[var phi]. Now, if [var phi]0 = id and dφεdε=v at ε = 0, we have for any w′,

ddε(Adφεw)ε=0=ddε((dφε)wφε1)ε=0=dv(w)dw(v)=advw.
(19)

Applying this to w′ = Ad[var phi]t−1 and v = vt, we get

f(0)=(dLvtdt,Adφtw)+(Lvt,advtAdφtw)=(dLvtdt+advt,Adφtw)=0

by Eq. (15). This completes the proof,

(dMtldt,w)=0.

Although the conservation of momentum has only been derived from formal arguments, we can check that, when it is satisfied, the generated deformation paths do provide extremal curves of the kinetic energy. The perturbation of the end point of the path ( φ0tv, t [set membership] [0, 1]) at time 1 under a perturbation vtε of vt is given by [22]:

ddεφ01vε(x)=01dφ0svφs1v(hsφ0sv)ds.
(20)

with hs(x)=dvsε(x)dε, the derivative being taken at ε = 0 (we have used the notation of Eq. (1)). Assume that this expression vanishes (so that the end point φ01vε remains unchanged). The first variation of the kinetic energy is given by 01vt,htLdt=01(Lvt,ht)dt=01(Lv0,dφ0tφt0vtφ0t)dt. Now, using (20) and the fact that d[var phi]0t [var phi]t1 = d[var phi]00 [var phi]01d[var phi]0t [var phi]t0, we get easily that 01dφ0tφt0vtφ0tdt=0 so that, by linearity,

01vt,htLdt=(Lv0,01dφ0tφt0htφ0tdt)=0.
(21)

4.3. Coadjoint Transport of Structures Along a Geodesic

For M [set membership] g*, the evolution tAdφt1M is called coadjoint transport. The fact that the momentum evolves by coadjoint transport along a geodesic implies the conservation of several properties whenever they are initially true, for Lv0. These properties will turn out to be of main importance for image processing applications.

In this section, we assume that M0 = Lv0 is given, and that the coadjoint transport Mt=Lvt=Adφt0(M0) is well defined at all considered times.

4.3.1. Coadjoint Evolution of the Support

Let Supp(M) denote the support of a momentum M [set membership] g*. It is defined by the complementary of the union of all open sets Ω′ [subset or is implied by] Ω which are such that (M, w) = 0 whenever w [set membership] g vanishes outside Ω′. We have the property:

Statement 3

If Mt=Adφt0(M0), then

Supp(Mt)=φ0tv(Supp(M0))

Indeed, assume that M0 vanishes Ω′ [subset or is implied by] Ω. Let w have its support included in [var phi]0t (Ω′). Then (Mt, w) = (M0, d[var phi]0t)−1 w * [var phi]0t and w * [var phi]0t vanishes outside Ω′, which implies that (Mt, w) = 0. Thus Supp(Mt) [subset or is implied by] [var phi]0t (Supp(M0)), and the reverse inclusion is true by inverting the roles of M0 and Mt.

As a first example, consider the case when M0 is finitely supported, and more precisely a sum of Dirac measures. This is legitimate since our hypotheses on L imply that Dirac measures belong to g*, therefore have the form Lv0 for some v0 [set membership] g. So, we assume that

(M0,w)=iai,w(xi)Rd,
(22)

where (xi)1≤in is a finite family of points in Ω (landmarks) and (ai)1≤in is a finite family of vectors in Rd. We write M0=i=1nδxiai, where, by definition

δxa:gRwa,w(x)
(23)

Denoting xi (t) [equals, single dot above] [var phi]t (xi), we obtain the fact that Mt is supported on {x1(t), …, xN (t)}. More precisely, a rapid computation shows that

Mt=i=1nδxi(t)ai(t)
(24)

with

ai(t)=(dxi(t)φt0)ai
(25)

so that the momentum remains a sum of Dirac measures. This is a special case of the property considered in the next section.

4.3.2. Coadjoint Transport of Measure

Measure-based momenta are given by

(M,w)=Ων0,wdμ0
(26)

where μ0 is a measure on Ω and ν0 is measurable and μ0-integrable. They generate a large class of geodesic evolutions, and have the attractive property that the momentum Lvt can be explicitly computed from the momentum at the origin.

Statement 4

Assume that (M0, w) = Ωleft angle bracketν0, wright angle bracket0 then the linear momentum functional evolves according to

(Mt,w)=Ωνt,wRddμtwhereνt(x)=(dxφt0)ν0φt0(x),μtμφt0,
(27)

i.e. μt (A) = μ([var phi]t0(A)) for any measurable set A.

The statement follows straightforwardly from the substitutions

(Mt,w)=Ων0(x),((dφ0t)1wφ0t)(x)Rddμ(x)=Ωdφ0t(x)φt0ν0(x),wφ0t(x)Rddμ0(x)=Ωνt(x),w(x)Rddμt(x)
(28)

Point-supported momentum evolution considered in the previous section, clearly is a particular case of this statement. As another illustration, consider the case of measures which are supported by submanifolds of Ω. In this case, the initial momentum is concentrated along the boundary Σ0 of a k-dimensional C 1 sub-manifold in Ω [subset or is implied by] Rd.

Let σ0 by the surface measure (given as the induced volume form on the sub-manifold) and let μ0 be supported by Σ0, such that for any smooth function on Ω, Ωf dμ0 = Σ0f α00 for some density α0 (not necessary positive) on the surface. Let ν0: Ω → Rd (the values of ν0 outside Σ0 will not be important) and define

(Lv0,w)=0ν0,wRdα0dσ0.
(29)

Using Statements 3 and 4, we get that the tranported measure μt is located on the transported sub-manifold Σt [equals, single dot above] [var phi]t0) (whose smoothness is preserved by the regularity of the diffeomorphisms in G) and can be written as μt = αt σt where σt is the k-dimensional volume measure on Σt. Moreover, if νt = d([var phi]t0)*γ0 * [var phi]t0, Statement 4 gives us the evolution of the momentum

(Lvt,w)=tνt,wαtdσt(y).
(30)

In the case where the sub-manifold is Ω itself, then σt = σ0 is the Lebesgue’s measure λ on Ω, and αt = α0 * [var phi]t,0|d[var phi]t,0|.

4.3.3. Coadjoint Transport of Orthogonality

The last property transported by geodesic evolution which is considered here is the normality with respect to a smooth submanifold of Ω. Since normality will be extensively studied in the next section, we here provide an illustration in a particular case.

Assume that ν0, in Eq. (26) can be expressed as

ν0=i=1rb0if0i
(31)

where (b0i)1ir and ( f0i) are two families of functions on Ω and 1 ≤ rd. Then, we get from Statement 4 that

νt=i=1rbtifti
(32)

where bti=b0iφt,0 and fti=f0iφt,0. Equation (32) can be interpreted as a normality property of the geodesic motion under initial condition (31). Indeed, let

0=i=1r(f0i)1({0}).

Assume that Σ0 is not empty and denotes N0(x)=Span{f0i(x)1ir} for any x [set membership] Σ0. Under appropriate transversality conditions, mainly

rdimN0,

Σ0 can be equipped with a structure of (dr)-dimensional C1 manifold and N0(x) is exactly the space of vectors normal to Σ at location x.

We then deduce easily that

t=i=1r(fti)1({0})

and equality (32) implies for any x [set membership] Σt

νt(x)Nt(x)
(33)

where Nt(x)=Span{fti(x)1ir} is the set of normal vectors to Σt at location x.

We deduce that if the momentum is normal to some k-dimensional sub-manifold Σ0, this normality property is preserved by coadjoint transport along a geodesic.

In the case of r = 1 and f01=f0, Σ0 is exactly the level set for threshold value 0 of f0 and the normality of the initial momentum to the level sets is preserved under geodesic motion. Since the threshold value is arbitrary, we deduce that the property is true for all level sets.

5. The Normal Momentum Motion Constraint

5.1. Heuristic Analysis

The conservation of momentum is a general property of geodesics in a group of diffeomorphisms with a right invariant metric. More can be said in the situation when diffeomorphisms are associated to deformations of geometric structures or images, which is the situation of interest for our applications. In this setting we are still looking for curves with shortest length in G, but we partially relax the fixed end-point condition by the constraint that the initial template is correctly mapped to the target: because there is a whole range of diffeomorphisms which deform one given structure into another, this condition is weaker than the fixed end-point condition, which means that there are more degrees of freedom for the optimization, and therefore more constraints on the minimum. For image matching, these additional constraints may essentially be summarized by the statement the momentum along the geodesic path is at all times normal to the level sets of the image. This is what we call the normal momentum constraint, which is described in this section.

We start with a heuristic analysis, for which I0, the image, is a smooth function defined on Ω. Let I1 be in the orbit of I0 for the group G of diffeomorphisms: there exists ψ [set membership] G such that I0 * ψ−1 = It. By compactness and semi-continuity arguments, one can prove the existence of a geodesic path [var phi] = ([var phi]t) such that

dG(Id,φ1)=infϕGdG(Id,ϕ)andI0=I1ϕ.
(34)

Let φ=(φtv) be such a solution and consider a first order expansion around t = 0, [var phi]t(x) [similar, equals] x + tv0(x) so that It=I0φt1(x)I0tI0,v0Rd. By definition, the cost to go from I0 to It is (still at first order) t|v0|L. However, any u [set membership] g such that left angle bracket[nabla]I0, uright angle bracketRd = left angle bracket[nabla]I0, v0right angle bracketRd, will lead to the same It so that the least deformation cost from I0 to It should be t|pI0 (v0)| where pI0 (v0) the unique solution of the minimization problem:

(P0):infuguLsubjectto:(v0u)(x),I0(x)Rd=0,xΩ
(35)

Since ( φtv) is a geodesic path minimizing the deformation cost from I0 to I1, it minimizes also the deformation cost from I0 to It yielding

v0=pI0(v0).
(36)

Introduce the set gI0 = {h [set membership] g | left angle bracket[nabla] I0(x), h(x)right angle bracketRd = 0, ∀x [set membership] Ω}: the constraints in An external file that holds a picture, illustration, etc.
Object name is nihms189824ig2.jpg can be restated as uv0 [set membership] gI0 so that pI0 (v) is the orthogonal projection of v on gI0, the space orthogonal to gI0. Hence, equality (36), translates to

hgsuchthatforallxΩ,I0(x),h(x)=0,wehavev0,hL=0.
(37)

Now, the fact that left angle bracket[nabla]I0(x), h(x)right angle bracket [equivalent] 0 means that h is a vector field which is tangent to the level sets of I0, and since left angle bracketv0, hright angle bracketL = (Lv0, h), we see that Lv0 vanishes when applied to any such vector field, or, that Lv0 is a linear form which is normal to the level sets of I0.

5.2. Rigorous Result

We now pass to a rigorous derivation of this property. Since it will be interesting to also consider images which are not smooth, we provide a new definition of the set gI0. Obviously, when I0 is smooth, h [set membership] gI0 is equivalent to the fact that, for any function f which is C1 on Ω, one has

ΩxI0,h(x)Rkf(x)dx=0.

Applying the divergence theorem (we assume that [partial differential]Ω is smooth enough and take advantage on the fact that elements of g vanish on [partial differential]Ω), we get

ΩI0(x)div(hf)(x)dx=0.

Since this has a meaning when I0 [set membership] L2(Ω), we now define

Definition 2

When I [set membership] L2(Ω), we denote

gI={hg:I,div(hf)L2=0forallfC1(Ω)}.

We still denote by pI the orthogonal projection on gI. The group G is assumed to be built as described in Section 2.1 (in particular g is continuously embedded in C1(Ω, Rd)).

Theorem 3. (Normal Momentum Constraint)

Assume that I0 [set membership] L2(Ω) and let φ=(φtv) be a geodesic path solution of (34). Then, for almost all t [set membership] [0, 1]

vtgIt.

The proof is given in Appendix A.

5.3. Examples

Consider again the case of smooth I, so that the condition h [set membership] gI is equivalent to h [set membership] g and for all x [set membership] Ω, left angle bracket[nabla]x I, h(x)right angle bracketRd = 0. Using notation (23), we get that I0(x),h(x)Rd=δx(I0(x)),hL so that gI0=(Span{δx(I0(x))xΩ}) and finally,

gI=Span{δx(I(x))xΩ}¯.
(38)

Vector fields v such that

(Lv,u)=Ωv(x),u(x)Rddμ(x),
(39)

where ν is normal2 to the level sets of I belong to gI and it can be shown that they form a dense subset.

For non-smooth I, we can similarly introduce ωf (I), as the unique element of g such that, for h [set membership] V, left angle bracketωf (I), hright angle bracket = left angle bracket I, div(h f)right angle bracketL2(Ω), and conclude that

gI=Span{ωf(I)fC1(Ω,R)}¯.
(40)

This implies that any element v [set membership] gI is such that Lv can be expressed as a limit

(Lv,u)=limNΩI(x)div(ufN)(x)dx

where f N [set membership] C1(Ω, R). The particular case when I is the indicator function of a smooth domain B [subset or is implied by] Ω (which can be interpreted as a smooth shape) is quite interesting. For x [set membership] [partial differential]B, let ν(x) be the outward normal to [partial differential]Ω and denote σB be the uniform measure on [partial differential]B. Then

ΩI(x)div(uf)(x)dx=Bdiv(uf)(x)dx=Bf(x)u(x),ν(x)RddσB(x).

In this case, we obtain a dense subspace of gI by considering elements v [set membership] g such that

(Lv,u)=Bv(x),u(x)Rddμ(x).
(41)

for some measure μ on [partial differential]B (the boundary of the shape).

Remark

We close this section with a technical, but important, remark. We have called normal momentum constraint the property that vt [perpendicular] gIt at almost all times. We have shown that this property is always true for geodesics minimizing (34). But there is another important issue, which is how much it constrains vt, or, in other terms, how big gI is for a given image I. That this is relevant, and sometimes non-trivial, may be seen from the following example: assume that we are in 2 dimensions (d = 2) and that I is a C1 image, with a non-vanishing gradient, at least on a dense subset of Ω. Then, on any point x such that [nabla] x I ≠ = 0, we can de-fine in a unique way a positively oriented orthonormal frame (τ (x) ν (x)) such that ν (x) = [nabla]x I/|[nabla]x I|. Then, if h [set membership] gI and h(x) ≠ = 0, we must have h/|h| = ± τ in a small ball around x. Now h, as an element of g must be smooth (depending on the choice made for L, and h/|h| has the same smoothness as h: this is impossible to achieve when τ itself is not smooth enough, which in such a case forces h(x) = 0. We thus get the property that h vanishes whenever τ (x) does not meet the smoothness requirements of g, which may very well be everywhere on Ω (or on a dense subset, which is the same since h is continous), in which case gI = {0} and the constraint is void, contrary to our intuition that the momentum should be aligned with ν. We see that, for the constraint to really be effective, we need some smoothness requirement on I. Fortunately, as illustrated by the example above, this smoothness is only required for the level sets of I, which must have a smooth boundary. With such an assumption, for example, one can show that if v [perpendicular] gI and

(Lv,u)=Ωξ(x),u(x)dμ(x)

for some measure μ on Ω and some vector field on ξ Ω, then ξ must be (μ-almost everywhere) orthogonal to the level sets of I. From a practical point of view smoothness of level sets may easily be obtained using algorithms such as mean curvature motion ([27]).

5.4. Conservation and Normality Property Check for Inexact Matching

Here, we give a brief account of situations in which proofs of conservation of momentum and the normality property can be carried on in a well-defined context, and retrieve the evolutions described in the previous section.

It is hard to make rigorous, in full generality, the variational argument we have used in the proof of Eq. (15). Notice that the well-definiteness of the conservation of momentum Eq. (17) is an issue by itself, since, when w [set membership] g, and [var phi]t is the diffeomorphism generated by a geodesic, there is a priori no reason for (d[var phi]t)−1w * [var phi]t) to belong to g: one must be able to define Lv0 on spaces which are bigger than g, which means that Lv0 needs to be somewhat smoother (as a distribution) than a generic element of g*.

However, there is a setting in which such a fact is true and easy to obtain: it is when the search for the geodesic is done with an approximation of the target, with some L2 penalty term added to control the error. We summarize this setting in the case of landmark matching, shape matching and image matching. In these three situations, we will retrieve the coadjoint transport of measure-based momenta. In all cases, results in [11, 35] ensure the existence of minimizers of the variational problem.

5.4.1. Inexact Landmark Matching

In this section, we assume that a measured space (An external file that holds a picture, illustration, etc.
Object name is nihms189824ig1.jpg, ρ), together with two measurable functions x, y: An external file that holds a picture, illustration, etc.
Object name is nihms189824ig1.jpg → Ω are given. The diffeomorphism [var phi] is searched to minimize

U(φ)=E(φ)+1σ2Jyiφ(xi)2dρ(i)

When ρ is discrete, this relates to point-based matching, x representing the landmark original positions and y giving the landmark target positions. If we express U in function of v, this requires the minimization of

U(v)=01||vt||L2dt+1σ2Jyiφ01v(xi)2dρ(i)

The main point here is to notice that the optimal solution v generates a geodesic in G between id and φ1v.

Proposition 1

Denote δxa the linear form on g such that (δxa,v)=a,v(x). Let v be a minimizer of Ũ. Then, letting xi(t)=φ0tv(xi)

Lvt=1σ2Jδxi(t)[(dxi(t)φt1v)(yixi(1))]dρ(i)
(42)

Proof

The proof of this result is a direct consequence of the identity, valid for s, t [set membership] [0, 1], v, h [set membership] L2([0, 1], g),

ddεφstv+εh(x)|ε=0=stdφsuvφutv(huφsuv)du.
(43)

the proof of which being carried on with usual ODE arguments and being omitted here. It is then straightforward to obtain (42), using the definition of the linear forms δxa, for x [set membership] Ω and a [set membership] Rd.

Equation (42) is a conservation of momentum equation for

Lv0=1σ2Jδxi[(dxiφ01v)(yixi(1))]dρ(i).

When An external file that holds a picture, illustration, etc.
Object name is nihms189824ig1.jpg is finite, this is equation (31) with ai=1σ2(dxiφ01v)(yixi(1)). Equation (42) now is exactly (24), since

ai(t)=(dxi(t)φt0v)ai=1σ2[dxi(0)φ01vdxi(t)φt,0v](yixi(1))=1σ2(dxi(t)φt1v)(yixi(1))

5.4.2. Inexact Shape Matching

We now consider the comparison of a binary set-indicator function, I0 = 1Ω0 ([Omega]0 [subset or is implied by] Ω having smooth boundaries) and a smooth function I1, through the minimization of

U(φ)=E(φ)+1σ2Ω1Ω0φ1(x)I1(x)2dx

over GL. We have

Proposition 2

Let v be a minimizer of U(φ01v) over L2([0, 1], g). Then

Lvt=1σ2Ω1(12I1)δφ1tv(x)[(dφ1tvφt1v)ν1]dσ1(x)
(44)

where Ω1=φ01v(Ω0), ν1 is the outward normal to [partial differential] Ω1 and σ1 is the volume measure on [partial differential] Ω1.

Proof

Taking a variation v + εh, the main issue is to compute the derivative of

12Ω1Ω0φ10v+εh(x)I1(x)2dx

This integral can be rewritten

Ω01v+εh(Ω0)(12I1(x))dx+12ΩI1(x)2dx

Since the last term is constant, we see that the problem boils down to the computation of the derivative of the first term, which can be written, after a change of variable and letting f1 = 1/2 − I1,

Ω0f1φ01v+εh(x)dxφ01v+εhdx

Define u by uφ01v(x)=ddεφ01v+εh(x). Simple computations, which can be, for example, found in [10], yields the fact that

ddεΩ0f1φ01v+εh(x)dxφ01v+εhdx=Ω1div(f1u)dx

Now the conclusion is a direct consequence of Eq. (43) and of the divergence theorem.

Here again, one straighforwardly checks that the conservation of momentum is satisfied. We have in particular

(Lv0,w)=1σ2Ω1(12I1)(dφ10vφ01v)ν1,wφ10v(x)dσ1(x)

Letting μ0=φ01v(σ1), and ν0(x)=(dxφ01v)ν1φ10v(x) this can be rewritten

(Lv0,w)=1σ2Ω0(12I1φ01v(x))v0,wdμ0(x)

which is under the general form of a measure-based momentum.

5.4.3. Inexact Image Matching

In this section, we let I0 and I1 be two smooth enough (say C1) functions defined on Ω (images). We consider the image matching problem which corresponds to minimizing, over G,

U(φ)=E(φ)+1σ2ΩI0φ1(x)I1(x)2dx

This problem is equivalent to minimizing

01||vt||L2dt+1σ2ΩI0(φ01v)1(x)I1(x)2dx

This matching problem has been studied, in particular in [4], to show that the optimal solution should satisfy, at each time t,

Lvt=1σ2dφt,1v(I0tvI1tv)I0tv
(45)

in which we have introduced the notation: I0tv=I0φt0v,I1tv=I1Φt1v, and |d[var phi]| for the Jacobian of [var phi]. This equation is in fact an equation of conservation of momentum, with

Lv0=1σ2dφ01v(I0I10v)I0

as can be deduced from Eq. (30), with α=1σ2dφ01v(I0I10v). Moreover, we can check also the normality property (31) which holds here with with r = 1 and f0 = I0. This allows us to conclude that for the geodesic path in the image space generated by inexact matching, the lifting of the path in G defines a geodesic for which the momentum stays normal to the level sets of the current image I0 * [var phi]t,0 at time t.

6. Geodesic Evolution in the Orbit

Thus far we have concentrated on the evolution of the flow of diffeomorphisms and its conservation of momentum. For all of our image understanding work we use the flow ([var phi]t, t [set membership] [0, 1]) to act on the elements in the orbits An external file that holds a picture, illustration, etc.
Object name is nihms189824ig1.jpg of a given template I = Itemp. Now we examine the geodesic flows in the orbit {It = I * [var phi]t, t [set membership] [0, 1]}, I [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms189824ig1.jpg, and provide the associated evolution equations.

6.1. Geodesic Evolution Equation of Landmark Points

Here we examine the finite dimensional landmark orbit denoted An external file that holds a picture, illustration, etc.
Object name is nihms189824ig1.jpgn, consisting of n-shapes IN = (x1, …, xn), each landmark (xi)1≤in is in Ω [subset or is implied by] Rd; correspondingly (ai)1≤in are a finite family of vectors in Rd. Denoting by xi (t) =̇ [var phi]t (xi), the trajectory in Ω of the point xi under the flow [var phi]t gives

Lvt=i=1nδxi(t)ai(t),
(46)

where ai (t) = (dxt (t) [var phi]t,0)*ai. From the identity

ddt((dxi(t)φt,0))=dxi(t)vt(dxi(t)φt,0),
(47)

we deduce that daidt(t)+(dxi(t)vt)ai(t)=0. To prove (47), one needs to remark that

dxi(t)φt,0=dφ0,t(xi)φt,0=(dxiφ0,t)1

Now,

ddt(dxiφ0,t)1=(dxiφ0,t)1ddt(dxiφ0,t)(dxiφ0,t)1=(dxiφ0,t)1(dφ0,t(xi)vtdxi,φ0,t)×(dxiφ0,t)1=dxi(t)φt,0dxi(t)vt

Hence we get the following geodesic evolution in the orbit of landmarks.

Proposition 3 (Landmark Transport)

The landmarks are transported along the geodesic according to the following equations with velocity vector field satisfying

vt=L1(i=1nδxi(t)ai(t))=i=1nK(xi(t),·)ai(t)

where K(x, y) is the Green kernel associated with L:3

{daidt(t)+(dxi(t)vt)ai(t)=0,dxidt(t)vt(xi(t))=0.
(49)

Note that the expression of vt from Eq. (48) can be introduced into the system (49), yielding an evolution equation which only depends on the landmarks in the orbit.

We notice the reduction of the vector field to the range space of the Green’s kernels travelling over the landmark trajectories is as in [17].

There is a straightforward extension of this result to geodesic curve evolution, in which x(0) is a parametrized curve σ [mapsto] x(σ) for σ [set membership] [0, L] and

Lv0=0Lδx(0,σ)ν0(σ)dσ

where ν0(σ) is normal to x(0). In this case, we have

Lvt=0Lδx(t,σ)νt(σ)dσ

with

  • νtt(t,σ)+(dx(t,σ)vt)ν(t,σ)=0,
  • xt(t,σ)vt(x(t,σ)=0.

6.2. Geodesic Image Evolution

Assume here that = α (x)dx has a density with respect to Lebesgue’s measure on Ω. In this case, t = |d[var phi]t0|α * [var phi]t0dx and

Lvt=dφt0αφt0It
(50)

From the conservation of momentum in Lagrangian coordinates for image based motion, we get for Lv0 = α0 [nabla] I0 that Lvt = αt |d[var phi]t,0|[nabla]It where αt = α0 * [var phi]t,0, It = I0 * [var phi]t,0. Let zt = αt |d[var phi]t,0| so that Lvt = zt [nabla] It.

Since ddt(ztφ0,t)=ddt(dφ0,tφt,0α0)=div(vt)φ0,tdφ0,tφt,0α0 we get dztdt+dzt(vt)+div(vt)zt. Moreover, we get easily Itt+It,νt=0. Hence we get the following geodesic evolution equation in image space.

Proposition 4 (Image Transport)

The image is transported along the geodesic according to the following equations: with vector field vt = L−1(zt [nabla]It):

  • dztdt+dzt(vt)+div(vt)zt=0,
  • Itt+It,vt=0.

Notice that these equations appear as a limit case of the evolution equations which have been studied in [34] for image comparison.

As illustrated above, the pair (I0, μ) provides a device for modeling deformations. In the cases we have studied, I0 was representing some geometrical structure (a curve, an image), which evolved with time accordi to the generated deformation, and μ, essentially quantifies the speed and direction of the deformation.

We get from this a natural way to represent the deformation of a template. Using Grenander’s original terminology, I0 would precisely be the template and μ is the generator of the deformation. Thus, fixing I0, and letting μ vary, we get a model which represents perturbations of the template.

An example of deformations of an image is provided in Fig. 2. The images have been obtained by solving Eq. (50) from an initial image g of a slice of macaque brain, and taking

Figure 2
Three random deformations of an image.

α(x)=X(x)1+xI0

where X is a Gaussian process.

7. Computational Results

The following results illustrate the computation of the momentum Lv0 = α0[nabla]I0 (as described in Section 6.2) from geodesic paths between two images. These geodesics are computed using F. Beg’s implementation of image matching, as described in [4]. In these experiments, the operator L is ([nabla] + c Identity),2 implemented via fast Fourier transform. The shooting algorithm solves the equation provided in Proposition 4 with initial condition (I0, z0 = α0).

Figure 3 shows the three objects studied, a smooth Gaussian bump for shift, circles for scale, and two mitochondria examin g both forward and inverse shooting.

Figure 3
Figure shows objects under translation (column 1), scale (column 2), and mitochondria 1 and 2.

Shown in Fig. 46 are examples illustrating the image based momentum and the diffeomorphisms generated via geodesic shooting. Figure 4 shows the results of the translation experiment. Panels 1 and 2 show I0 and I1; panel 3 shows the diffeomorphism generated via geodesic shooting applied to I0, and illustrates how solving the conservation of momentum equation allows to recover I1 from I0 and Lv0.

Figure 4
Rows 1, 2: Results from translation experiment. Panels 1, 2, and 3 show I0, I1, I0 * [var phi]10. Panels 4, 5 and 6 show α0, Lv0 versus α0[nabla]I0, and V0 and L−1α[nabla] I0.
Figure 6
Results from mitochondria 2. Panels 1, 2, and 3 show I1, I0, I1 * [var phi]10. Panels 4, 5 and 6 show α0, LV0 versus α0[nabla] I0, and V0.

Shown in panel 4 is the density α showing the concentration near the boundary of I0. Superposed in panel 5 are the predicted directions of the momentum, given by α0[nabla]I0, and the value Lv0 obtained from Beg’s algorithm. In almost all case they appear as one line indicating a good accuracy of the algorithm. Panel 6 indicates that the vector field V0 demonstrating that while α and the momentum Lv0 are highly localized, the velocity of motion extends over the entire object.

Shown in Fig. 5 are similar results for the scale experiment.

Figure 5
Rows 1, 2: Results from scale experiment. Panels 1, 2, and 3 show I0, I1, I0 *; [var phi]10. Panels 4, 5 and 6 show α0, Lv0 versus α0[nabla]I0, and V0 and L−1α[nabla] I0.

Shown in Fig. 6 are two sets of results for the geodesic shooting of the mitochondria. The organization of the results are the same as for the translation and scale experiments. Shown in Fig. 6 are two sets of results for the geodesic shooting of the mitochondria. The organization of the results are the same as for the translation and scale experiments.

Figure 7
Results from mitochondria 1. Panels 1, 2, and 3 show I0, I1, I0 * [var phi]10. Panels 4, 5 and 6 show α0, LV0 versus α0[nabla] I0, and V0.

Acknowledgments

Michael I. Miller was supported by grants from the National Institute of Wealth numbers P41-RR15241-01A1, U24-RR0211382-01, P20-MH071616-01.

Appendix A: Proof of Theorem 3

Proof

Since gIt is closed, we have to show that for almost all t, vt = pIt (vt). Denote htvtpIt(vt). For ε [set membership] [0, 1], let vtεvt+εht, and φtεφtvε (one can check, but we skip the argument, that tht is measurable and belongs to L2([0, 1], g), so that this variation is valid).

The proof essentially consists in showing that, for all 0 ≤ t ≤ 1

I0φt0=I0φt0ε.
(51)

Indeed, assume that this result is proved. Considering ε = 1 and t = 1, we deduce that I1=I0φ1,01. However, since left angle bracketht, vtright angle bracketL = left angle bracketvtpIt(vt), vtright angle bracketL = 0, we get vt+htL2=vtL2htL2. Since tvt corresponds to paths with lowest kinetic energy from I0

to I1, we deduce that 01htL2dt=0 and the proof is ended.

We now return to Eq. (51). Using the formula

dφt0(x)dt=dxφt0vt(x)

and letting qtε(x)=φt0φ0tε(x), we obtain

qtε(x)t=dφ0tε(x)φt0vtφ0tε(x)+dφ0tε(x)φt0(vtφ0tε(x)+εhtφ0tε(x))=εdφ0tε(x)φt0htφ0tε(x)

We first prove Eq. (51) under the assumption that I0 is C1. From the computation above, we have

t(I0qtε)(x)=εφt0(φ0tε(x))I0,dφ0tε(x)φt0htφ0tε(x)Rd=εφ0tε(x)It,ht(φt0ε(x)))Rd=0

since by definition of the projection pIt (vt), we have for any x [set membership] Ω

ht(x),It(x)Rd=vt(x),It(x)RdpIt(vt),It(x)Rd=0.

This implies I0φt,0φ0,tε=I0 which yields Eq. (51) in this case. When I0 is not smooth, the proof goes by showing that

tΩI0qtε(x)f(x)dx=ΩIt(x)div(htgtε)(x)dx

for smooth f on Ω and gtεφtεdφtε=f which can be done either by a direct (heavy) computation, or by using a density argument, based on the fact that, by the divergence theorem, this is true for smooth I0 (we skip the details).

Footnotes

1These arguments are purely formal since ht includes the Lie bracket which cannot be guaranteed to belong to g (in which case our variation would not be justified).

2When I has smooth level sets, we say that a vector field ν is normal to its level sets when, denoting by Ωi the set {Ii}, v(x) is normal to [partial differential] Ωi if x [set membership] [partial differential]Ωi for some i and x = 0 otherwise.

3The explicit form for L−1 depending on the kernel K is defined as follows. For any x, y [set membership] Ω, the bilinear form Kx,y on Rd × Rd defined by

Kxy(a,b)δxa,δybg=(δyb,L1(δxa))=L1(δxa)(y),bRd
(49)

Contributor Information

MICHAEL I. MILLER, Center of Imaging Science & Department of Biomedical Engineering, The John Hopkins University, 301 Clark Hall, Baltimore, MD 21218, USA.

ALAIN TROUVÉ, CMLA, ENS de Cachan, 61 Avenue du Président Wilson, 94235 Cachan CEDEX, France.

LAURENT YOUNES, Center for Imaging Science & Department of Applied Math and Statistics, The Johns Hopkin, 3245 Clark Hall, Baltimore, MD 21218, USA.

References

1. Arnold VI. Mathematical Methods of Classical Mechanics. 1. Springer; 1978. 1989.
2. Arnold V, Khesin B. Topological methods in hydrodynamics. Ann Rev Fluid Mech. 1992;24:145–166.
3. Arnold VI. Sur la géometrie differentielle des groupes de Lie de dimension infinie et ses applications à l’hydrodynamique des fluides parfaits. Ann Inst Fourier (Grenoble) 1966;1:319–361.
4. Beg MF, Miller MI, Trouvé A, Younes L. Computing large deformation metric mappings via geodesics flows of diffeomorphisms. Int J Comp Vis. 2005 February;61(2):139–157.
5. Camion V, Younes L. Geodesic interpolating splines. In: Figueiredo MAT, Zerubia J, Jain AK, editors. EMM-CVPR 2001, Vol. 2134 of Lecture Notes in Computer Sciences. Springer-Verlag; 2001. pp. 513–527.
6. Chen T, Metaxas D. Medical Image Computing and Computer-Assisted Intervention—Miccai 2000, Vol. 1935 of Lecture Notes in Computer Science. Springer-Verlag; 2000. Image segmentation based on the integration of Markov random fields and deformable models; pp. 256–265.
7. Christensen GE, Rabbitt RD, Miller MI. Deformable templates using large deformation kinematics. IEEE Trans Image Processing. 1996;5(10):1435–1447. [PubMed]
8. Cohen I, Cohen L, Ayache N. Using deformable surfaces to segment 3-d images and infer differential structures. Computer Vision Graphics Image Processing. 1992;56(2):242–263.
9. Cootes T, Taylor C, Cooper D, Graham J. Active shape models—Their training and application. Comp Vision and Image Understanding. 1995;61(1):38–59.
10. Delfour MC, Zolésio J-P. Shapes and Geometries. Analysis, Differential Calculus and Optimization. SIAM; 2001.
11. Dupuis P, Grenander U, Miller M. Variational problems on flows of diffeomorphisms for image matching. Quart App Math. 1998;56:587–600.
12. Garrido A, De la Blancà NP. Physically-based active shape models: Initialization and optimization. Pattern Recognition. 1998;31(8):1003–1017.
13. Grenander U. General Pattern Theory. Oxford Univeristy Press; 1994.
14. Grenander U, Chow Y, Keenan D. HANDS: A Pattern Theoretic Study of Biological Shapes. Springer-Verlag; New York: 1990.
15. Grenander U, Miller MI. Representations of knowledge in complex systems. J Roy Stat Soc B. 1994;56(3):549–603.
16. Grenander U, Miller MI. Computational anatomy: An emerging discipline. Quart App Math. 1998;56:617–694.
17. Joshi S, Miller MI. Landmark matching via large deformation diffeomorphisms. IEEE Trans Image Processing. 2000;9(8):1357–1370. [PubMed]
18. Kass M, Witkin A, Terzopolous D. Snakes: Active contour models. International Journal of Computer Vision. 1988;1(4):321–331.
19. Marsden JE, Ratiu TS. Introduction to Mechanics and Symmetry. Springer; 1994.
20. Mignotte M, Meunier J. A multiscale optimization approach for the dynamic contour-based boundary detection issue. Computerized Medical Imaging and Graphics. 2001;25(3):265–275. [PubMed]
21. Miller M, Joshi S, Maffitt DR, McNally JG, Grenander U. Mitochondria, membranes and amoebae: 1, 2 and 3 dimensional shape models. In: Mardia K, editor. Statistics and Imaging. II. Carfax Publishing; Abingdon, Oxon: 1994.
22. Miller M, Trouvé A, Younes L. On the metrics and Euler-Lagrange equations of computational anatomy. Annual Review of Biomedical Engineering. 2002;4:375–405. [PubMed]
23. Miller M, Younes L. Group actions, homeomorphisms, and matching: A general framework. International Journal of Computer Vision. 2001;41(1/2):61–84.
24. Montagnat J, Delingette H. Cvrmed-Mrcas’97, Vol. 1205 of Lecture Notes in Computer Science. Springer-Verlag; 1997. Volumetric medical images segmentation using shape constrained deformable models; pp. 13–22.
25. Montagnat J, Delingette H, Ayache N. A review of deformable surfaces: Topology, geometry and deformation. Image and Vision Computing. 2001;19(14):1023–1040.
26. Mumford D. Questions Mathématiques En Traitement Du Signal et de L’Image. Chap 3. Institut Henri Poincaré; 1998. Pattern theory and vision; pp. 7–13.
27. Osher S, Sethian JA. Front propagating with curvature dependent speeds: Algorithms based on Hamilton-Jacobi formulation. Journal of Comp Physics. 1988;79:12–49.
28. Pham D, Xu C, Prince J. Current methods in medical image segmentation. Ann Rev Biomed Engng. 2000;2:315–337. [PubMed]
29. Schultz N, Conradsen K. 2d vector-cycle defonnable templates. Signal Processing. 1998;71(2):141–153.
30. Sclaroff S, Liu LF. Deformable shape detection and description via model-based region grouping. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2001;23(5):475–489.
31. Staib L, Duncan J. Boundary finding with parametrically deformable models. IEEE Trans Pattern Analysis and Machine Intelligence. 1992;14:1061–1075.
32. Staib L, Duncan J. Model-based deformable surface finding for medical images. IEEE Trans Medical Imaging. 1996;15(5):1–13. [PubMed]
33. Terzopoulos D, Metaxas D. Dynamic models with local and global deformations: Deformable superquadrics. IEEE Trans Patt Anal Mach Intell. 1991;13:703–714.
34. Trouvé A, Younes L. Local geometry of deformable template. SIAM Journal of Mathmatical Analysis. to appear.
35. Trouvé A. Action de groupe de dimension infinie et reconnaissance de formes. CR Acad Sci Paris, Série I, No 321. 1995:1031–1034.
36. Trouvé A. Diffeomorphisms groups and pattern matching in image analysis. Int J Computer Vision. 1998;28:213–221.
37. Vaillant M, Davatzikos C. Finding parametric representations of the cortical sulci using an active contour model. Medical Image Analysis. 1997;1(4):295–315. [PubMed]
38. Westin CF, Lorigo LM, Faugeras O, Grimson WEL, Dawson S, Norbash A, Kikinis R. Medical Image Computing and Computer-Assisted Intervention—MICCAI 2000, Vol. 1935 of Lecture Notes in Computer Science. Springer-Verlag; 2000. Segmentation by adaptive geodesic active contours; pp. 266–275.
39. Xu C, Pham DL, Prince JL. Finding the brain cortex using fuzzy segmentation, isosurfaces and deformable surface models. XVth Int. Conf. on info Proc. in Medical Imaging; June 1997.
40. Xu C, Pham DL, Prince JL. Medical image segmentation using deformable models. In: Fitzpatrick J, Sonka M, editors. SPIE Handbook on Medical Imaging - Volume III: Medical Image Analysis. SPIE; Bellingham, WA: 2000. pp. 129–174.
41. Xu C, Prince JL. Gradient vector flow: A new external force for snakes. CVRP. 1997 Nov;
42. Xu C, Prince JL. Gradient vector flow deformable models. In: Bankman I, editor. Handbook of Medical Imaging. Academic Press; San Diego, CA: 2000.
43. Yezzi A, Tsai A, Willsky A. A fully global approach to image segmentation via coupled curve evolution equations. Journal of Visual Communication and Image Representation. 2002;13(1–2):195–216.