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 = 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 .

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∘φt−1, t [0, 1], I. 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 Ω d. 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 v which depends on time and space, such that

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

and (0, x) = x for all x Ω. We will also use the notation
φtv(x) for (t, x) and

φsv−1=φtv∘φsv−1.

(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 is given, the element of which being smooth enough vector fields on Ω, and denote the norm and inner product on this space by || ||, ,. 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], ). We will always assume that can be embedded in the space of (
C01(Ω,Rd), || ||1, ∞), containing vector fields on Ω, which vanish on Ω, where

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

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

Fixing v, one can define the linear form wv, wg, which will be denoted Lv. We therefore have the identity

(Lv,w)≐〈v,w〉g

(we use the standard notation (M, w) for the linear form M applied to w). By definition, Lv belongs to the dual, * of , and L can be seen as an operator L: → * (this is the canonical duality operator of 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 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 *, under our running assumption that is embedded in the space of C1 functions:

L1 vector fields of Ω: if ψ: Ω d is integrable, define

(M,v)=∫Ω〈ψ(x),v(x)〉Rddx

Let now μ be any measure on Ω, and ψ be μ integrable. Define

(M,v)=∫Ω〈ψ(x),v(x)〉Rddμ(x)

Dirac measures: as a particular case of the previous, define, for x Ω and ad,

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 vL1([0, 1], ), and let (
φ0tv, t [0, 1]) be the associated flow, defined in the previous section. Along time, each point x Ω, 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 Ω, 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, yvt(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)=12∫01||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 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* 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(φt−1(φ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 ), denoted Adv. We have obtained the relation

vt=Adφtvtl.

To interpret the adjoint action pictorially, the new vector field under the adjoint action v → (d)v−1 has to be interpreted as the transformation of v under the deformation generated by . Figure 1 shows how the field vl at location x is transported by the flow to the value v(y) at location y = (x) by pushing forward (using ) 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)).

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φt∗Lvt,vtl).

This leads to the definition
Mtl=Adφt∗Lvt 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 Gd() be a Lie group of d × d matrices with Lie algebra . 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, is the vector space of antisymmetric 3 × 3 matrices. Let tgt be a trajectory in this group. Then, the angular velocity at 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: → *. 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 and *). The total energy is

E(g)=12∫01(Jat,at)dt.

We retrieve again the analogy with the diffeomorphisms by letting

