Lie Algebras for Fast Diffeomorphic Image Registration

. An image will be a square-integrable function on $\varOmega$ , that is, an element of $L^2(\varOmega , \mathbb {R})$ . Let $\mathrm{Diff}^\infty (\varOmega )$ denote the Lie group of smooth diffeomorphisms on $\varOmega$ . The Lie algebra for $\mathrm{Diff}^\infty (\varOmega )$ is the space $V = \mathfrak {X}^\infty (T\varOmega )$ of smooth vector fields on $\varOmega$ , which is the tangent space at the identity transform. We equip V with the weak Sobolev metric:

$\begin{aligned} \langle v, w \rangle _V = \int _{\varOmega } (Lv(x), w(x)) dx, \end{aligned}$

(1)

where $L = (I - \alpha \varDelta )^s$ is a symmetric, positive-definite, differential operator for some scalar, $$$\alpha > 0$$” src=”/wp-content/uploads/2016/09/A339424_1_En_19_Chapter_IEq10.gif”>, and integer power, s. The dual to the tangent vector v is a momentum, <IMG alt=$ . The notation (m, v) denotes the pairing of a momentum vector $m \in V^*$ with a tangent vector $v \in V$ . The metric L is an invertible operator, and $K = L^{-1}$ maps the momentum vector $m \in V^*$ back to the tangent vector $v = Km \in V$ .

LDDMM: In the LDDMM framework, we will consider diffeomorphisms that are generated by flows of time-varying velocity fields from V. More specifically, consider a time-varying velocity field, $v : [0,1] \rightarrow V$ , then we may define the flow $t \mapsto \phi _t \in \mathrm{Diff}^s(\varOmega )$ as a solution to the equation

$\frac{d \phi _t}{dt}(x) = v_t \circ \phi _t(x).$

The registration between two images, $I_0, I_1 \in L^2(\varOmega , \mathbb {R})$ , is the minimizer of the energy,

$\begin{aligned} E(v_t) = \int _0^1 \Vert v_t\Vert ^2_V \, dt + \frac{1}{2\sigma ^2} \Vert I_0 \circ \phi _1^{-1} - I_1\Vert ^2_{L^2}. \end{aligned}$

(2)

Geodesic Shooting: Given an initial velocity, $v_0 \in V$ , at

, the geodesic path $t \mapsto \phi _t \in \mathrm{Diff}^\infty (\varOmega )$ under the right-invariant Riemannian metric (1) is uniquely determined by the Euler-Poincaré equations (EPDiff) [1, 9],

$\begin{aligned} \frac{\partial v}{\partial t} = -\mathrm{ad}_{v}^\dagger v&= -K\mathrm{ad}_{v}^*m \nonumber \\&=-K\left[ (Dv)^Tm + Dm\, v + m\,{\text {div}}\,v\right] , \end{aligned}$

(3)

where D denotes the Jacobian matrix, and the operator $\mathrm{ad}^*$ is the dual of the negative Lie bracket of vector fields,

$\begin{aligned} \mathrm{ad}_v w = -[v,w] = Dv w - Dw v. \end{aligned}$

(4)

By integrating Eq. (3) forward in time, we generate a time-varying velocity $v_t : [0,1] \rightarrow V$ , which itself is subsequently integrated in time by the rule $(d\phi _t(x) /dt) = v_t \circ \phi _t(x)$ to arrive at the geodesic path, $\phi _t(x) \in \mathrm{Diff}^s(\varOmega )$ . Details are found in [14, 15]. Noting that the geodesic $\phi _t$ is fully determined by the initial condition

, we can rewrite the LDDMM image matching objective function (2) in terms of the initial velocity

$\begin{aligned} E(v_0) = ( Lv_0, v_0 )\, + \frac{1}{2\sigma ^2} \Vert I_0 \circ \phi _1^{-1} - I_1\Vert ^2. \end{aligned}$

(5)

Fig. 1.

Fourier coefficients of the discretized K operator on a $128 \times 128$ grid, with parameters $\alpha = 3$ , $$s = 3$$

3 Finite-Dimensional Lie Algebras of Diffeomorphism

While the Lie algebra $V = \mathfrak {X}^\infty (\varOmega )$ is an infinite-dimensional vector space, we must of course approximate smooth vector fields in V with finite-dimensional discretizations in order to represent them on a computer. In this section, we show that careful choices for the discretization of vector fields and their corresponding Lie brackets leads to a representation that is itself a finite-dimensional Lie algebra. The key observation is that the K operator is a low-pass filter, and as such, suppresses high frequency components in the velocity fields (see, for example, Fig. 1). As the K operator appears as the last operation on the right-hand side of the EPDiff Eq. (3), we can see that the velocity fields in the geodesic evolution do not develop high frequency components. This suggests that a standard implementation of geodesic shooting, using high-resolution velocity fields, wastes a lot of effort computing the high frequency components, which just end up being forced to zero by K. Instead, we propose to use a low-dimensional discretization of velocity fields as bandlimited signals in the Fourier domain.

