Reduction by Lie Group Symmetries in Diffeomorphic Image Registration and Deformation Modelling
Abstract
: We survey the role of reduction by symmetry in the large deformation diffeomorphic metric mapping framework for registration of a variety of data types (landmarks, curves, surfaces, images and higher-order derivative data). Particle relabelling symmetry allows the equations of motion to be reduced to the Lie algebra allowing the equations to be written purely in terms of the Eulerian velocity field. As a second use of symmetry, the infinite dimensional problem of finding correspondences between objects can be reduced for a range of concrete data types, resulting in compact representations of shape and spatial structure. Using reduction by symmetry, we describe these models in a common theoretical framework that draws on links between the registration problem and geometric mechanics. We outline these constructions and further cases where reduction by symmetry promises new approaches to the registration of complex data types.1. Introduction
Registration, the task of establishing correspondences between multiple instances of objects, such as images, landmarks, curves and surfaces, plays a fundamental role in a range of computer vision applications, including shape modelling [1], motion compensation and optical flow [2], remote sensing [3] and medical imaging [4]. In the subfield of computational anatomy, establishing inter-subject correspondences between organs allows the statistical study of organ shape and shape variability [5]. Examples of the fundamental role of registration include quantifying developing Alzheimer’s disease by establishing correspondences between brain tissue at different stages of the disease [6]; measuring the effect of chronic obstructive pulmonary disease on lung tissue after removing variability caused by the respiratory process [7]; and correlating the shape of the hippocampus to schizophrenia after inter-subject registration [8].
In this paper, we survey the role of reduction by symmetry in diffeomorphic registration and deformation modelling, linking symmetry as seen from the field of geometric mechanics with the image registration problem. All of our calculations will be formal and void of functional analytic detail, although citations will be used when available. We focus on large deformations modelled in subgroups of the group of diffeomorphic mappings on the spatial domain in the context of large deformation diffeomorphic metric mapping (LDDMM) [1,9–11]. Connections with geometric mechanics [12,13] have highlighted the role of symmetry, and properties that were previously known to be connected with the specific data types have been described in a common theoretical framework [14]. We wish to describe these connections in a form that highlights the role of symmetry and points towards future applications of the ideas.
1.1. Symmetry and Information
One of the main reasons symmetry is useful in data analysis and numerics is its ability to reduce the complexity of information that represents data. Lower information complexity can lead to more stable statistical analysis and the reduced need of computational resources.
As a toy example, consider a spinning top. Upon choosing a reference configuration, the orientation of the top is given by a rotation matrix, i.e., an element R ∈ SO(3) (see Figure 1). To describe the direction of the tip of the top, it suffices to provide the orientation matrix R. However, R is contained in SO(3), a three-dimensional space, while the space of possible directions is the two-sphere, S2, which is only of dimension two. Therefore, providing the full matrix R is an over-representation of the tip direction. It suffices to solely provide the vector R · k ∈ S2 where k = (0, 0, 1) is the direction of the tip in a reference configuration. Note that if , then . Therefore, given only the direction k′ = R · k, we can only reconstruct R up to an element , which preserves k. The subgroup of rotations that preserve k can be identified with SO(2). Specifically, this identification comes from perceiving a rotation about k as a rotation of the plane, which is perpendicular to k. This insight allows us to express the space of directions S2 as a homogeneous space S2 ≡ SO(3)/SO(2). In terms of information, we can cartoonishly express this by:
This picture is typical for many group quotients. Generally speaking, if X is a manifold and G acts freely and properly on X, then:
When X is infinite dimensional, this formula is less insightful. However, X/G is smaller than X in that there exists a surjective map from X to X/G with non-trivial level sets. Reduction by symmetry can be implemented when a problem posed on X has G symmetry and can be rewritten as a problem posed on X/G. As X/G is smaller than X, reduction by symmetry can yield more stable subsequent statistical analysis of observed data and more tractable algorithms, as will be shown later in the article. This reduction is particularly dramatic when dim(X) = ∞ and dim(X/G)<∞.
1.2. Symmetry in Registration
Registration of objects contained in a spatial domain, e.g., the volume to be imaged by a scanner, can be formulated as the search for a deformation that transforms both domain and objects to establish an inter-object match. The data available when solving a registration problem is generally insufficient for determining the displacement of every point of the domain. This is the case when images to be matched have areas of constant intensity and no derivative information can guide the registration. For example, the “best” deformation for matching the two discs in Figure 2 is ambiguous, except at the boundary of the discs, where the images are non-constant. Similarly, when 3D shapes are matched based on the similarity of their surfaces, the deformation of the interior cannot be derived from the available information. In these cases, the deformation model is over-complete, and a range of deformations can provide equally adequate matches for the data. The registration problem or the registration cost-function is thus symmetric with respect to the subset of transformations just described. When the deformation model is a Lie group, the deformations for which the registration is symmetric form a subgroup. The quotient by this subgroup of symmetries of the registration cost-function can provide vastly more compact representations. This situation arises in several cases with the LDDMM framework: when registering images, only displacements orthogonal to the level lines of an image are needed, and when registering shapes, the information left in the quotient is supported on the shape surface only.
1.3. Content and Outline
Although a degree of comfort with differential geometry will be assumed, it is the aim of this paper to make the role of symmetry in registration and deformation modelling clear to the non-expert in geometric mechanics and symmetry groups in image registration. We begin the paper by presenting the background for the registration problem and the large deformation approach before outlining some necessary notions from differential geometry. For more information on the Riemannian geometry behind the LDDMM approach to image registration, we refer the reader to [5]. We continue by describing how reduction by symmetry leads to an Eulerian formulation of the equations of motion when reducing to the Lie algebra. The symmetry of dissimilarity measures allows additional reductions, and we use isotropy subgroups to reduce the complexity of the registration problem further. Lastly, we survey the effect of symmetry in a range of concrete registration problems. The paper ends with concluding remarks.
2. Registration
The registration problem consists of finding correspondences between objects that are typically point sets (landmarks), curves, surfaces, images or more complicated spatially-dependent data, such as diffusion weighted images (DWI). The problem can be approached by letting M be a spatial domain containing the objects to be registered. M can be a compact finite dimensional differentiable manifold without boundary or ℝd itself with d = 2, 3. It is common to consider manifolds with boundaries, as well. In such cases, care must be taken with regards to boundary conditions. For example, vector-fields must be tangential to the boundary. Here, we consider only manifolds without a boundary.
A map φ : M → M can deform or warp the domain by mapping each x ∈ M to φ(x). The deformation encoded in the warp will apply to the objects in M, as well as the domain itself. For example, if the objects to be registered consist of point sets {x1,…, xN}, xi ∈ M, the set will be mapped to {φ(x1),…, φ(xN)}. For surfaces S ⊂ M, φ similarly results in the warped surface φ(S). Because those operations are associative, the mapping φ acts on {xi} or S, and we write φ · {xi} and φ · S for the warped objects. An image is a function I : M → ℝ, and φ acts on I as well, in this case by composition with its inverse φ · I = I ○ φ−1; see Figure 3. For this, φ must be is invertible, and commonly, we restrict to the set of invertible and differentiable mappings Diff(M). For various other types of data objects, the action of a warp on objects can be defined in a manner similar to that of point sets, surfaces and images. This fact relates a range registration problems to the common case of finding appropriate warps φ, which bring the objects into correspondence. Different shape instances can be realized by letting warps act on a base shape, and a class of shape models can thereby be obtained by using deformations as shape representations [1].
The search for appropriate warps can be formulated in a variational formulation with an energy:
The regularization term can take various forms often modelling physical properties, such as elasticity [15], and derivatives of φ are often penalized to ensure smoothness. For some choices of R, existence and analytical properties of minimizers of Equation (1) have been derived [16]; however, in general, it is difficult to ensure that solutions are diffeomorphic by penalizing φ in itself. The free-form-deformation (FFD; [17]) and related approaches model the deformation by a displacement vector field u on M = ℝd, so that φ(x) = x + u(x). Smoothness is here ensured by the choice of basis functions, e.g., B-splines, or by applying a regularization term on u. Smooth and invertible mappings can be obtained by integrating flows [9,11] to obtain one-parameter families or paths of mappings φt, t ∈ [0, 1]. The warp φ0 at t = 0 is here the identity mapping id ∈ Diff(M), and the dissimilarity is measured at the endpoint φ1. The time evolution of φt can be described by the differential equation:
Here, the time-dependent diffeomorphism φt is related to ut through Equation (2). If the norm ‖·‖V that measures the irregularity of ut is sufficiently strong (e.g., Hk with ), then φt will be a diffeomorphism for all t [19]. This approach thus gives a direct way of enforcing properties of the generated warp: instead of regularizing φ directly, the analysis is lifted to a normed space V that is much easier control. The energy E in Equation (3) has the same minimizers as the geometric formulation of LDDMM used in the next section.
Direct approaches to solving the optimization problem in Equation (3) must handle the fact that the problem of finding a warp is now expanded to that of finding a time-dependent family of warps. This is a huge increase in dimensionality. This formulation of registration is thus very difficult to represent numerically and to optimize and analyse statistically. For several data types, it has been shown how optimal paths for Equation (3) have specific properties that reduce the dimensionality of the problem, making practical solutions feasible. In the next section, we outline the geometric framework that is needed when we, in the later sections, use reduction by symmetry to describe these data-dependent results in a common theoretical framework.
3. Notions from Differential Geometry
In this section, we will introduce a number of notions from differential geometry in a fairly informal manner. We will use conventions from [20,21] where a more rigorous understanding of differential geometry can be found.
We will assume that the reader has at least an intuitive picture of the notion of a smooth manifold M. For the purpose of this paper, M will either be assumed to be compact without a boundary or ℝn. The tangent bundle of M is the space of velocity vectors tangential to M. Notationally, the tangent bundle of M, denoted TM, is the set of pairs (x, v) where x ∈ M and v is a vector tangential to M at the point x (see Figure 4). A vector-field is a continuous map u : M → TM, such that u(x) ∈ TM is a vector above x for all x ∈ M. We will denote the space of vector-fields by .
Given a vector-field , we may consider the initial value problem:
Conversely, let φt ∈ Diff(M) be a time-dependent diffeomorphism. For any x ∈ M, we observe that φt(x) is a curve in M. If this curve is differentiable we may consider its time-derivative, , which is a vector above the point φt(x) ∈ M. From these observations, it follows that is a vector above x. Therefore, the map ut : M → TM, given by , is a vector-field called the Eulerian velocity field of φt.
As will be described shortly, the Eulerian velocity field contains less data than . This reduction in data can be viewed from the perspective of symmetry. Given any ψ ∈ Diff(M), the curve φt can be transformed to the curve φt ○ ψ. We observe that:
Thus, φt and φt ○ ψ both have the same Eulerian velocity fields. In other words, the Eulerian velocity field, ut, is invariant under particle relabellings. More precisely, we may view Diff(M) as a manifold in its own right, and view as a vector in the infinite-dimensional tangent bundle T Diff(M) above the “point” φt ∈ Diff(M). Thus, the vector contains both velocities and a base diffeomorphism φt. Given ut and φt, we can reconstruct via As has been shown, we can also construct ut from by its own definition. However, we cannot reconstruct from ut, which is why ut contains less data.
Finally, we will denote some linear operators on the space of vector-fields. Let Φ ∈ Diff(M), and let . The push-forward of u by Φ is the vector-field given by:
In local coordinates (x1,…, xn), this looks like:
By inspection, we see that Φ∗ is a linear operator on . One can view Φ∗u as “u in a new coordinate system”, because any geometric property of u is also inherited by Φ∗u. For example, if S is invariant under u, then Φ(S) is invariant under Φ∗u. We define the pull-back by Φ∗u := (Φ−1)∗u. Note that Φ∗(Φ∗(u)) = u.
As Φ∗ is a linear operator, a well-defined operator exists, which is dual to Φ∗. Let denote the dual space to , i.e., the set of linear maps , which are continuous with respect to a chosen vector-space topology on . Given , we define by the identity:
In local coordinates, we may represent m as a one-form density, given by mi(x)dxi⊗(dx1∧…∧dxn) with components m1(x),…, mn(x). In this local coordinate description, the i-th component of the push-forward looks like:
Finally, we define the Lie derivative. Let . The Lie derivative operator, with respect to w, is the linear operator defined by:
Again, as £w is a linear operator on , we can define a dual operator on . If , we can define £w[m] ∈ X(M)∗ by the equation:
We conclude the section with a table of notation for the reader’s convenience, see Table 1.
4. Reduction by Symmetry in LDDMM
In this section, we will present necessary conditions satisfied by local extremizers of the variational problem Equation (3). The resulting conditions will first involve the computation of a curve in Diff(M), as well as its time-derivative in T Diff(M). We then invoke a Diff(M) symmetry of the problem to reduce this computation to a computation on instead of T Diff(M) ≅ Diff(M) ×. Secondly, we describe how the symmetry of the dissimilarity measure allows further reductions.
4.1. Reduction to the Lie Algebra
The variational formulation Equation (3) of LDDMM is equivalent to minimizing the energy:
Example 1. Given images I0, I1 ∈ L2(M), we consider the dissimilarity measure:
In this article, we will consider the metric on connected components of Diff(M) given by:
Given P, minimizers of E Equation (6) must satisfy:
That Equation (7) is a necessary condition satisfied by the minimizers of Equation (6) (under certain analytic assumptions) is stated and proven in [1] (Proposition 11.6) for the case of M = ℝn. The proof is based on the definition of the curve energy Equation (6) as a sum a Riemannian distance d and the dissimilarity measure F. Minimization of the Riemannian distance d yields a search for a geodesic in Diff(M), because geodesics are distance minimizing. The corresponding geodesic equations are given by the first two lines of Equation (7). The term F only penalizes the end-point of the geodesic, and the minimization condition manifests as the third line of Equation (7).
Issues regarding the well-posedness of Equation (7) are non-trivial, because P is merely injective, but not bijective, and so, there is no guarantee that P can be inverted on a given at each time in order to obtain a vector-field . Fortunately, safety guards for well-posedness are known (e.g., [19], Theorem 1, or [25]).
Using Equation (7) for computational purposes is difficult because Diff(M) is a non-linear infinite dimensional space. Moreover, the dissimilarity measure F only comes into play at time t = 1, and the distance function is an integral over the vector-space . It would be beneficial if we could rewrite the extremizers in terms of the Eulerian velocity field u and the flow at t = 1. In fact, this is often the case. One (formally) must take the time-derivative of the term “ ” and apply Equation (5). Explicitly, this computation is performed as follows. Let , and observe:
As w is arbitrary, we find ∂tmt + £ut[mt] = 0. This allows us to reformulate Equation (7) as:
The advantage of this formulation is that the bulk of the computation occurs on the vector-space rather than on the space T Diff(M). Registration algorithms based on Equation (8) differ from the algorithm proposed by Beg et al. in [26]. In [26], a gradient descent on the time-dependent Eulerian vector field ut is used to minimize the curve energy. The algorithm is posed on the velocity field ut instead of the momentum field mt as Equation (8) suggests. The momentum at time t is retrieved by a transport equation similar to the first equation in Equation (7). The evolution Equation (8) effectively allows one to search over the space of initial conditions, , rather than over the larger space of curves, C1([0, 1]; ).
This reduction of the problem from T Diff(M) given in Equation (7) to the problem over space of vector-fields, given in Equation (8) is the first instance of reduction by symmetry. In particular, this corresponds to the fact that the space of vector-fields is identifiable as a quotient space:
Additionally, the map is the quotient projection.
4.2. Isotropy Subgroups
The reduction to dynamics on T Diff(M) to dynamics on occurs primarily because the distance function is Diff(M) invariant. However, one cannot completely abandon Diff(M), because the solution requires one to compute the Time 1 flow, . Fortunately, there is a second reduction, which allows us to avoid computing in its entirety. This second reduction corresponds to the invariance properties of the dissimilarity measure F. Let GF ⊂ Diff(M) denote the set of diffeomorphisms, which leave F invariant, i.e.:
One can readily verify that GF is a subgroup of Diff(M), and so, we call GF the isotropy subgroup of F.
Having defined GF, we can now consider the homogeneous space Q = Diff(M)/GF, which is the quotient space induced by the action of the right composition of GF on Diff(M). This quotient space is “smaller” in the sense of the amount of data required to describe an element of it. In terms of maps, this can be seen by defining the map φ ∈ Diff(M) ↦q = [φ]/GF ∈ Q, where [φ]/GF denotes the equivalence class of φ. We call this mapping the quotient projection, because it sends Diff(M) to Q surjectively. While these notions are theoretically quite complicated, often they manifest less so in practice.
Example 2. In this example, we consider a simple aspect of the landmark matching problem. Let M ⊂ ℝn be the closure of some open set. Let x1, x2, y1, y2 ∈ M with x1 ≠ x2, and consider the dissimilarity measure:
We see that:
Note that Diff(M) is infinite dimensional, while Q is of dimension 2 dim(M). Two examples of diffeomorphisms contained in GF are shown in Figure 5.
Example 3. In this example, we consider the matching problem for greyscale images. Let I0, I1 ∈ Hk(M) be images. There is a natural action of Diff(M) on Hk given by sending each image I ∈ Hk(M) to I ○ φ ∈ Hk(M). We could consider the matching function F : Diff(M) → ℝ given by:
This function measures the difference between a deformed version of I0 and a fixed image, I1. The isotropy subgroup GF is the group of images that preserve 00. Such a diffeomorphism would preserve each of the level lets of I0, but could permute the points within a given level set.
If one is able to understand Q, then one can use this insight to reformulate the dissimilarity measure F as a function on Q, rather than Diff(M). In particular, there exists a unique function FQ : Q → ℝ defined by the property FQ([φ]/GF ) = F(φ). Again, this is useful in the sense of data, as illustrated in the following example.
Example 4. Consider the dissimilarity measure F of Example 2. The function, FQ : Q → ℝ is:
Finally, note that Diff(M) acts upon Q by the left action:
Usually, we will simply write ψ · q for the action of ψ ∈ Diff(M) on a given q ∈ Q. This means that acts upon Q infinitesimally, as it is the Lie algebra of Diff(M).
Example 5. Consider the setup of Example 2. Here, Q = M × M − ΔM×M, and the left action of Diff(M) is given by:
These constructions allow us to rephrase the initial optimization problem using a reduced curve energy. Minimization of E is equivalent to minimization of:
A geometric derivation of this formula can be found in [13] (Lemma 2.8, (2.12), and (2.13)). Again, the solution only depends on the Eulerian velocity and q(1). For this reason, we see that the GF symmetry of F provides a second reduction in the data needed to solve our original problem.
4.3. Orthogonality
In addition to reducing the amount of data that we must keep track of, there is an additional consequence to the GF-symmetry of F. In particular, there is a potentially massive constraint satisfied by the Eulerian velocity u.
To describe this, we must introduce an isotropy algebra. Given , we can define the (time-dependent) isotropy algebra:
This is nothing, but the “Lie algebra” associated with the isotropy group Gq(t) = {ψ ∈ Diff(M) | ψ·q(t) = q(t)}. The use of quotes here is deliberate. If we let denote an RKHS obtained from the space of vector-fields, then some of these are permitted to be non-smooth, which means that the standard Lie bracket of vector-fields does not close.
It turns out that the velocity field ut that minimizes E (or EQ) is orthogonal to with respect to the chosen inner-product. Intuitively, this is quite sensible, because velocities that do not change q(t) do not alter the data and simply waste control effort. Equivalently, said from the perspective Lagrange-multipliers, we know that the Lagrange-multipliers used to enforce optimality should point in a direction orthogonal to those which leave the cost functional unaltered. These variations on the same statement are formalized below.
Proposition 1. Let u satisfy Equation (8) or Equation (9). Then, m = P[u] annihilates.
Proof. Let u be the solution to Equation (9). We will first prove that u1 is orthogonal to . Let . We observe:
However, w1 leaves q(1) fixed, so . Therefore, 〈P [u1], w1〉 = 0. Let . In coordinates, this means:
One can directly verify that for all t ∈ [0, 1]. Denoting mt = P [ut], as in Equation (9), we find:
The last equality follows from Equation (5). Thus, 〈P [ut], wt〉 is constant. We have already verified that at t = 1, this inner-product is zero, and therefore, 〈P [ut], wt〉 = 0 for all time. That w1 is an arbitrary element of makes wt an arbitrary element of at each time. Thus, ut is orthogonal to for all time. □
At this point, we should return to our example to illustrate this idea.
Example 6. Again, consider the setup of Example 2. In this case, q(t) = (q1(t), q2(t)) ∈ M × M − ΔM×M. The spaceis the space of vector-fields that vanish at q1(t) and q2(t). Therefore, ut is orthogonal to q(t) if and only if mt = P [ut] satisfies:
This orthogonality constraint allows one to reduce the evolution equation on to an evolution equation on Q, which might be finite dimensional if GF is large enough. In particular, there is a horizontal lift, , uniquely defined by the conditions , and with respect to the chosen inner-product on vector-fields.
Example 7. Consider the setup of Example 2 with M = ℝn. Then, . Let K : ℝn×ℝn→ℝn×n be the reproducing kernel of P. Then, is given by:
One can immediately observe that h↑ is injective and linear in . In other words, is an injective linear map for fixed q ∈ Q. Because the optimal ut is orthogonal to , we may invert h↑(q(t), ·) on ut. In particular, we may often write the equation of motion on TQ, rather than on . This is a massive reduction if Q is finite dimensional. In particular, the inner-product structure on induces a Riemannian metric on Q given by:
The equations of motion in Equations (8) and (9) map to the geodesic equations on Q.
Proposition 2. Let u extremize E or EQ. Then, there exists a unique trajectory q(t) ∈ Q, such that. Moreover, q(t) is a geodesic with respect to the metric g.
Proof. Let u minimize E. Thus, u satisfies Equation (9). By the previous proposition, ut is orthogonal to . As is injective on , there exists a unique , such that Note that E can be written as:
Thus, minimizers of E correspond to geodesics in Q with respect to the metric g. □
If we let H: T∗Q → ℝ be the Hamiltonian induced by the metric on Q, we obtain the most data-efficient form or Equations (8) and (9). Minimizers of E (or EQ) are:
We see that this is a boundary value problem posed entirely on Q. If Q is finite dimensional, this is a massive reduction in terms of data requirements.
Example 8. Consider the setup of Example 2 with M = ℝn. The metric onis most easily expressed on the cotangent bundle T∗Q. If K is the matrix valued kernel of P, the metric on T∗Q takes the form:
4.4. Descending Group Action
A related approach to defining distances on a space of objects to be registered consists of defining an object space upon which Diff(M) acts transitively (this means that for any o1, o2 ∈ , there exists a φ ∈ Diff(M) such that φ·o1 = o2) with distance:
Here, the distance on is defined directly from the distance in the group acting on the objects; see for example [1,5]. With this approach, the Riemannian metric descends from Diff(M) to a Riemannian metric on , and geodesics on lift by horizontality to geodesics on Diff(M). The quotient spaces Q obtained by reduction by symmetry and their geometric structure correspond to the object spaces and geometries defined with this approach. Intuitively, reduction by symmetry can be considered a removal of redundant information to obtain compact representations, while letting the metric descend to the object space constitutes an approach to defining a geometric structure on an already known space of objects. The resulting solutions are equivalent to the ones presented in this article, because ≅ Diff(M)/Go, where Go = {ψ ∈ Diff(M) | ψ(o) = o} for some fixed reference object o ∈ .
5. Examples
Here, we present a number of concrete examples of how reduction by symmetry can reduce the infinite dimensional registration problem over Diff(M) to lower, in some cases finite, dimensional problems. In all examples, the symmetry of the dissimilarity measure with respect to a subgroup of Diff(M) gives a reduced space by quotienting out the isotropy subgroups.
5.1. Landmark Matching
The space Q used in the examples in Section 4 constitutes a special case of the landmark matching problem, where sets of landmarks Q = {(x1,…, xN)| xi ∈ M, xi ≠ xj ∀i ≠ j} are placed into spatial correspondence trough the left action φ·(x1,…, xN) = (φ(x1),…, φ(xN)} of Diff(M) by minimizing the dissimilarity measure . The landmark space Q arises as a quotient of Diff(M) from the isotropy group GF, as in Example 2.
Reduction from Diff(M) to Q in the landmark case has been used in a series of papers starting with [27]. Hamilton’s equations (Equation (10)) take the form:
5.2. Curve and Surface Matching
The space of smooth non-intersecting closed parametrized curves in ℝn is also known as the space of embeddings, denoted Emb(S1, ℝn). The parametrization can be removed by considering the right action of Diff(S1) on Emb(S2, ℝn) given by:
Then, the quotient space Gr(S1, ℝn):=Emb(S1, ℝn)/Diff(S1) is the space of unparameterized curves. The space Gr(S1, ℝn) is a special case of a non-linear Grassmannian [29], and it has a manifold structure under certain conditions on the space of embeddings and the space of diffeomorphisms [30].
When the parametrization is not removed, embedded curves and surfaces can be matched with the current dissimilarity measure [31,32]. If M is a volume manifold, then the objects are considered elements of the dual space of Ωk(M), the space of differential k-forms on M. In the surface case, a bounded submanifold S ⊂ M can be seen as an element of [Ωk(M)]∗ by its evaluation on a k-form, w ∈ Ω2(M), given by:
The isotropy groups for curves and surfaces generalize the isotropy groups of landmarks by consisting of warps that keep the objects fixed, i.e.,
The momentum field will be supported on the transported curves/surfaces φ(t).S for optimal paths for E in Diff(M).
5.3. Image Matching
Images can be registered using either the L2-difference defined in Example 1 or with other dissimilarity measures, such as mutual information or the correlation ratio [33,34]. The similarity will be invariant to any infinitesimal deformation orthogonal to the gradient of dissimilarity measure. In the L2 case, this is equivalent to any infinitesimal deformation orthogonal to the level sets of the moving image [35]. The momentum field thus has the form p(t) = α(t)∇φ(t).I0 for a smooth function α(t) on M (see Figure 6), and the registration problem can be reduced to a search over the scalar field α(t) instead of vector field p(t).
Minimizers for E follow the PDE [5]:
In particular, the isotropy group of a source image f0 ∈ C∞ (M) is the subgroup of diffeomorphisms, which preserve the level sets of f0. The quotient space Diff(M)=Iso(f0) can be identified with the orbit of f0 by diffeomorphisms, i.e., Diff(M)=Iso(f0) ≅ Orb(f0):= {φ*f0 | φ ∈ Diff(M)}. This orbit is difficult to identify with a more concrete object, in contrast to, e.g., the case of landmark matching. However, it can be characterized by various properties. For example, for a function f ∈ Orb(f0) and any c ∈ f0 (M) ⊂ ℝ, the level sets f−1(c) and f0−1(c) have the same topology.
5.4. Jet Matching
In [14,36], an extension of the landmark case has been developed where higher-order spatial information is advected with the landmarks. The spaces of jet-particles arise as extensions of the reduced landmark space Q by quotienting out smaller isotropy subgroups known as jet-groups. A thorough account of jet-groups, including Lie group and algebra structures, can be found in [37] (Chapter 4). We provide a brief introduction here. Let G(0) be the isotropy subgroup for a single landmark:
Let now k be a positive integer. For any k-differentiable map f from a neighbourhood of q, the k-jet of f is denoted . In coordinates, consists of the coefficients of the k-th order Taylor expansions of f aboutx. The higher-order isotropy subgroups are then given by:
That is, the elements of G(k) fix the Taylor expansion of the deformation φ up to order k. The definition naturally extends to finite numbers of landmarks, and the quotients Q(k) = G/G(k) can be identified as the sets:
When the dissimilarity measure F is dependent not just on positions, but also on higher-order information around the points, reduction by symmetry implies that optimal solutions for E will be parametrized by k-jets in the same way as Q(0) parametrizesoptimal paths for E in the landmark case. The higher-order jets can thus be used for landmark matching when the dissimilarity measure is dependent on the local geometry around the landmarks. For example, matching of the first order structure, such as image gradients, leads to first-order jets, and matching of local curvature leads to second-order jets.
5.5. Discrete Image Matching
The image matching problem can be discretized by evaluating the L2-difference at a finite number of points. In practice, this always happens when the integral is evaluated at finitely many pixels of the image. In [36,38], it is shown how this reduces the image matching PDE (12) to a finite dimensional system on Q when the integral is approximated by pointwise evaluation at a grid Λh:
In [38], a discretization scheme with higher-order accuracy is in addition introduced with an O(hd+2) approximation F(2) of F. The increased accuracy results in the entire energy E being approximated to order O(hd+2). The solution space in these cases becomes the jet-space Q(2). For a given order of approximation, a corresponding reduction in the number of required discretization points is obtained. The reduction is countered by the increased information encoded in each two-jet. The momentum field thus encodes both local displacement, local linear deformation and second order deformation; see Figure 7. The discrete solutions will converge to solutions of the non-discretized problem as h → 0.
5.6. DWI/DTI Matching
Image matching is invariant with respect to variations parallel to the level lines of the images. With diffusion weighted images (DWI) and the variety of models for the diffusion information (e.g., diffusion tensor imaging (DTI) [39], Gaussian mixture fields [40]), first or higher-order information can be reintroduced into the matching problem. In essence, by letting the dissimilarity measure depend on the diffusion information, the isotropy subgroup of the matching problem becomes smaller.
The exact form of the of DWI matching problem depends on the diffusion model and how Diff(M) acts on the diffusion image. In [41], the diffusion is represented by the principal direction of the diffusion tensor, and the data objects to be matched are thus vector fields. The action by elements of Diff(M) is defined by:
The action rotates the diffusion vector by the Jacobian of the warp, keeping its length fixed. Similar models can be applied to DTI with the preservation of principle direction scheme (PPD, [42,43]) and to GMF-based models [44]. The dependency on the Jacobian matrix implies that a reduced model must carry first order information in a similar fashion to the one-jet space Q(1); however, any irrotational part of the Jacobian can be removed by symmetry. The full effect of this has yet to be explored.
As in the case of image matching, the quotient can be identified with the orbit of the source data under diffeomorphisms.
5.7. Fluid Mechanics
Incidentally, the equation of motion:
In [48], one of the kernel parameters in [46], which controls the compressibility of the u, was taken to the incompressible limit. This allowed a realization of the particle methods described in [47]. The constructions of [48] are the same as presented in this article, but with Diff(M) replaced by the group of volume-preserving diffeomorphisms of ℝd. Velocity fields induced by first order jet-particles are visualized in Figure 8.
6. Discussion and Conclusions
The information available for solving the registration problem is in practice not sufficient for uniquely encoding the deformation between the objects to be registered. Symmetry arises in both particle relabelling symmetry that gives the Eulerian formulation of the equations of motion and in symmetry groups for specific dissimilarity measures.
For landmark matching, reduction by symmetry reduces the infinite dimensional registration problem to a finite dimensional problem on the reduced landmark space Q. For matching curves and surfaces, symmetry implies that the momentum stays concentrated at the curves and surfaces allowing a reduction by the isotropy groups of warps that leave the objects fixed. In image matching, symmetry allows reduction by the group of warps that do not change the level sets of the image. Jet-particles arise from smaller isotropy subgroups and, hence, larger reduced spaces Q(1) and Q(2) that encode locally affine and second order information.
Reduction by symmetry allows these cases to be handled in one theoretical framework. We have surveyed the mathematical construction behind the reduction approach and its relation to the above-mentioned examples. As data complexity rises both in terms of resolution and structure, symmetry will continue to be an important tool for removing redundant information and achieving compact data representations.
Acknowledgments
Henry O. Jacobs would like to thank Darryl Holm for providing him with a bridge from geometric mechanics to the field of image registration algorithms. Henry O. Jacobs is supported by the European Research Council Advanced Grant 267382 FCCA. Stefan Sommer is supported by the Danish Council for Independent Research with the project “Image-Based Quantification of Anatomical Change”.
Author Contributions
Henry O. Jacobs and Stefan Sommer wrote the paper.
Conflicts of Interest
The authors declare no conflict of interest.
References
- Younes, L. Shapes and Diffeomorphisms; Springer-Berlin: Heidelberg, Germany, 2010. [Google Scholar]
- Brox, T.; Bruhn, A.; Papenberg, N.; Weickert, J. High Accuracy Optical Flow Estimation Based on a Theory for Warping. In Computer Vision-ECCV 2004; aPajdla, T., Matas, J., Eds.; Springer-Berlin: Heidelberg, Germany, 2004; pp. 25–36. [Google Scholar]
- Dawn, S.; Saxena, V.; Sharma, B. Remote Sensing Image Registration Techniques: A Survey. In Image and Signal Processing; Elmoataz, A., Lezoray, O., Nouboud, F., Mammass, D., Meunier,, J., Eds.; Springer-Berlin: Heidelberg, Germany, 2010; pp. 103–112. [Google Scholar]
- Sotiras, A.; Davatzikos, C.; Paragios, N. Deformable Medical Image Registration: A Survey. IEEE Trans. Med. Imag 2013, 32, 1153–1190. [Google Scholar]
- Younes, L.; Arrate, F.; Miller, M.I. Evolutions equations in computational anatomy. NeuroImage 2009, 45, S40–S50. [Google Scholar]
- Boyes, R.G.; Rueckert, D.; Aljabar, P.; Whitwell, J.; Schott, J.M.; Hill, D.L.G.; Fox, N.C. Cerebral atrophy measurements using Jacobian integration: Comparison with the boundary shift integral. NeuroImage 2006, 32, 159–169. [Google Scholar]
- Gorbunova, V.; Jacobs, S.S.A.M.; Lo, P.; Dirksen, A.; Nielsen, M.; Bab-Hadiashar, A.; de Bruijne, M. Early detection of emphysema progression. In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2010; Springer-Berlin: Heidelberg, Germany, 2010; Volume 13, pp. 193–200. [Google Scholar]
- Joshi, S.C.; Miller, M.I.; Grenander, U. On the Geometry and Shape of Brain Sub-Manifolds. Int. J. Patt. Recogn. Artif. Intell 1997, 11, 1317–1343. [Google Scholar]
- Dupuis, P.; Grenander, U.; Miller, M.I. Variational Problems on Flows of Diffeomorphisms for Image Matching. Q. Appl. Math. 1998, 56, 587–600. [Google Scholar]
- Trouvé, A. An Infinite Dimensional Group Approach for Physics Based Models in Pattern Recognition, 1995. Available online: http://cis.jhu.edu/publications/papers_in_database/alain/trouve1995.pdf accessed on 28 April 2015.
- Christensen, G.; Rabbitt, R.; Miller, M. Deformable templates using large deformation kinematics. IEEE Trans. Image Process 1996, 5, 1435–1447. [Google Scholar]
- Holm, D.D.; Ratnanather, J.T.; Trouvé, A.; Younes, L. Soliton Dynamics in Computational Anatomy. NeuroImage 2004, 23, S170–S178. [Google Scholar]
- Bruveris, M.; Gay-Balmaz, F.; Holm, D.D.; Ratiu, T.S. The momentum map representation of images. J. Nonlinear Sci 2011, 21, 115–150. [Google Scholar]
- Jacobs, H. Symmetries in LDDMM with higher order momentum distributions, 2013, arXiv, p. 1306.3309. arXiv.org e-Print archive. Available online: http://arxiv.org/abs/1306.3309 accessed on 28 April 2015.
- Pennec, X.; Stefanescu, R.; Arsigny, V.; Fillard, P.; Ayache, N. Riemannian Elasticity: A Statistical Regularization Framework for Non-linear Registration. In MICCAI 2005; Springer-Berlin: Heidelberg, Germany, 2005; pp. 943–950. [Google Scholar]
- Derfoul, R.; Le Guyader, C. A. Relaxed Problem of Registration Based on the Saint Venant–Kirchhoff Material Stored Energy for the Mapping of Mouse Brain Gene Expression Data to a Neuroanatomical Mouse Atlas. SIAM J. Imag. Sci. 2014, 7, 2175–2195. [Google Scholar]
- Rueckert, D.; Sonoda, L.I.; Hayes, C.; Hill, D.L.; Leach, M.O.; Hawkes, D.J. Nonrigid registration using free-form deformations: Application to breast MR images. IEEE Trans. Med. Imag. 1999, 18, 712–721. [Google Scholar]
- Horn, B.K.; Schunck, B.G. Determining optical flow. Artif. Intell. 1981, 17, 185–203. [Google Scholar]
- Trouvé, A.; Younes, L. Local geometry of deformable templates. SIAM J. Math. Anal. 2005, 37, 17–59. [Google Scholar]
- Abraham, R.; Marsden, J.E. Foundations of Mechanics; Benjamin/Cummings Publishing Co. Inc.: Reading, MA, USA, 1978. [Google Scholar]
- Kobayashi, S.; Nomizu, K. Foundations of Differential Geometry, Volume I; Wiley Classics Library, John Wiley & Sons, Inc: New York, NY, USA, 1996. [Google Scholar]
- Taylor, M. Pseudo Differential Operators; Springer-Verlag: Berlin, Germany; New York, NY, USA, 1974. [Google Scholar]
- Holm, D.D.; Marsden, J.E. Momentum maps and measure-valued solutions (peakons, filaments, and sheets) for the EPDiff equation. In The Breadth of Symplectic and Poisson Geometry; Birkhäuser Boston: Boston, MA, USA, 2005; Volume 232, pp. 203–235. [Google Scholar]
- Micheli, M.; Glaunès, J.A. Matrix-valued Kernels for Shape Deformation Analysis. Geometry Imag. Comput 2014, 1, 57–39. [Google Scholar]
- Dupuis, P.; Grenander, U.; Miller, M.I. Variational problems on flows of diffeomorphisms for image matching. Quart. Appl. Math. 1998, 56, 587–600. [Google Scholar]
- Beg, M.F.; Miller, M.I.; Trouvé, A.; Younes, L. Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms. IJCV 2005, 61, 139–157. [Google Scholar]
- Joshi, S.; Miller, M. Landmark matching via large deformation diffeomorphisms. IEEE Trans. Image Process 2000, 9, 1357–1370. [Google Scholar]
- Vaillant, M.; Miller, M.; Younes, L.; Trouvé, A. Statistics on diffeomorphisms via tangent space representations. NeuroImage 2004, 23, S161–S169. [Google Scholar]
- Gay-Balmaz, F.; Vizman, C. Dual pairs in fluid dynamics. Ann. Glob. Anal. Geometry. 2012, 41, 1–24. [Google Scholar]
- Bauer, M.; Bruveris, M.; Michor, P.W. Overview of the Geometries of Shape Spaces and Diffeomorphism Groups. J. Math. Imag. Vis. 2014, 50, 60–97. [Google Scholar]
- Vaillant, M.; Glaunès, J. Surface Matching via Currents. In Information Processing in Medical Imaging; Springer-Berlin: Heidelberg, Germany, 2010; pp. 381–392. [Google Scholar]
- Glaunès, J. Transport Par Difféomorphismes de Points, de Mesures et de Courants Pour La Comparaison de Formes et L’anatomie Numérique. Ph.D. Thesis, Université Paris 13, Villetaneuse, France, 2005. [Google Scholar]
- Wells, W.M., III; Viola,, P.; Kikinis,, R. Multi-Modal Volume Registration by Maximization of Mutual Information. Med. Image Anal 1996, 1, 35–51. [Google Scholar]
- Roche, A.; Malandain, G.; Pennec, X.; Ayache, N. The Correlation Ratio as a New Similarity Measure for Multimodal Image Registration. In Proceedings of the First International Conference on Medical Image Computing and Computer-Assisted Intervention. In MICCAI 1998; Springer-Berlin: Heidelberg, Germany, 1998; pp. 1115–1124, ACM ID: 709612. [Google Scholar]
- Miller, M.I.; Trouvé, A.; Younes, L. Geodesic Shooting for Computational Anatomy. J. Math. Imag. Vis. 2006, 24, 209–228. [Google Scholar]
- Sommer, S.; Nielsen, M.; Darkner, S.; Pennec, X. Higher-Order Momentum Distributions and Locally Affine LDDMM Registration. SIAM J. Imag. Sci. 2013, 6, 341–367. [Google Scholar]
- Kolář, I.; Michor, P.W.; Slovák, J. Natural Operations in Differential Geometry; Springer Verlag: Berlin, Germany, 1999. [Google Scholar]
- Jacobs, H.O.; Sommer, S. Higher-order Spatial Accuracy in Diffeomorphic Image Registration. J. Geom. Imaging Comput, 2015. in press. Available online: http://arxiv.org/abs/1412.7504 accessed on 28 April 2015. [Google Scholar]
- Basser, P.J.; Mattiello, J.; LeBihan, D. Estimation of the effective self-diffusion tensor from the NMR spin echo. J. Magn. Resonance Ser. B 1994, 103, 247–254. [Google Scholar]
- Jian, B.; Vemuri, B.C.A. Unified Computational Framework for Deconvolution to Reconstruct Multiple Fibers From Diffusion Weighted MRI. IEEE Trans. Med. Imag 2007, 26, 1464–1471. [Google Scholar]
- Cao, Y.; Miller, M.I.; Winslow, R.L.; Younes, L. Large deformation diffeomorphic metric mapping of vector fields. IEEE Trans. Med. Imag. 2005, 24, 1216–1230. [Google Scholar]
- Alexander, D.C.; Gee, J.C.; Bajcsy, R. Strategies for Data Reorientation during Non-rigid Warps of Diffusion Tensor Images. In Medical Image Computing and Computer-Assisted Intervention—MICCAI’99; Taylor, C., Colchester, A., Eds.; Number 1679 in Lecture Notes in Computer Science; Springer-Verlag: Berlin, Germany, 1999; pp. 463–472. [Google Scholar]
- Alexander, D.C.; Pierpaoli, C.; Basser, P.J.; Gee, J.C. Spatial transformations of diffusion tensor magnetic resonance images. IEEE Trans. Med. Imag. 2001, 20, 1131–1139. [Google Scholar]
- Cheng, G.; Vemuri, B.C.; Carney, P.R.; Mareci, T.H. Non-rigid Registration of High Angular Resolution Diffusion Images Represented by Gaussian Mixture Fields. In Medical Image Computing and Computer-Assisted Intervention—MICCAI 2009; Springer-Verlag: Berlin, Germany, 2009; pp. 190–197. [Google Scholar]
- Arnold, V.I. Sur la géométrie différentielle des groupes de Lie de dimension infinie et ses applications à l’hydrodynamique des fluides parfaits. Ann. l’Instit. Fourier 1966, 16, 316–361. [Google Scholar]
- Mumford, D.; Michor, P.W. On Euler’s equation and “EPDiff”. J. Geometr. Mech 2013, 5, 319–344. [Google Scholar]
- Jacobs, H.O.; Ratiu, T.S.; Desbrun, M. On the coupling between an ideal fluid and immersed particles. Phys. D Nonlinear Phenome 2013, 265, 40–56. [Google Scholar]
- Cotter, C.J.; Holm, D.D.; Jacobs, H.O.; Meier, D.M. A jetlet hierarchy for ideal fluid dynamics. J. Phys. A Math. Theor. 2014, 47, 352001. [Google Scholar]
Notion | Notation |
---|---|
manifold | M |
tangent bundle | TM |
the group of diffeomorphisms | Diff(M) |
the space of vector-fields | |
the push-forward of u u(M) with respect to Φ ∈ Diff(M) | Φ*u |
the push-forward of m ∈ (M)* with respect to Φ 2∈ Diff(M) | Φ*m |
The Lie derivative (w.r.t. u ∈ (M)) | |
the dual space to (M) | (M)* |
Evaluation of m ∈ (M)* on u ∈ (M) | 〈m, u〉 |
© 2015 by the authors; licensee MDPI, Basel, Switzerland This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/4.0/).
Share and Cite
Sommer, S.; Jacobs, H.O. Reduction by Lie Group Symmetries in Diffeomorphic Image Registration and Deformation Modelling. Symmetry 2015, 7, 599-624. https://doi.org/10.3390/sym7020599
Sommer S, Jacobs HO. Reduction by Lie Group Symmetries in Diffeomorphic Image Registration and Deformation Modelling. Symmetry. 2015; 7(2):599-624. https://doi.org/10.3390/sym7020599
Chicago/Turabian StyleSommer, Stefan, and Henry O. Jacobs. 2015. "Reduction by Lie Group Symmetries in Diffeomorphic Image Registration and Deformation Modelling" Symmetry 7, no. 2: 599-624. https://doi.org/10.3390/sym7020599
APA StyleSommer, S., & Jacobs, H. O. (2015). Reduction by Lie Group Symmetries in Diffeomorphic Image Registration and Deformation Modelling. Symmetry, 7(2), 599-624. https://doi.org/10.3390/sym7020599