A method for large diffeomorphic registration via broken geodesics

Computational anatomy is an area of research focused on developing computational models of biological organs to study the anatomical variabilities in the deformation space. Anatomical variations arise due to structural differences across individuals and changes due to growth or atrophy in an individual. A common approach for quantifying the variability is to register the different images/shapes and parameterize the deformations with velocity fields. The registration algorithms typically optimize an energy functional based on a similarity function computed between the fixed and moving images.
Many initial image registration attempts use energy functionals inspired by physical processes to model the deformation as an elastic deformation [1], or viscous flow [2] or diffusion [3]. The diffusion-based approaches have been explored for 3D medical images in general [4] and with deformations constrained to be diffeomorphic [5] to ensure preservation of the topology. The two main approaches used to capture diffeomorphisms are parametric and nonparametric methods. The Free Form Deformation (FFD) model [6, 7] is a widely used parametric deformation model for medical image registration, where a rectangular grid with control points is used to model the deformation. Large diffeomorphic deformations [7] are handled by concatenating multiple FFDs. Deformable Registration via Attribute Matching and Mutual-Saliency Weighting (DRAMMS) [8] is a popular FFD-based method, which also handles inter-subject registration. DRAMMS matches Gabor features and prioritizes the reliable matching between images while performing registration. The main drawback of the deformations captured by FFD models is that they do not guarantee invertibility. The non-parametric methods represent the deformation with stationary or time varying velocity vector field. The diffeomorphic log-demon [5] is an example of the former while the Large Deformation Diffeomorphic Metric Mapping (LDDMM) [9] inspired from [10] is an example of the latter approach. In LDDMM, deformations are defined as geodesics on a Riemannian manifold, which is attractive; however, the methods based on this framework are computationally complex. The diffeomorphic log-demon framework [5], on the other hand, assigns a Lie group structure and assumes a stationary velocity field (SVF) which leads to computationally efficient methods, which is of interest to the community for practical purposes. This has motivated the exploration of a stationary LDDMM framework [11] that leverages the SVF advantage. The captured deformations are constrained to be symmetric in time-varying LDDMM [12] and log-demon [13] methods. The log-demon framework is of interest to the community for practical purposes because of its computational efficiency and simplicity.

The Lie group stucture gives a locally defined group exponential map to map the SVF to the deformation. Thus log-demon framework is meant to capture only neighboring elements in the manifold, i.e., only a limited degree of deformations can be captured. This will be referred to as the limited coverage issue of the SVF methods in this paper.

Notwithstanding the limited coverage, several SVF based methods have been reported for efficient medical image registration with different similarity metrics, sim, such as local correlation between the images [14], spectral features [15], modality independent neighborhood descriptors [16] and wavelet features [17, 18].

Lie group exponential map need not have any Riemannian metric associated with it as it is derived from the Lie group properties, and the log-demon framework is known to be metric free. The framework however, helps to analyze the underlying manifold and develop a representative group element called a template to study the variabilities among the manifold elements. The log-Euclidean mean [19] and the bi-invariant mean [20] are commonly used alternative approaches in diffeomorphic demons-based methods for template creation. Defining a metric space with SVF methods embedding complex deformations, is a challenging problem. A shape metric approximation method has been explored in [21] for the Lie group exponential map, with the shape metric taken to be the minimal length group exponential path connecting two shapes.

Defining a metric space with broader coverage for diffeomorphic deformations with Lie group structure is still unaddressed in the literature. The large deformations of interest are those that arise in the registration of inter-subject images/shapes. In this paper, we propose an image registration framework which can capture large diffeomorphisms. Any of the SVF based registration algorithms can be used in this framework. SVF based algorithms cannot handle complex deformations because the deformations are constrained to be smooth for the entire image and thus constrain the possible degree of deformation to some extent. We address this drawback by splitting the large deformation as compositions of smaller deformations.

Let $G$ be a finite-dimensional Lie group with Lie algebra $\mathfrak{g}$. Recall that $\mathfrak{g}$ is the tangent space $T_{e}G$ at the identity $e$ of $G$. The exponential map $\exp:\mathfrak{g}\rightarrow G$ is defined as follows: Let $v\in\mathfrak{g}$. Then $\exp(v)=\gamma_{v}(1)$, where $\gamma$ is the unique one-parameter subgroup of the Lie group $G$ with $v$ being its tangent vector at $e$. The vector $v$ is called the infinitesimal generator of $\gamma$. The exponential map is a diffeomorphism from a small neighborhood containing 0 in the Lie algebra $\mathfrak{g}$ to a small neighborhood containing $e$ of $G$.

Due to the fact that a bi-invariant metric may not exist for most of the Lie groups considered in medical image registration, the deformations considered here are elements of a Lie group with the Cartan-Schouten Connection [22]. This is the same as the one considered in the log-demon framework [5]. This is a left invariant connection [19] in which geodesics through the identity are one-parameter subgroups. The group geodesics are the geodesics of the connection. Any two neighboring points can be connected with a group geodesic. That is, if the stationary velocity field $v$ connecting two images in the manifold $\mathcal{G}$ is small enough, then its group exponential map forms a geodesic. Similarly every $\mathfrak{g}\in G$ has a geodesically convex open neighbourhood [19].