〈a,b〉g=(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 is [a, b] = ab − ba.

Theorem 1

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

∂Ja∂t−ada∗(Ja)=0.

(2)

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

Proof

Let (tg0(t)) be an extremal curve for the kinetic energy and ((t, h) 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

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

∂Ja∂t−ada∗(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
Mtb≐Jat and the Euler equation is then:

∂Mb∂t−ada∗(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:

Mtb≐Adgt−1∗(Mtb)

(8)

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

Proof

Thus, from the conservation of the momentum to the space,
Mts≡M0s, we deduce that

Jat=Adgt∗(La0),

(11)

or equivalenty, for any b:

(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 and the associated adjoint action ad of 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,(adv∗f,w)=(f,advw).

(14)

with G, w, f*.

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 , it cannot be guaranteed that it is still so for [w, v] = (dw)v − (dv)w: in situations of interest, is, in fact, not a Lie algebra: Adw and advw do not necessarily belong to . As a consequence, the definition of
adv∗f 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 . 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 LvH (i.e. it is a function), one has

adv∗Lv=div(Lv⊗v)+dv∗Lv.

(16)

where div(uv) = 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 (tt) be extremal and ((t, ε) t,ε) be a smooth deformation around ε = 0, with the abuse of notation t,0 = 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

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

(adv∗Lv,w)=(Lv,dvw−dwv)=(dv∗Lv,w)−(Lv,dwv)

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

(div(Lv⊗v),w)=−(Lv,dw.v).

It appears that the Euler equation (15) with
adv∗Lv=divvLv+dv∗Lv 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 MtLvt; 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,

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

(18)

Formal Derivation

Indeed, fix w 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′= Ad Ad′. Now, if 0 = id and
dφεdε=v at ε = 0, we have for any w′,

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 [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
∫01〈vt,ht〉Ldt=∫01(Lvt,ht)dt=∫01(Lv0,dφ0tφt0vt∘φ0t)dt. Now, using (20) and the fact that d0tt1 = d0001d0tt0, we get easily that
∫01dφ0tφt0vt∘φ0tdt=0 so that, by linearity,

∫01〈vt,ht〉Ldt=(Lv0,∫01dφ0tφt0ht∘φ0tdt)=0.

(21)

4.3. Coadjoint Transport of Structures Along a Geodesic

For M*, the evolution
t↦Adφt−1∗M 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*. It is defined by the complementary of the union of all open sets Ω′ Ω which are such that (M, w) = 0 whenever w vanishes outside Ω′. We have the property:

Statement 3

If
Mt=Adφt0∗(M0), then

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

Indeed, assume that M0 vanishes Ω′ Ω. Let w have its support included in 0t (Ω′). Then (Mt, w) = (M0, d0t)−1w0t and w0t vanishes outside Ω′, which implies that (Mt, w) = 0. Thus Supp(Mt) 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 *, therefore have the form Lv0 for some v0. So, we assume that

(M0,w)=∑i〈ai,w(xi)〉Rd,

(22)

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

δx∗a:g→Rw↦〈a,w(x)〉

(23)

Denoting xi (t) 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,w〉dμ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) = ∫Ων0, wdμ0 then the linear momentum functional evolves according to

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 C1 sub-manifold in Ω d.

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 α0dσ0 for some density α0 (not necessary positive) on the surface. Let ν0: Ω → d (the values of ν0 outside Σ0 will not be important) and define

(Lv0,w)=∫∑0〈ν0,w〉Rdα0dσ0.

(29)

Using Statements 3 and 4, we get that the tranported measure μt is located on the transported sub-manifold Σtt (Σ0) (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(t0)*γ0t0, 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 = α0t,0|dt,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.

where
Nt(x)=Span{∇fti(x)∣1≤i≤r} 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 ψG such that I0ψ−1 = It. By compactness and semi-continuity arguments, one can prove the existence of a geodesic path = (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, t(x) x + tv0(x) so that
It=I0∘φt−1(x)≃I0−t〈∇I0,v0〉Rd. By definition, the cost to go from I0 to It is (still at first order) t|v0|L. However, any u such that I0, ud = I0, v0d, 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):infu∈g∣u∣Lsubjectto:〈(v0−u)(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 I0 = {h | I0(x), h(x)d = 0, ∀x Ω}: the constraints in can be restated as u − v0I0 so that pI0 (v) is the orthogonal projection of v on
gI0⊥, the space orthogonal to I0. Hence, equality (36), translates to

∀h∈gsuchthatforallx∈Ω,〈∇I0(x),h(x)〉=0,wehave〈v0,h〉L=0.

(37)

Now, the fact that I0(x), h(x) 0 means that h is a vector field which is tangent to the level sets of I0, and since v0, hL = (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 I0. Obviously, when I0 is smooth, hI0 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 Ω is smooth enough and take advantage on the fact that elements of vanish on Ω), we get

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

Since this has a meaning when I0L2(Ω), we now define

Definition 2

When IL2(Ω), we denote

gI={h∈g:〈I,div(hf)〉L2=0forallf∈C1(Ω)}.

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 is continuously embedded in C1(Ω, d)).

Theorem 3. (Normal Momentum Constraint)

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

5.3. Examples

Consider again the case of smooth I, so that the condition hI is equivalent to h and for all x Ω, x I, h(x)d = 0. Using notation (23), we get that
〈∇I0(x),h(x)〉Rd=〈δx∗(∇I0(x)),h〉L 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 such that, for hV,ωf (I), h = I, div(h f)L2(Ω), and conclude that

gI⊥=Span{ωf(I)∣f∈C1(Ω,R)}¯.

(40)

This implies that any element vI is such that Lv can be expressed as a limit

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

where f NC1(Ω, ). The particular case when I is the indicator function of a smooth domain B Ω (which can be interpreted as a smooth shape) is quite interesting. For xB, let ν(x) be the outward normal to Ω and denote σB be the uniform measure on 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 I by considering elements v such that

(Lv,u)=∫∂B〈v(x),u(x)〉Rddμ(x).

(41)

for some measure μ on 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 vtIt 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 I 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 x I ≠ = 0, we can de-fine in a unique way a positively oriented orthonormal frame (τ (x) ν (x)) such that ν (x) = x I/|x I|. Then, if hI and h(x) ≠ = 0, we must have h/|h| = ± τ in a small ball around x. Now h, as an element of 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 , which may very well be everywhere on Ω (or on a dense subset, which is the same since h is continous), in which case I = {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 vI 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, and t is the diffeomorphism generated by a geodesic, there is a priori no reason for (dt)−1wt) to belong to : one must be able to define Lv0 on spaces which are bigger than , which means that Lv0 needs to be somewhat smoother (as a distribution) than a generic element of *.

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 (, ρ), together with two measurable functions x, y: → Ω are given. The diffeomorphism is searched to minimize

U(φ)=E(φ)+1σ2∫J∣yi−φ(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σ2∫J∣yi−φ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
δx∗a the linear form on such that
(δx∗a,v)=〈a,v(x)〉. Let v be a minimizer of Ũ. Then, letting
xi(t)=φ0tv(xi)

Lvt=−1σ2∫Jδxi(t)∗[(dxi(t)φt1v)∗(yi−xi(1))]dρ(i)

(42)

Proof

The proof of this result is a direct consequence of the identity, valid for s, t [0, 1], v, hL2([0, 1], ),

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
δx∗a, for x Ω and ad.

5.4.2. Inexact Shape Matching

We now consider the comparison of a binary set-indicator function, I0 = 1Ω0 (0 Ω 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], ). Then

Lvt=1σ2∫∂Ω1(12−I1)δφ1tv(x)∗[(dφ1tvφt1v)∗ν1]dσ1(x)

(44)

where
Ω1=φ01v(Ω0), ν1 is the outward normal to Ω1 and σ1 is the volume measure on Ω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)(12−I1(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+εh∣dx

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+εh∣dx=∫Ω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(12−I1)〈(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(12−I1∘φ01v(x))〈v0,w〉dμ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σ2∣dφt,1v∣(I0tv−I1tv)∇I0tv

(45)

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

Lv0=−1σ2∣dφ01v∣(I0−I10v)∇I0

as can be deduced from Eq. (30), with
α=−1σ2∣dφ01v∣(I0−I10v). 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 I0t,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 (t, t [0, 1]) to act on the elements in the orbits of a given template I = Itemp. Now we examine the geodesic flows in the orbit {It = It, t [0, 1]}, I, and provide the associated evolution equations.

6.1. Geodesic Evolution Equation of Landmark Points

Here we examine the finite dimensional landmark orbit denoted n, consisting of n-shapes IN = (x1, …, xn), each landmark (xi)1≤i≤n is in Ω d; correspondingly (ai)1≤i≤n are a finite family of vectors in d. Denoting by xi (t) =̇ t (xi), the trajectory in Ω of the point xi under the flow t gives

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

(46)

where ai (t) = (dxt (t) 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

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 σx(σ) for σ [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

∂νt∂t(t,σ)+(dx(t,σ)vt)∗ν(t,σ)=0,

∂x∂t(t,σ)−vt(x(t,σ)=0.

6.2. Geodesic Image Evolution

Assume here that dμ = α (x)dx has a density with respect to Lebesgue’s measure on Ω. In this case, dμt = |dt0|αt0dx and

Lvt=∣dφt0∣α∘φt0∇It

(50)

From the conservation of momentum in Lagrangian coordinates for image based motion, we get for Lv0 = α0I0 that Lvt = αt |dt,0|It where αt = α0t,0, It = I0t,0. Let zt = αt |dt,0| so that Lvt = ztIt.

Since
ddt(zt∘φ0,t)=ddt(∣dφ0,tφt,0∣α0)=−div(vt)∘φ0,t∣dφ0,tφt,0∣α0 we get
dztdt+dzt(vt)+div(vt)zt. Moreover, we get easily
∂It∂t+〈∇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(ztIt):

dztdt+dzt(vt)+div(vt)zt=0,

∂It∂t+〈∇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

7. Computational Results

The following results illustrate the computation of the momentum Lv0 = α0I0 (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 ( + 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.

Shown in Fig. 4–6 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.

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 α0I0, 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. 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.

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 ht =̇ vt − pIt(vt). For ε [0, 1], let
vtε≐vt+εht, and
φtε≐φtvε (one can check, but we skip the argument, that t → ht is measurable and belongs to L2([0, 1], ), 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 ht, vtL = vt − pIt(vt), vtL = 0, we get
∣vt+ht∣L2=∣vt∣L2−∣ht∣L2. Since t → vt corresponds to paths with lowest kinetic energy from I0

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∫ΩI0∘qtε(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 (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 {I ≤ i}, v(x) is normal to Ωi if xΩ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 Ω, the bilinear form Kx,y on d × d defined by