More specifically, let $\tilde{V}$ denote the space of bandlimited velocity fields on $\varOmega$ , with frequency bounds $N_1, N_2, \ldots , N_d$ in each of the dimensions of $\varOmega$ . It will be convenient to represent an element $\tilde{v} \in \tilde{V}$ in the Fourier domain. That is, let $\tilde{v} \in \tilde{V}$ be a multidimensional array: $\tilde{v}_{k_1, k_2, \ldots , k_d} \in \mathbb {C}^d$ , where $k_i \in 0, \ldots , N_i - 1$ is the frequency index along the ith axis. Note that to ensure $\tilde{v}$ represents a real-valued vector field in the spatial domain, we have the constraint that $\tilde{v}_{k_1, \ldots , k_d} = \tilde{v}^*_{N_1 - k_1, \ldots , N_d - k_d}$ , where

denotes the complex conjugate. There is a natural inclusion mapping, $\iota : \tilde{V} \rightarrow V$ , of $\tilde{V}$ into the space V of smooth vector fields, given by the Fourier series expansion:

$\begin{aligned} \iota (\tilde{v})(x_1, \ldots , x_d) = \sum _{k_1 = 0}^{N_1} \cdots \sum _{k_d = 0}^{N_d} \tilde{v}_{k_1, \ldots , k_d} e^{2 \pi k_1 x_1} \cdots e^{2 \pi k_d x_d}. \end{aligned}$

(6)

Next, we define an operator that is a discrete analog of (4), the Lie bracket on continuous vector fields. First, we will denote by $\tilde{D} \tilde{v}$ the central difference Jacobian matrix of a discrete vector field, $\tilde{v} \in \tilde{V}$ . This can be computed in the discrete Fourier domain as a tensor product $\tilde{D} \tilde{v} = \eta \otimes \tilde{v}$ , where $\eta \in \tilde{V}$ is given by

$\eta _{k_1, k_2, \ldots , k_d} = (i \sin (2 \pi k_1), \ldots , i \sin (2 \pi k_d)).$

Second, we note that pointwise multiplication of matrix and vector field in the spatial domain corresponds to convolution in the Fourier domain. Because convolution between two bandlimited signals does not preserve the bandlimit, we must follow a convolution operation by truncation back to the bandlimits,

, in each dimension. We denote this truncated convolution operator between a matrix and a vector field as $\star$ . Now we are ready to define our discrete bracket operator for any two vectors $\tilde{v}, \tilde{w} \in \tilde{V}$ as

$\begin{aligned}{}[\tilde{v}, \, \tilde{w}] = (\tilde{D} \tilde{v}) \star \tilde{w} - (\tilde{D} \tilde{w}) \star \tilde{v}. \end{aligned}$

(7)

The next theorem proves that this operation satisfies the properties to be a Lie bracket on $\tilde{V}$ .

Theorem 1

The vector space $\tilde{V}$ , when equipped with the bracket operation (7), is a finite-dimensional Lie algebra. That is to say, $\forall \, \tilde{x}, \, \tilde{y}, \, \tilde{z} \in \tilde{V}$ and $a, \, b \in \mathbb {R}$ , the following properties are satisfied:

(a)

Linearity: $[a \tilde{x} + b \tilde{y}, \, \tilde{z}] = a \, [\tilde{x}, \, \tilde{z}] + b \, [\tilde{y}, \, \tilde{z}]$ ,

(b)

Anticommutativity: $[\tilde{x}, \, \tilde{y}]=-[\tilde{y}, \, \tilde{x}]$ ,

(c)

Jacobi identity: $[\tilde{x}, \, [\tilde{y}, \, \tilde{z}]] + [\tilde{z}, \, [\tilde{x}, \, \tilde{y}]] + [\tilde{y}, \, [\tilde{z}, \, \tilde{x}]] = 0$ .

Proof

Linearity and anticommutativity are immediate. We have

(a) Linearity:

$\begin{aligned}{}[a \tilde{x} + b \tilde{y}, \, \tilde{z}]&= \tilde{D} (a\tilde{x} + b\tilde{y})\star \tilde{z}- \tilde{D} \tilde{z} \star (a\tilde{x} + b\tilde{y}) \\&= a(\tilde{D} \tilde{x}\star \tilde{z}- \tilde{D} \tilde{z} \star \tilde{x}) + b(\tilde{D} \tilde{y} \star \tilde{z} - \tilde{D} \tilde{z} \star \tilde{y}) \\&= a [\tilde{x}, \, \tilde{z}] + b [\tilde{y}, \, \tilde{z}], \end{aligned}$