Fast Diffeomorphic Image Registration via Fourier-Approximated Lie Algebras

Abstract

This paper introduces Fourier-approximated Lie algebras for shooting (FLASH), a fast geodesic shooting algorithm for diffeomorphic image registration. We approximate the infinite-dimensional Lie algebra of smooth vector fields, i.e., the tangent space at the identity of the diffeomorphism group, with a low-dimensional, bandlimited space. We show that most of the computations for geodesic shooting can be carried out entirely in this low-dimensional space. Our algorithm results in dramatic savings in time and memory over traditional large-deformation diffeomorphic metric mapping algorithms, which require dense spatial discretizations of vector fields. To validate the effectiveness of FLASH, we run pairwise image registration on both 2D synthetic data and real 3D brain images and compare with the state-of-the-art geodesic shooting methods. Experimental results show that our algorithm dramatically reduces the computational cost and memory footprint of diffemorphic image registration with little or no loss of accuracy.

Appendix A

1.1 Properties of truncated convolution

We compute the truncated convolution in a bandlimited space by padding a sufficient number of zeros, $n_i-1$, at the end of each dimension, and then truncate the signal back to its original space. Note that (1) the convolution is modulo-$q_i$ ($q_i = 2n_i-1$) circular convolution, (2) the zero frequency component is shifted to the center of the domain.

We first introduce the kth element of a truncated convolution of any tangent vector field $\tilde{v}, \tilde{w} \in \tilde{V}$ as

$$\begin{aligned}{}[\tilde{v} *\tilde{w}]_k = \sum _{0 \leqslant l < q} \tilde{v}_{k-l}\tilde{w}_l, \end{aligned}$$

(18)

where $l=(l_1, \ldots , l_d)\in \mathbb {Z}^d$ are the padded frequency indices with $l_i \in \{0, \ldots , q_i\}$ along the ith axis. Let $M_i = \lfloor \frac{n_i-1}{2}\rfloor$, then $k=(k_1, \ldots , k_d) \in \mathbb {Z}^d$ denotes truncated frequency indexes with $k_i \in \{ M_{i}, \ldots , M_i + n_i - 1\}$. Notice here the $k-l$ represents $k-l\,(\mathrm{mod}\,q)$ with cyclic boundary conditions, and the vector field $\tilde{v}$ can also be replaced with a matrix field.

We then prove the commutativity and associativity of this truncated convolution. For any $\tilde{u}, \tilde{v}, \tilde{w} \in \tilde{V}$,

• Commutativity: $\tilde{u} *\tilde{v}= \tilde{v} *\tilde{u}$

Proof

The kth element of $\tilde{u} *\tilde{v}$ is

$$\begin{aligned}{}[\tilde{u} *\tilde{v}]_k = \sum _{0 \leqslant l < q} \tilde{u}_{k-l}\tilde{v}_l. \end{aligned}$$

By changing the coordinates using $k-l=k'$ where $k \in [M, M+n-1]$ and $k' \in [0, q]$ due to the cyclic condition, we rewrite the equation above as