SVF based registration methods capture only a limited degree of deformation because exponential mappings are only locally defined. In order to perform registration of a moving image towards a fixed image, SVF is computed iteratively by updating it with a smoothed velocity field. This update is computed via a similarity metric that measures the correspondence between the moving and fixed images. The spatial smoothing has a detrimental effect as we explain next. A complex deformation typically consists of spatially independent deformations in a local neighbourhood. Depending on the smoothing parameter value, only major SVF updates in each region is considered for registration. Thus, modeling complex deformations with a smooth stationary velocity field is highly dependent on the similarity metric and the smoothing parameter in a registration algorithm. Finding an ideal similarity metric and an appropriate smoothing parameter applicable for any registration problem, irrespective of the complexity of the deformation and the type of data, is difficult. We propose to address this issue as follows: Deform the moving image toward the fixed image by sequentially applying an SVF based registration. The SVF based algorithm chooses the major or the predominant (correspondence-based) deformation component among the spatially independent deformations in all the neighbourhoods to register along these predominant directions. The subsequent steps in the algorithm captures the next set of predominant directions sequentially. These sequentially captured deformations has a decreasing order of degree of pixel displacement caused by the deformations. Mathematically speaking, the discussion above can be summarised as follows. Consider complex deformations as a composition of finite group geodesics and use a registration metric approximation to quantify the deformation between two images in terms of the length of a broken geodesic connecting them; a broken geodesic is a piecewise smooth curve, where each curve segment is a geodesic.

In the proposed method, the similarity-based metric selects the predominant here deformation in each sequential step. The deformation that can bring the moving image in a step maximally closer to the target is selected from the one-parameter subgroup of deformations. In the manifold $\mathcal{G}$ every geodesic is contained in a unique maximal geodesic. Hence the maximal group geodesic $\gamma_{i}$ computed using log-demon registration framework deforms the sequential image $S_{i-1}$ in the previous step maximally closer to $S_{N}$. The maximal group geodesic paths are composed to get the broken geodesic path. As the deformation segments are diffeomorphic, the composed large deformation of the segments is also diffeomorphic.

In the proposed method, the coverage of the SVF method and the degree of deformation determines the number of subgroups $N$ needed to cover the space. The feature based SVF methods in general, give more coverage for a single such subgroup and reduce the value of $N$. The local deformations on the manifold $\mathcal{G}$ form a one-parameter subgroup, which is a finite-dimensional path connected subgroup. The topological space formed by composing the group geodesics also forms a path connected manifold, as shown in Figure 1. So a broken geodesic is a good choice to capture complex deformations as it can be computed between any two points in the manifold along with the metric associated with the path.

A broken geodesic $\gamma:[0,T]\rightarrow M$ has finite number of geodesic segments $\gamma_{i}$ for partitions of the domain $0<t_{1}<t_{2}<\cdots<t_{i}<\cdots t_{N}=T$ where $i=1,....N$. The proposed algorithm to deform $S_{0}$ towards $S_{N}$ is given in Algorithm 1. We have chosen the registration algorithm from [13] to compute SVF, $u_{i}$, in Algorithm 1. The Energy term is as defined below where the first term is a functional of the similarity measure, which captures the correspondence between images, with $\mbox{sim}(S_{N},S_{i})=S_{N}-S_{i-1}\circ\exp(v_{i})$. The second term is a regularization term, with Reg$(\gamma_{i})=\lVert\bigtriangledown\gamma_{i}\rVert^{2}$.

The proposed method was implemented using a simple intensity based diffeomorphic log-demon technique [23] for illustrating the concept which is openly available at: http://dx.doi.org/10.17632/29ssbs4tzf.1. This choice also facilitates understanding the key strengths of the method independently. Two state-of-the-art (SOTA) methods are considered for performance comparison with the proposed method: the symmetric LDDMM implementation in ANTs [12] and DRAMMS which is a feature based, free-form deformation estimation method [8]. These two methods are considered to be good tools for inter-subject registration [24]. Publicly available codes were used for the SOTA methods with parameter settings as suggested in [24] for optimal performance. Both methods were implemented with B-spline interpolation, unless specified. 3D registration was done and the images used in the experiments were sourced from [25] and [26] unless specified.

The proposed method and registration metric were validated with both controlled and natural deformations.

Group exponential map based methods, with simple similarity registration metrics, fail to capture large deformations as the map is local in nature. We have addressed this issue in this paper by modelling large deformations with broken geodesic paths with the path length taken to be the associated registration metric. The proposed framework/method makes the underlying manifold path-connected. The results of implementation with a simple log-demon method show the performance to be superior to SOTA methods for complex/large deformations. Efficient implementation of the proposed approach are currently being explored.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.