$$\begin{aligned}{}[\tilde{u} *\tilde{v}]_k = \sum _{0 \leqslant k' < q} \tilde{u}_{k'}\tilde{v}_{k-k'} = [\tilde{v} *\tilde{u}]_k. \end{aligned}$$

$\square$

• Associativity: $(\tilde{u} *\tilde{v}) *\tilde{w} = \tilde{u} *(\tilde{v} *\tilde{w})$

Proof

The kth element of $(\tilde{u} *\tilde{v}) *\tilde{w}$ is

$$\begin{aligned}{}[(\tilde{u} *\tilde{v}) *\tilde{w}]_k =\sum _{0 \leqslant l'< q} \sum _{0 \leqslant l < q} \tilde{u}_{k-l'-l}\tilde{v}_l\tilde{w}_{l'}, \end{aligned}$$

where $k \in [M, M+n-1]$. We then rewrite the equation above by changing the coordinates using $k-l'-l=a$, where $a \in [0, q]$, as

$$\begin{aligned}{}[(\tilde{u} *\tilde{v}) *\tilde{w}]_k =\sum _{0 \leqslant l'< q} \sum _{0 \leqslant a < q} \tilde{u}_a\tilde{v}_{k-a-l'}\tilde{w}_{l'} = [\tilde{u} *(\tilde{v} *\tilde{w})]_k. \end{aligned}$$

$\square$

1.2 Deriving $\mathrm{ad}^*$ operator

Before showing the derivation of $\mathrm{ad}^*$ in a bandlimited space, we first introduce the pairing between a result of the truncated convolution (18) and a momentum vector field $\tilde{m} \in \tilde{V}^*$ is

$$\begin{aligned} (\tilde{m}, \tilde{v} *\tilde{w})&= \sum _{M \leqslant k< M+n}(\tilde{m}_{k-M}, \sum _{0 \leqslant l< q} \tilde{v}_{k-l}\tilde{w}_l) \\&= \sum _{M \leqslant k< M+n} \sum _{0 \leqslant l < q} \left( \tilde{m}_{k-M}, \tilde{v}_{k-l}\tilde{w}_l\right) \end{aligned}$$

We are now ready to derive the $\mathrm{ad}^*$ operator from its definition

$$\begin{aligned} (\mathrm{ad}_{\tilde{v}}^* \tilde{m}, \tilde{w}) = (\tilde{m}, \mathrm{ad}_{\tilde{v}} \tilde{w}). \end{aligned}$$

Since $\mathrm{ad}_{\tilde{v}} \tilde{w} = [\tilde{v}, \tilde{w}]$, after plugging in the definition of Lie bracket defined in (10) we have

$$\begin{aligned} (\mathrm{ad}_{\tilde{v}}^* \tilde{m}, \tilde{w})&= \left( \tilde{m}, (\tilde{D} \tilde{v}) *\tilde{w} - (\tilde{D} \tilde{w}) *\tilde{v}\right) \\&= \left( \tilde{m}, (\tilde{D} \tilde{v}) *\tilde{w}\right) - \left( \tilde{m}, (\tilde{D} \tilde{w}) *\tilde{v}\right) \\&= \sum _{M \leqslant k< M+n} \sum _{0 \leqslant l < q} (\tilde{m}_{k-M}, \tilde{v}_{k-l} \tilde{\eta }^T_{k-l} \tilde{w}_l)\\&\qquad - (\tilde{m}_{k-M}, \tilde{w}_l\tilde{\eta }^T_l \tilde{v}_{k-l}) \\ \end{aligned}$$

To separate $\tilde{w}$, we change coordinates by defining $k-M=k'$, $l-M=l'$ where $k' \in [0, n-1]$ and $l' \in [-M, q-M]$. For the purpose of notation simplicity, we drop the range of variables in the following equations as

$$\begin{aligned} (\mathrm{ad}_{\tilde{v}}^* \tilde{m}, \tilde{w})&= \sum _{k'} \sum _{l'} (\tilde{m}_{k'}, \tilde{v}_{k'-l'} \tilde{\eta }^T_{k'-l'} \tilde{w}_{l'+M}) \nonumber \\&\qquad - (\tilde{m}_{k'}, \tilde{w}_{l'+M}\tilde{\eta }^T_{l'+M} \tilde{v}_{k'-l'}) \nonumber \\&= \sum _{k'} \sum _{l'} \left( (\tilde{v}_{k'-l'} \tilde{\eta }^T_{k'-l'})^{T*} \tilde{m}_{k'}, \tilde{w}_{l'+M}\right) \nonumber \\&\qquad - (\tilde{v}^*_{k'-l'} \tilde{m}^T_{k'}, \tilde{w}_{l'+M}\tilde{\eta }^T_{l'+M} ) \nonumber \\&= \sum _{k'} \sum _{l'} \left( (\tilde{v}_{k'-l'} \tilde{\eta }^T_{k'-l'})^{T*} \tilde{m}_{k'}, \tilde{w}_{l'+M}\right) \nonumber \\&\qquad - (\tilde{v}^*_{k'-l'} \tilde{m}^T_{k'}\tilde{\eta }^*_{l'+M}, \tilde{w}_{l'+M}) \nonumber \\&= \sum _{k'} \sum _{l'} \left( (\tilde{v}_{k'-l'} \tilde{\eta }^T_{k'-l'})^{T*} \tilde{m}_{k'}, \tilde{w}_{l'+M}\right) \nonumber \\&\qquad + (\tilde{v}^*_{k'-l'} \tilde{m}^T_{k'}\tilde{\eta }_{l'+M}, \tilde{w}_{l'+M}) \nonumber \\&=\left( (\tilde{D} \tilde{v})^T \star \,\tilde{m}, \tilde{w}\right) + \left( \tilde{\Gamma } (\tilde{m} \otimes \tilde{v}), \tilde{w}\right) \nonumber \nonumber \\&=\left( (\tilde{D} \tilde{v})^{T} \star \,\tilde{m} + \tilde{\Gamma } (\tilde{m} \otimes \tilde{v}), \tilde{w}\right) . \end{aligned}$$

(19)

Finally, we have

$$\begin{aligned} \mathrm{ad}_{\tilde{v}}^* \tilde{m}&= (\tilde{D} \tilde{v})^{T} \star \, \tilde{m} + \tilde{\Gamma } (\tilde{m} \otimes \tilde{v}), \\ \mathrm{ad}_{\tilde{v}}^\dagger \tilde{v}&= \tilde{K} \mathrm{ad}_{\tilde{v}}^* \tilde{m} = \tilde{K} [(\tilde{D} \tilde{v})^{T} \star \,\tilde{m} + \tilde{\Gamma } (\tilde{m} \otimes \tilde{v})]. \end{aligned}$$