Information-Theoretic Registration with Explicit Reorientation of Diffusion-Weighted Images

TL;DR

This paper introduces LORD, an information-theoretic registration method for diffusion-weighted images that explicitly optimizes over directional scales, improving alignment of complex fiber structures.

Contribution

It extends hierarchical scale-space models to include directional information and non-rigid deformations, enhancing DWI registration accuracy.

Findings

Successfully deforms orientation distribution functions (ODFs)

Handles complex fiber configurations like crossings and fanning

Improves deformation retrieval over scalar-based methods

Abstract

We present an information-theoretic approach to the registration of images with directional information, and especially for diffusion-Weighted Images (DWI), with explicit optimization over the directional scale. We call it Locally Orderless Registration with Directions (LORD). We focus on normalized mutual information as a robust information-theoretic similarity measure for DWI. The framework is an extension of the LOR-DWI density-based hierarchical scale-space model that varies and optimizes the integration, spatial, directional, and intensity scales. As affine transformations are insufficient for inter-subject registration, we extend the model to non-rigid deformations. We illustrate that the proposed model deforms orientation distribution functions (ODFs) correctly and is capable of handling the classic complex challenges in DWI-registrations, such as the registration of fiber-crossings along with kissing, fanning, and interleaving fibers. Our experimental results clearly illustrate a novel promising regularizing effect, that comes from the nonlinear orientation-based cost function. We show the properties of the different image scales and, we show that including orientational information in our model makes the model better at retrieving deformations in contrast to standard scalar-based registration.

Equations

$$h_{\beta \alpha \sigma }(i,j|\phi ,\mathbf{x})=$$

$${\int }_{{\mathbb{R}}^3}{P}_{\beta }({\mathbf{I}}_{\sigma }(\phi (\mathbf{y}))-i)P_{\beta }({\mathbf{J}}_{\sigma }(\mathbf{y})-j)W_{\alpha }(\mathbf{y}-\mathbf{x})d\mathbf{y}$$

$$p_{\beta \alpha \sigma }(i,j|\phi ,\mathbf{x})\simeq \frac{h_{\beta \alpha \sigma }(i,j|\phi ,\mathbf{x})}{{\int }_{{\Lambda }^2}{h}_{\beta \alpha \sigma }(k,l|\phi ,\mathbf{x})dkdl}$$

$${\mathcal{F}}_{lin}={\int }_{{\Lambda }^2}f(i,j)p(i,j)didj$$

$${\mathcal{F}}_{non-lin}={\int }_{{\Lambda }^2}f(p(i,j))didj.$$

$${\mathbf{I}}_{\sigma \kappa }(\mathbf{x},\mathbf{v})$$

$${\Gamma }_{\kappa }(\mathbf{\nu },\mathbf{v})=C{e}^{\kappa ({\mathbf{\nu }}^{\top }\mathbf{v}{)}^2}$$

$$C=M\left(\frac{1}{2},\frac{1}{d},\kappa \right)$$

$$\Phi :{\mathbb{R}}^3\times {\text{SS}}^2\to {\mathbb{R}}^3\times {\text{SS}}^2,\Phi (\mathbf{x},\mathbf{v})=\left(\phi (\mathbf{x}),{\psi }_{\mathbf{x}}(\mathbf{v})\right).$$

$$\begin{gathered} h_{\mathbf{I}\circ \Phi ,\mathbf{J};\sigma \kappa \beta \alpha }(i,j|\mathbf{x})= \\ {\int }_{{\mathbb{R}}^3\times {\text{SS}}^2}{P}_{\beta }({\mathbf{I}}_{\sigma \kappa }(\phi (\mathbf{y}),\psi (\mathbf{v}))-i)P_{\beta }({\mathbf{J}}_{\sigma \kappa }(\mathbf{y},\mathbf{v})-j)\times \\ {W}_{\alpha }(\mathbf{y}-\mathbf{x})d\mathbf{y}d\mathbf{v} \end{gathered}$$

$$p_{\mathbf{I}\circ \Phi ,\mathbf{J};\sigma \kappa \beta \alpha }(i,j|\mathbf{x})=\frac{h_{\mathbf{I}\circ \Phi ,\mathbf{J};\sigma \kappa \beta \alpha }(i,j|\mathbf{x})}{{\int }_{{\Lambda }^2}{h}_{\mathbf{I}\circ \Phi ,\mathbf{J};\sigma \kappa \beta \alpha }(l,k|\mathbf{x})dldk}.$$

$$\begin{gathered} h_{\mathbf{I}\circ \Phi ,\mathbf{J};\sigma \kappa \beta }(i,j)= \\ {\int }_{{\mathbb{R}}^3\times {\text{SS}}^2}{P}_{\beta }({\mathbf{I}}_{\sigma \kappa }(\phi (\mathbf{x}),\psi (\mathbf{v}))-i)P_{\beta }({\mathbf{J}}_{\sigma \kappa }(\mathbf{x},\mathbf{v})-j)d\mathbf{x}d\mathbf{v} \end{gathered}$$

$$p_{\mathbf{I}\circ \Phi ,\mathbf{J};\sigma \kappa \beta }(i,j)=\frac{h_{\mathbf{I}\circ \Phi ,\mathbf{J};\sigma \kappa \beta }(i,j)}{{\int }_{{\Lambda }^2}{h}_{\mathbf{I}\circ \Phi ,\mathbf{J};\sigma \kappa \beta }(l,k)dldk}.\text{hfil}$$

$$B_{{\mathbf{c}}^r}^r(\mathbf{x})=\sum _{l,m,n=0}^3{B}_l(u)B_m(v)B_n(w){\mathbf{c}}_{i+l,j+m,k+l}^r$$

$${\phi }_{\mathbf{c}}(\mathbf{x})=\mathbf{x}+\sum _{r=1}^R{B}_{{\mathbf{c}}^r}^r(\mathbf{x}).$$

$${\psi }_{\mathbf{c}}(\mathbf{v})=\frac{\mathbf{v}+\left(J_{\mathbf{x}}\mathbf{B}\cdot \mathbf{v}\right)\mathbf{c}}{|\mathbf{v}+\left(J_{\mathbf{x}}\mathbf{B}\cdot \mathbf{v}\right)\mathbf{c}|}$$

$$\mathcal{F}(\mathbf{c})=\text{NMI}(\mathbf{I}\circ {\Phi }_{\mathbf{c}},\mathbf{J};\sigma \kappa \beta )=\frac{H_{\sigma \kappa \beta }(\mathbf{I}\circ {\Phi }_{\mathbf{c}})+H_{\sigma \kappa \beta }(\mathbf{J})}{H_{\sigma \beta \kappa }(\mathbf{I}\circ {\Phi }_{\mathbf{c}},\mathbf{J})}.$$

$${\mathcal{S}}^r(\mathbf{c})=-\frac{{\lambda }_r}{2}\sum _{i=1}^{p_r}\Vert c_i^r-\frac{1}{|N^r(i)|}\sum _{j\in N^r(i)}{c}_j^r{\Vert }^2$$

$${\mathbf{c}}^{\ast }={\mathrm{argmax}}_{\mathbf{c}}\mathcal{M}(\mathbf{c})={\mathrm{argmax}}_{\mathbf{c}}\left(\mathcal{F}(\mathbf{c})+\mathcal{S}(\mathbf{c})\right).$$

$$J_{\mathbf{c}}NMI(\mathbf{I}\circ {\Phi }_{\mathbf{c}},J)=d_{\mathbf{I}}\mathcal{L}({\mathbf{I}}_{\sigma \kappa }\circ {\Phi }_{\mathbf{c}})J_{\mathbf{c}}\left({\mathbf{I}}_{\sigma \kappa }\circ {\Phi }_{\mathbf{c}}\right)$$

$$NMI(\mathbf{I}\circ {\Phi }_{\mathbf{c}},\mathbf{J};\sigma \kappa \beta )$$

$${\mathbf{I}}_{\sigma \kappa }$$

$$\sum _{n=1}^N{\int }_{{\mathbb{R}}^3}\mathbf{I}(\mathbf{y},{\nu }_n)K_{\sigma }({\phi }_{\mathbf{c}}(\mathbf{x})-\mathbf{y}){\bar{\Gamma }}_{\kappa }({\nu }_n,{\psi }_{\mathbf{c}}(\mathbf{v}))d\mathbf{y}$$

$${\bar{\Gamma }}_{\kappa }({\nu }_n,{\psi }_{\mathbf{c}}(\mathbf{v}))=\frac{e^{\kappa ({\nu }_n^{\top }{\psi }_{\mathbf{c}}(\mathbf{v}){)}^2}}{{\sum }_{i=1}^N{e}^{\kappa ({\nu }_i^{\top }{\psi }_{\mathbf{c}}(\mathbf{v}){)}^2}}.$$

$$\begin{gathered} J_{\mathbf{c}}{\mathbf{I}}_{\sigma \kappa }({\phi }_{\mathbf{c}}(\mathbf{x}),{\psi }_{\mathbf{c}}(\mathbf{v}))= \\ \sum _{n=1}^N\left({\int }_{{\mathbb{R}}^3}\mathbf{I}(\mathbf{y},{\nu }_n)J_{\mathbf{c}}{K}_{\sigma }({\phi }_{\mathbf{c}}(\mathbf{x})-\mathbf{y})d\mathbf{y}\right){\bar{\Gamma }}_{\kappa }({\nu }_n,{\psi }_{\mathbf{c}}(\mathbf{v}))+ \\ \sum _{n=1}^N\left({\int }_{{\mathbb{R}}^3}\mathbf{I}(\mathbf{y},{\nu }_n)K_{\sigma }({\phi }_{\mathbf{c}}(\mathbf{x})-\mathbf{y})d\mathbf{y}\right)J_{\mathbf{c}}{\bar{\Gamma }}_{\kappa }({\nu }_n,{\psi }_{\mathbf{c}}(\mathbf{v})) \end{gathered}$$

$$J_{\mathbf{c}}{K}_{\sigma }({\phi }_{\mathbf{c}}(\mathbf{x})-\mathbf{y})=-\frac{K_{\sigma }({\mathbf{\phi }}_{\mathbf{c}}(\mathbf{x})-\mathbf{y})}{{\sigma }^2}{\left({\phi }_{\mathbf{c}}(\mathbf{x})-\mathbf{y}\right)}^{\top }\mathbf{B}(\mathbf{x})$$

$$J_{\mathbf{c}}f({\nu }^{\top }{\psi }_{\mathbf{c}}(\mathbf{v}))=2\kappa f({\nu }^{\top }{\psi }_{\mathbf{c}}(\mathbf{v})){\psi }_{\mathbf{c}}(\mathbf{v}{)}^{\top }\nu {\nu }^{\top }{J}_{\mathbf{c}}{\psi }_{\mathbf{c}}(\mathbf{v}).$$

$$J_{\mathbf{c}}{\bar{\Gamma }}_{\kappa }({\nu }_n,{\psi }_{\mathbf{c}}(\mathbf{v}))=$$

$$-2\kappa \frac{{\bar{\Gamma }}_{\kappa }({\nu }_n,{\psi }_{\mathbf{c}}(\mathbf{v}))}{{\sum }_{i=1}^N{f}_i}{\psi }_{\mathbf{c}}(\mathbf{v}{)}^{\top }\left(\sum _{i\ne n}{f}_i{\nu }_i{\nu }_i^{\top }\right)J_{\mathbf{c}}{\psi }_{\mathbf{c}}(\mathbf{v}).$$

$$J_{\mathbf{c}}{\psi }_{\mathbf{c}}(\mathbf{v})=\frac{1}{|\mathbf{V}|}{\pi }_{{\mathbf{V}}^{\perp }}{J}_{\mathbf{x}}\mathbf{B}\cdot \mathbf{v}$$

$$\begin{gathered} J_{\mathbf{c}}{\mathbf{I}}_{\sigma \kappa }({\phi }_{\mathbf{c}}(\mathbf{x}),{\psi }_{\mathbf{c}}(\mathbf{v}))= \\ -\sum _{n=1}^N{\int }_{{\mathbb{R}}^3}\mathbf{I}(\mathbf{y},{\nu }_n)K_{\sigma }({\phi }_c(\mathbf{x})-\mathbf{y}){\bar{\Gamma }}_{\kappa }({\nu }_n,{\psi }_{\mathbf{c}}(\mathbf{v}))\times \\ \{\frac{({\phi }_{\mathbf{c}}(\mathbf{x})-\mathbf{y}{)}^{\top }\mathbf{B}(\mathbf{x})}{{\sigma }^2}+ \\ \frac{2\kappa {\mathbf{V}}^{\top }\left({\sum }_{i\ne n}{f}_i{\nu }_i{\nu }_i^{\top }\right)}{|\mathbf{V}|{\sum }_{i=1}^N{f}_i}\mathbf{\pi }{}_{\mathbf{V}{}^{\perp }}\left(\frac{J_{\mathbf{x}}\mathbf{B}\bullet \mathbf{v}}{|\mathbf{V}|}\right)\}d\mathbf{y}. \end{gathered}$$

$$\nabla \mathcal{F}(\mathbf{c})={\left(d_{\mathbf{I}}\mathcal{L}({\mathbf{I}}_{\sigma \kappa }\circ {\Phi }_{\mathbf{c}})J_{\mathbf{c}}{\mathbf{I}}_{\sigma \kappa }\circ {\Phi }_{\mathbf{c}}\right)}^T.$$

$${\mathcal{S}}^r({\mathbf{c}}^r)=-\frac{{\lambda }_r}{2}\sum _{i=1}^{p_r}\frac{1}{|N^r(i){|}^2}\Vert T_i^r{\mathbf{c}}^r{\Vert }^2$$

Full text

1 Introduction

In this work, we study registration problems for images of the type $I:\Omega\to{\mathcal{F}}({\mathbb{P}}^{2})$ with $\Omega$ an open, bounded domain of ${\mathbb{R}}^{3}$, the image domain, ${\mathbb{P}}^{2}$ is the projective plane of directions (vector lines) in ${\mathbb{R}}^{3}$ and ${\mathcal{F}}({\mathbb{P}}^{2})$ is a suitable function space on ${\mathbb{P}}^{2}$: each $I(\bm{x})$ is a function $I(\bm{x}):{\mathbb{P}}^{2}\to{\mathbb{R}}$ which provides a measure of a phenomenon along a three-dimensional direction. We may and will also reinterpret such an image $I$ as an image ${\bf I}:\Omega\times{\mathbb{P}}^{2}\to{\mathbb{R}}$, ${\bf I}(\bm{x},\bm{v})=I(\bm{x})(\bm{v})$. These functions can be used to represent orientation distribution functions used in Diffusion-Weighted Imaging (DWI).

DWI is a non-invasive Magnetic Resonance Imaging (MRI) protocol that can be used to infer microstructures of biological tissues by tracking the movement of water molecules, otherwise invisible in structural MRI. However, the spatio-directional geometry of DWI signals makes it a challenge in image registration, a key tool for comparing and analyzing medical image data. Also, DWI acquired on different scanners or with different protocols have a non-linear relationship, and defining the similarity between two DWIs is inherently difficult johansen2013diffusion . Information-theoretic similarity measures such as mutual information (MI) and normalized mutual information (NMI) have been shown to handle both linear and non-linear relationships very well darkner2013locally ; studholme1999overlap ; darkner2018collocation and are therefore attractive for use in the registration of DWI.

Our contribution is a scale-space formulation of density estimation for image similarity for DWI in a nonrigid registration setting. The model encompasses explicit reorientation providing the computational framework for estimating non-linear similarity measures well-suited for DWI. The model allows for image similarities such as NMI by extending the Locally Orderless Registration for Diffusion-Weighted Images (LOR-DWI) jensen2015locally to nonrigid registration as affine transformations provide insufficient flexibility to perform an accurate inter-subject registration. LORD includes directional information in the image similarity, which allows explicit optimization over the reorientation of diffusion gradients. This extends from diffusion tensor images (DTI) to the raw high-angular resolution diffusion images (HARDI) or the topographically inverted Orientation DensityFunctions (ODF).

To demonstrate the capabilities of the framework, we first build an objective function over a spline base parametric deformation field, from NMI and a simple regularization, in a multi-resolution setting. We partially compute the gradient of the objective – as some of the computations are standard – and indicate how it is optimized. Then we experimentally study some of its properties, such as the effect of the reorientation on the ODFs, the effect of the scale-space on optimization and regularization on simulated DWI data and on a synthetic deformation of data from a subject of the Human Connectome Project (HCP) van2013wu . We show that LORD preserves the ODFs and produces excellent mappings for crossing, kissing and curving fibers, while empirically providing an inherent regularizing effect. A example PyTorch implementation is provided in https://github.com/sunedarkner/Registration.

The density formulation also allows us to optimize over the isocurves of the DWI signal $I:\mathbb{R}^{3}\times\mathbb{P}^{2}\rightarrow\mathbb{R}$. In principle, this framework could be extended to other non-trivial geometry of other image modalities where scale-space structures can be defined. In this work, we focus on DWI data where we use a classical scale-space structure on this non-linear geometry.

2 Related Work

This work addresses two of the major challenges in the voxel-based registration of DWI: the reorientation of DWI in image registration, and the non-linear similarity between DWI.

Image registration refers to a process that transforms data into a shared coordinate system. For DWI, the common way to register two images is to use scalar-based methods on quantitative measures, such as the fractional anisotropy (FA) or the mean diffusivity (MD) o2017advances . However, as such methods disregard most of the directional information in DWI, methods have been developed to also account for the reorientation of the diffusion profile. Most of these methods are created on top of scalar-based methods and iteratively reorients the gradients after changing the deformation field. Some of the most popular approaches can be found in avants2009advanced ; tournier2012mrtrix ; jenkinson2012fsl ; fischl2012freesurfer . However, registration frameworks have also been designed with an objective function that explicitly optimizes over the reorientation of the gradients, such as DTI-TK zhang2006deformable , DT-REFineD yeo2009dt , and the more recent DR-TAMAS irfanoglu2016dr . These frameworks have been shown to generally outperform scalar-based frameworks for DWI wang2011dti ; zhang2014large ; wang2016evaluation ; wang2017evaluations .

For scalar based approaches to image similarity, registration is straightforward as any popular scalar-based measure can be used, e.g. the sum-of-squared difference (SSD) yeo2009dt or mutual information (MI) wells1996mmv ; wang2017evaluations . Explicit reorientation strategies inherently define the similarity over both position and orientation, and once the full diffusion profile is part of the similarity, these similarity-measures can be formulated in a well-suited way for the non-linear relationship between pairs of DWI. Both zhang2006deformable and yeo2009dt used variations of SSD in the objective function, while MI was used in jensen2015locally through an extension of the LOR framework in darkner2013locally ; darknersporring2011ipmi . As argued in jensen2015locally , the invariant and statistical properties of MI makes it a logical choice for DWI, where multiple factors result in a less functional, more statistical relationship, such as variations in $b$-values, non-monoexponential behavior in biological tissues, and inter-scanner variability johansen2013diffusion . MI is often used in the standard registration of complex modalities and, as such, in scalar-based registration of DWI van2007nonrigid , and in the pre-processing of DWI, treiber2016characterization . MI and normalized MI (NMI) studholme1999overlap provide a non-linear statistical measure for DWI. It is nevertheless possible that more functional measures can be well-suited, among others, cross-correlation (CC) and normalized CC (NCC), both of which can be defined from the LOR density-formulation sporring2011jacobians ; darknersporring2012pami . The density-based DWI comes from the generalized way of estimating image similarity measures based on Locally Orderless Images (LOI) koenderink1999structure . The first mention of LOI in the context of image registration was in hermosillo2002variational where a variational approach to image registration was presented. The LOR framework darknersporring2011ipmi ; darknersporring2012pami generalized a range of similarity measures as linear and non-linear functions of density estimates for scalar-valued images.

3 Locally Orderless Directional Imaging

3.1 Notations

$\Omega\subset\mathbb{R}^{3}$ is the spatial domain of the images under consideration. A scalar image is a function $\bm{I}:\Omega\to\mathbb{R}^{3}$. We assume that we can extend it on the whole $\mathbb{R}^{3}$, for instance by extending it with 0. The convolution of two functions $f,g:{\mathbb{R}}^{n}\to{\mathbb{R}}$ is defined as $f*g(\bm{x})=\int_{{\mathbb{R}}^{n}}f(\bm{y})g(\bm{x}-\bm{y})\,d\bm{y}$. The projective space of directions of $\mathbb{R}^{3}$, i..e, the set of vector lines of ${\mathbb{R}}^{3}$, is denoted by $\mathbb{P}^{2}$, and the unit sphere of $\mathbb{R}^{3}$ by $\mathbb{S}^{2}$. We will encounter spatio-directional images $\bm{I}:\Omega\times\mathbb{P}^{2}\to\mathbb{R}$, which we similarly assume to be extendable to $\mathbb{R}^{3}\times\mathbb{P}^{2}$. This is necessary in both cases so as to define their spatial smoothing via convolution. We use the following elementary property: As $\mathbb{P}^{2}$ is naturally identified as the quotient $\mathbb{S}^{2}/\{\pm 1\}$ by the antipodal symmetry, a function $f:\mathbb{P}^{2}\to\mathbb{R}$ can be lifted to an antipodal symmetric function $\tilde{f}:\mathbb{S}^{2}\to\mathbb{R}$. Conversely, any antipodal symmetric function $g:\mathbb{S}^{2}\to\mathbb{R}$ factors through $\mathbb{P}^{2}$. A spatio-directional image can and will be lifted to an antipodal symmetric image $\tilde{\bm{I}}:\Omega\times\mathbb{S}^{2}$: $\tilde{\bm{I}}(\bm{x},-\bm{v})=\tilde{\bm{I}}(\bm{x},\bm{v})$. Such a lift allows us to represent ${\mathbb{P}}^{2}$-functions as $\SS^{2}$ ones, with the benefit that $\SS^{2}$ is an oriented manifold, while ${\mathbb{P}}^{2}$ is not. This means in particular that ${\mathbb{P}}^{2}$ does not have a global volume form. On the other hand, $\SS^{2}$ has, and we use the one induced by its standard Riemannian metric.

We will denote by $\bm{I}$ both the spatio-directional image and its antipodal symmetric lifting in the following.

3.2 Recalling the LOR Framework

The LOR framework defines the density estimates over three scales: The image scale $\sigma$, the intensity scale $\beta$, and the integration scale $\alpha$. In the context of scalar registration, for a transformation, $\varphi:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$, the estimated histogram $h$ and the corresponding density $p$ is computed as

[TABLE]

where $i,j\in[a_{1},a_{2}]=:\Lambda$ are values in the image intensity range, $\bm{I}_{\sigma}(\bm{\varphi}(\bm{x}))=(\bm{I}*K_{\sigma})(\bm{\varphi}(\bm{x}))$ and $J_{\sigma}(\bm{x})=(\bm{J}*K_{\sigma})(\bm{x})$ are images convolved with the kernel $K_{\sigma}$ with standard deviation $\sigma$, $P_{\beta}$ is a Parzen-window of scale $\beta$ (bishop:2006:PRML , paragraph 2.5.1), and $W_{\alpha}$ is a Gaussian integration window of scale $\alpha$. The marginals are trivially obtained by integration over the appropriate variable. The LOR-approach to similarity lets us use a set of generalized similarity measures, linear or non-linear

[TABLE]

The linear measure $f(i,j)$ includes e.g. sum of squared differences and Huber. The non-linear $f(p(i,j))$ includes e.g. MI, normalized MI (NMI), see darknersporring2012pami for details.

3.3 The LORD framework

This work addresses the estimation of the image similarity $\mathcal{F}$ of DWI in the context of nonrigid registration an extension of our previous work jensen2015locally . In this context, ${\bf I}$ and ${\bf J}$ are spatio-directional signals. Specifically, DWI MR attenuation signals at location $\bm{x}$, for a gradient direction $\bm{v}$, are modeled by $S(\bm{x},\bm{v})=S_{0}(\bm{x})e^{-b\bm{I}(\bm{x},\bm{v})}$  tao2006method and apparent diffusion coefficients volumes are given by $\bm{I}(\bm{x},\bm{v})=-\frac{1}{b}\log\frac{S(\bm{x},\bm{v})}{S_{0}(\bm{x})}$. Gradient directions $\bm{v}$ belong to $\mathbb{S}^{2}$ but the diffusion is orientation-free, $\bm{I}(\bm{x},\bm{v})\approx\bm{I}(\bm{x},-\bm{v})$. This defines naturally a spatio-directional image $\Omega\times{\mathbb{P}}^{2}\to{\mathbb{R}}$. In order to apply LORD, the histogram and density estimates Section 3.2 and Equation 3.2 must be extended to spatio-directional data, and the action of the spatial transformation on the directions must be defined.

We introduce a kernel on the sphere as an extension to the density estimates of LOR to include directional information. This kernel accounts for directional smoothing and defines the LORD framework: we extend spatial smoothing to be spatio-directional. The directional smoothing preserves the antipodal symmetry, equivalently the projective structure, via a symmetric kernel $\Gamma_{\kappa}(\bm{\nu},\bm{v})$ on $\mathbb{S}^{2}$. We define the smoothed signal $\bm{I}_{\sigma\kappa}$ at scales $(\sigma,\kappa)$ by

[TABLE]

Here $K_{\sigma}$ is a Gaussian kernel with standard deviation $\sigma$. The integral over $\SS^{2}$ is sometimes referred to as a spherical convolution (though actually a pseudo-convolution GallierQuaintance:2020b ). $\Gamma_{\kappa}(\bm{\nu},\bm{v})$ is a Watson distribution jupp1989unified ; sra_karp:2013 , on $\mathbb{S}^{2}$

[TABLE]

where $M(\frac{1}{2},\frac{d}{2},\kappa)$ the confluent hypergeometric function also called the Kummer function ($d=3$ in our case) AbramowitzStegun1974 , $\pm\bm{v}$ the center of the distribution, and $\kappa$ the concentration parameter, which is roughly inverse proportional to the variance on the sphere. Because of the symmetry property of the Watson distribution and the antipodal symmetry of $\bm{I}(\bm{x},\bm{v})$, it is clear that $\bm{I}_{\sigma\kappa}(\bm{x},\bm{v})$ is antipodal symmetric too. As one alternative, a symmetrized von Mises-Fisher jupp1989unified distribution or a symmetrized heat kernel could be considered.

3.3.1 Action on Orientation

A diffeomorphic transformation $\varphi:\Omega\to\Omega$ acts on directions via its Jacobian: at $\bm{x}\in\Omega$, $J_{\bm{x}}\varphi$ sends the vector line ${\mathbb{R}}\bm{v}$, $\bm{v}\not=0$, to the vector line ${\mathbb{R}}J_{\bm{x}}\varphi(\bm{v})$, this is well-defined as $\det(J_{\bm{x}}\varphi)\not=0$. Thus, for each $\bm{x}\in\Omega$, $J_{\bm{x}}\varphi$ gives rise to projective linear transformation ${\mathbb{P}}J_{\bm{x}}\varphi$. Using the representation ${\mathbb{P}}^{2}\simeq\SS^{2}/\{\pm 1\}$, we can write ${\mathbb{P}}J_{\bm{x}}\varphi:\{\pm\bm{v}\}\mapsto\{\pm\frac{J_{\bm{x}}\varphi(\bm{v})}{|J_{\bm{x}}\varphi(\bm{v})|}\}$. We denote by $\psi$ the mapping $(\bm{x},\bm{v})\mapsto\frac{J_{\bm{x}}\varphi(\bm{v})}{|J_{\bm{x}}\varphi(\bm{v})|}$. Because it will only appear inside an antipodal symmetric kernel, $\psi_{\bm{x}}=\psi(\bm{x},\cdot)$ represents ${\mathbb{P}}{J_{\bm{x}}\varphi}$ without ambiguity.

3.3.2 Density, orientation and transformation

Because of spatial convolutions, we assume that the spatial domain $\Omega$ of $\varphi$ can be extended to ${\mathbb{R}}^{3}$, by taking it to be the identity out of a compact set $D$ containing $\Omega$. We set

[TABLE]

In the sequel, we in general omit $\bm{x}$ in $\psi_{\bm{x}}$, as the spatial location should be unambiguous. We write the joint histogram and density for similarity in image registration as

[TABLE]

[TABLE]

If one assumes that $W_{\alpha}$ is a Gaussian kernel with standard deviation $\alpha$ and we let $\alpha\rightarrow\infty$, we obtain full spatial integration, and we can write the joint histogram and joint probability density functions (PDF)

[TABLE]

[TABLE]

Similar formulas are obtained for single histograms$h_{\bm{K};\sigma\kappa\beta}$ and corresponding PDF $p_{\bm{K};\sigma\kappa\beta}$ with $K=\bm{I}\circ\Phi$ and $\bm{J}$, they are just the appropriate marginals of the joint distribution. This is the situation we consider in rest of this work.

3.4 LORD and Free Form Deformations

3.4.1 Free-form deformations.

In jensen2015locally , we assumed that $\varphi$ was a global affine transformation, with $\psi_{\bm{x}}\equiv\psi$ independent of $\bm{x}$, being the projectivization of the linear part of $\varphi$. In the present work, we assume $\varphi$ to be a non-rigid transformation, and we use Rueckert et al. framework rueckert1999nonrigid : the transformation $\varphi$ is given as a hierarchical spline representation affinely parameterized by a spatial grid ${\bm{c}}$ of control points. We follow rueckert1999nonrigid and describe briefly the construction and refer the reader to it and references therein for details. Our domain $\Omega$ is chosen to be a box $[0,N_{x}]\times[0,N_{y}]\times[0,N_{z}]$. At a given spatial resolution $r$, a 3D-grid $\bm{c}^{r}=(\bm{c}^{r}_{ijk})$ of size $(n^{r}_{x},n^{r}_{y},n^{r}_{z})$ of equally spaced (to start with) control points, with spacing distance $\delta_{r}$, is chosen. This defines a displacement function

[TABLE]

with $\bm{x}=(x,y,z)$, $i=\lfloor x/n^{r}_{x}\rfloor$-1, $j=\lfloor y/n^{r}_{y}\rfloor-1$, $k=\lfloor z/n^{r}_{z}\rfloor-1$, $u=x/n^{r}_{x}-i+1$, $v=y/n^{r}_{y}-j+1$, $w=z/n^{r}_{z}-k+1$ and $B_{l},B_{m}$ and $B_{n}$ are 1D cubic basis B-spline functions. We build $R$ grids of control points $\bm{c}=(\bm{c}^{1},\dots,\bm{c}^{R})$ from coarse to fine resolution, with corresponding $B^{r}_{\bm{c}^{r}}$. We set

[TABLE]

All the $B_{\bm{c}^{r}}^{r}$ depend linearly on $\bm{c}^{r}$, $\varphi_{\bm{c}}$ can be written as $\bm{x}+\bm{B}(\bm{x})\text{vec}(\bm{c})$, $\text{vec}(\bm{c})$ is a suitable vectorization of the control points, $\text{vec}(\bm{c})$ has dimensions $3k$, with $k$ the total number of control points, and $\bm{B}(\bm{x})\in{\mathbb{R}}^{3\times 3k}$. In the sequel, we will use $\bm{c}$ to denote both the grid of control points and its vectorization $\text{vec}(\bm{c})$, so as to avoid too cumbersome notations.

The spatial Jacobian $J_{\bm{x}}\varphi_{\bm{c}}(\bm{v})=\bm{v}+\left(J_{\bm{x}}\bm{B}\cdot\bm{v}\right)\bm{c}$. $J_{\bm{x}}\bm{B}$ is a tensor of size $3\times 3k\times 3$, which by (3.13) is composed of quadratic and cubic B-spline basis functions. The ’$\cdot$’ operation is the contraction ${\mathbb{R}}^{3\times 3k\times 3}\times{\mathbb{R}}^{3}$, $(r_{uvw},s_{w})\mapsto\sum_{w}r_{uvw}s_{w}$. We then have

[TABLE]

and set $\Phi_{\bm{c}}=(\varphi_{\bm{c}},\psi_{\bm{c}})$.

3.4.2 Registration objective.

We use the normalized mutual information (NMI) as a similarity measure

[TABLE]

The different terms are differential entropies of the joint density (3.12) and its marginals.

We add a regularizer in the form of a simple penalization on the non-uniformity of the control point grid by the squared difference between a point $\bm{c}_{i}$ and its direct neighbors. Recall that control points are organized as a family of $R$ grids, from coarse to fine resolution. The regularizer $\mathcal{S}(\bm{c})$ is the sum $\sum_{r=1}^{R}\mathcal{S}^{r}(\bm{c}^{r})$ at each resolution, where $\bm{c}^{r}=(c^{r}_{1},\dots c^{r}_{p_{r}})^{T}$ is the grid of control points at resolution level $r$. We denote by $N^{r}(i)$ the set of indices $j$ such that control point $c^{r}_{j}$ is neighbor to control point $c^{r}_{i}$ and by $|N^{r}(i)|$ its cardinal. We set

[TABLE]

$\lambda_{r}$ is a strictly positive parameter controlling the degree of smoothness. The negative sign in Equation 3.17 comes from the fact that we maximize the regularized NMI. We seek the transformation $\Phi_{\bm{c}^{*}}$ which maximizes the regularized NMI

[TABLE]

We use a quasi-Newton method to compute an optimum of objective function (3.18), L-BFGS from schmidt2005minfunc . We adapt the method to our multiresolution setting. At step $r\in\{1\dots,R\}$, all grid points except $c^{r}$ are kept fixed and we optimize our objective with respect to $c^{r}$.

This method requires the gradient $\nabla{\mathcal{M}}(\bm{c}^{r})$ with respect to the control points vector $\bm{c}^{r}$ for each resolution. The dependency of ${\mathcal{M}}$ in $\bm{c}^{r}$ appears, in the similarity term, deep inside its definition, this is reflected in its differential. In the next paragraph, we provide some details of its computation for the similarity term, and in Section 3.6, we derive the gradient of the regularization term.

3.5 Differential of the similarity term

In this paragraph , we omit the subscript $r$ and write $\bm{c}$ instead of $\bm{c}^{r}$. Computations are of the same nature across resolutions. The following diagram (Figure 1) exhibits the different PDFs that are used in the NMI term and illustrates its dependency in $\bm{c}$, which occurs only through $\Phi_{\bm{c}}$. Writing $NMI(\bm{I}\circ\Phi_{\bm{c}},\bm{J})$ as ${\mathcal{L}}(\bm{I}_{\sigma\kappa}\circ\Phi_{\bm{c}})$, its $\bm{c}$-differential can be written as

[TABLE]

(we use the notation $d_{\bm{I}}$ as there is in general no coordinate system on image spaces). A full formula for the Jacobian $J_{\bm{c}}NMI(\bm{I}\circ\Phi_{\bm{c}},\bm{J})$ (and thus gradient, as its transpose) is long and not very informative. Computations for $d_{\bm{I}}{\mathcal{L}}$, correspond to the red frame in Figure 1, and are obtained by application of standard derivation rules. They can actually be found in rueckert1999nonrigid ; darknersporring2012pami ; jensen2015locally . The computations in this paper which differ from the aforementioned correspond to the blue frame in Figure 1 and are now provided.

In our implementation, we replaced the continuous Watson kernel smoothing by a sum over a discrete set of $N$ directions $\nu_{1},\dots,\nu_{N}$ at each voxel $\bm{x}$. The normalization factor $M(\frac{1}{2},\frac{1}{3},\kappa)$ in the Watson kernel is replaced by its approximation from the discrete directions and is no longer a constant. The spatio-directional smoothing (3.5) is replaced by the semi-discrete smoothing

[TABLE]

where we have set

[TABLE]

The Jacobian of the spatio-directional smoothing with respect to the control point parameter $\bm{c}$ is:

[TABLE]

Because $\varphi_{\bm{c}}(\bm{x})=\bm{x}+{\bm{B}(\bm{x})}\bm{c}$, its Jacobian $J_{\bm{c}}\varphi_{\bm{c}}(\bm{x})$ is simply $\bm{B}(\bm{x})$ and, as $K_{\sigma}$ is a Gaussian kernel,

[TABLE]

The part involving directional kernel is somewhat more complex, as it involves $\psi_{\bm{c}}$ and thus the Jacobian $J_{\bm{x}}\varphi_{\bm{c}}$, as well as the normalizing factor in Equation 3.21. Set $f(x)=e^{\kappa x^{2}}$ and define $f_{i}=f(\nu_{i}^{\top}\psi_{\bm{c}}(\bm{v}))$. Since $f^{\prime}(x)=2\kappa xf(x)$, we have

[TABLE]

Then, by a straightforward calculation, we obtain

[TABLE]

The vector $\psi_{\bm{c}}(\bm{v})$ is the normalization of $\bm{V}:=J_{\bm{x}}\varphi_{\bm{c}}(\bm{v})$ $=\bm{v}+\left(J_{\bm{x}}\bm{B}\cdot\bm{v}\right)\bm{c}$. The differential of the normalization $\bm{V}\mapsto\frac{\bm{V}}{|\bm{V}|}$ is $\frac{1}{|\bm{V}|}\bm{\pi}_{\bm{V}^{\bot}}$ with $\bm{\pi}_{\bm{V}^{\bot}}$ the orthogonal projection on $\bm{V}^{\bot}$, the subspace orthogonal to $\bm{V}$. The Jacobian of $J_{\bm{c}}J_{\bm{x}}\varphi_{\bm{c}}(\bm{v})$ is simply $J_{\bm{x}}\bm{B}\cdot\bm{v}$ and

[TABLE]

Putting things together, the Jacobian with respect to the control point parameter $\bm{c}$ of the spatio-directional smoothing is given by

[TABLE]

The gradient of the similarity term is

[TABLE]

The operator $d_{\bm{I}}{\mathcal{L}}(\bm{I}_{\sigma\kappa}\circ\Phi_{\bm{c}})$ acts among other things by spatial, directional and intensity integration of the local changes described in $J_{\bm{c}}\left(\bm{I}_{\sigma\kappa}\circ\Phi_{\bm{c}}\right)$. It is rather complex and makes it difficult to provide an intuitive explanation for the different terms which compose Equation 3.27.

3.6 Regularization

In order to compute the gradient of $\mathcal{S}^{r}(\bm{c}^{r})$ defined in Equation 3.17, we introduce linear mappings $T^{r}_{i}:\bm{c}\mapsto|N^{r}(i)|c^{r}_{i}-\sum_{j\in N^{r}(i)}c^{r}_{j}$. The regularizer in Equation 3.17 can be rewritten as

[TABLE]

and by classical manipulation, we obtain that

[TABLE]

$T^{r*}_{i}$ is the adjoint of $T^{r}_{i}$, and given by $(T_{i}^{*}\bm{d})_{j}=N^{r}(i)\bm{d}$ if $j=i$, $-\bm{d}$ if $j\in N^{r}(i)$ and [math] otherwise. The resulting operator $\Delta_{r}=-\sum_{i=1}^{p_{r}}\frac{1}{|N^{r}(i)|^{2}}T_{i}^{r*}T^{r}_{i}$ is a discrete Laplacian. The sought gradient is

[TABLE]

In this work, we chose $\lambda_{1}=\dots=\lambda_{R}$.

3.6.1 Implementation details

To ensure that the deformation is a diffeomorphism, we check the Jacobian of the deformation for negative eigenvalues. In case of a negative eigenvalue of the Jacobian we perform backtracking, which in practice is performed by L-BFGS internally.

4 Experiments

To illustrate the properties of the LORD method, we conduct a series of experiments on simulated data and artificially warped real data.

4.1 Simulated Examples

Artificially generated samples enable visual inspection and validation of the framework in a highly controlled environment. While artificial experiments can be found in most DWI registration papers, artificial data is rarely available. To our knowledge, there are no open sources of simulated DWI data for comparing registration frameworks, although the DIPY project may be a good source for generating simulations garyfallidis2014dipy . To this end, we created our own simulated DWI data, which is freely available. The generation of the artificial data is described next.111Contact us or henrikgjensen@gmail.com for the code or examples.

4.1.1 Simulating DWI data

The simulated DWI (HARDI) samples were created by deforming a unit sphere of equally distributed directions semechko2012suite to a certain HARDI or ODF shape. Figure 2 shows two simulated HARDI samples and their corresponding ODFs; the first is a single fiber ODF, and the second is a crossing fiber ODF. DWI samples are antipodal symmetric and every ODF from the simplest to the most complex can be constructed through a combination of single fiber ODFs. The simulated data is visualized using the regularized q-ball imaging (QBI) algorithm, which uses a linear combination of real spherical harmonics to represent either the direct QBI sample or transformed ODF descoteaux2007regularized .

These models can be combined to form simulated DWI tracts in various DWI shapes, such as crossing fiber tracts Figure 3. A $20\times 20$ grid is used throughout this section to create blueprints of fiber tract constellations to apply LORD. The images are colored according to the generalized FA (GFA) value where dark blue regions represent free isotropic diffusion. To simulate a more realistic DWI scenario, random uniform noise has been added to the samples in Figure 3.

The simulated voxels with unit density are rescaled for the free diffusion to have a low density or mean diffusivity to resemble real $b_{0}$ normalized data.

4.1.2 Parametric Setup

For consistency, the same setup is used for all experiments on the simulated data, unless specifically stated otherwise.

Hierarchical mesh resolution.

In the B-spline deformation, the spacing between the control points is decreased in order to iteratively increase the degrees of freedom in the registration. We use 4 spacings $\delta_{1}=4$, $\delta_{2}=3.5$, $\delta_{3}=3$ and $\delta_{4}=2$. 10 iterations are used at all resolutions except the last, which terminates based on the optimal tolerance of $\varepsilon=10^{-6}$, or 90 iterations.

Watson concentration, directional resolution.

The concentration parameter is set to $\kappa=15$, which is sufficiently smooth to represent the 100 uniform directions used.

Spatial resolution.

A full spatial resolution is used with a B-spline smoothing at a near-Gaussian variance of $\sigma=0.6$ voxels.

Histogram size.

Due to the relatively small data sample, a low number of bins is used for the histogram (i.e. intensity resolution) $20\times 20$. This allows for some larger, more stable, but less refined deformations darknersporring2012pami .

HARDI registration, ODF visualization.

We register the raw HARDI models, but we visualize the ODFs of the warped data based on the Funk-Radon transform (FRT). We do this to illustrate that the warped raw data is correctly reoriented and do not visibly suffer from affine shearing.

Regularization.

We set the weight $\lambda$ to $10^{-4}$. Each of the examples below could be highly improved by an experiment-specific choice of parameters. However, for consistency, the same set of parameter was used for all experiments.

4.2 Experiment 1: Single Fiber Tracts

The first set of experiments is based on artificially generated distributions of HARDI shells and corresponding ODFs, each representing different fiber constellations. The experiment maps a straight fiber tract to a curved tract of the same width. This allows us to discuss some of the differences between a good reconstruction and a correct mapping.

4.2.1 Straight and curved

Three experiments illustrate different scenarios where straight fibers are registered to wavy fibers. These experiments illustrate the regularizing effects of using the full diffusion profile for registration. In the first experiment, the position of the end of the fibers is constrained by adding intersecting fibers at the end and beginning. Figure 5 shows the simulated moving (straight) and target (wavy) images.

The results of the first experiment are shown in Figure 6, where Figure 6a is the reconstructed warp, and Figure 6b shows the final spatial mapping from the moving image, overlaid on the original target image. As the figures show, the straight fibers are stretched in correspondence with the curvature, and the reconstructed ODFs from the HARDI-based registration are rotated correctly, although smoothed by the interpolation.

The second experiment is performed without features on the boundaries. The results are shown in Figure 7, where we observe that the length is preserved as the straight fiber tract is mapped to a sub-part of the curved tract. Intuitively, a correct registration in such a case should preserve the length of the straight fiber after deformation. The fact that this happens with no strong outside regularization, seems to indicate an inherent regularization effect of the cost function. In the third experiment, the proposed similarity is compared with the equivalent scalar-based registration by performing a mean diffusivity registration ($\kappa=0$) of the tracts with no signal on the boundaries. The mean diffusivity carries no directional information. The result are shown in Figure 8.

This pure scalar-based registration is driven by the edges of the simulated tracts only. The reconstruction appears to be correct, but the final spatial mapping indicates a lack of regularization as the fibers are stretched unevenly.

4.3 Experiment 2: Crossing Fiber Tracts

The second set of experiments are designed to investigate the registration of crossing fiber tracts.

4.3.1 Straight and shifted

The first experiment examines the framework’s ability to register two crossing tracts with a horizontal and vertical shift Figures 9a and 9b. Circular fibers have been added as fixed points in the image to illustrate the local shift of the crossing tracts. The result of the registration is shown in Figures 9c and 9d. The final spatial mapping, shown with arrows, is accurately including the reconstruction. Note that the reconstruction is subject to smoothing effect in the interpolation.

4.3.2 Two degrees of shearing

The second experiment involving crossing fibers shows three fiber tracts crossing under a varying amount of shear with a fixed horizontal tract (Figures 10a and 10b). The purpose is to investigate a relatively large deformation combined with a change to the complex crossing at the center The results are shown in Figures 10c and 10d. As the figure shows, the structure of the complex center-crossing closely matches the orientations of 45 degree crossing fibers. We remind the reader that the registration was not performed on the ODFs but directly on the simulated HARDI models and subsequently reconstructed.

4.4 Experiment 3: Fanning Fiber Tracts

The third set of experiments examines the most complex configurations namely the registration of kissing and fanning fibers. Note that the artificial cases constructed in these experiments are such that moving and target images are in a 1-1 correspondence.

4.4.1 Fanning and kissing fiber tracts in a crossing

The first experiment consist of two DWI images, that simulate both fanning (dispersing) and kissing (interleaving) fiber tracts. The moving image in Figures 11a and 11b contains a crossing with a few fibers fanning in and out along the vertical centerline. The target image contains two curved tracts fanning in and out, and merging at the central crossing. The results of the registration experiment are shown in Figures 11c and 11d, where both the reconstructed warp and spatial mapping show a registration that follows the lines of the original target image. The resulting registration can move, shrink, and turn the center-crossing to fit the target image.

4.4.2 Kissing fiber tracts

The second experiment involves two straight fiber tracts that are registered to two curving and interleaving tracts (Figures 12a and 12b). The result, shown in Figures 12c and 12d, displays a large and difficult deformation that by all accounts appears to be a successful registration.

4.5 Synthetic Deformations of Real Data

In this set of experiments, the registration framework is evaluated on real data obtained from the HCP van2013wu . By introducing a random synthetic warp on a subset of the brain data, ground truth is obtained where the goal is to register the warped image back to the original image. These experiments illustrate the effect of the parameters of the model in a realistic, but guaranteed diffeomorphic scenario.

4.5.1 DWI example data

The DWI data used in this experiment is shown in Figure 13, where the region of interest (ROI) is the blue overlay on the subject’s $b_{0}$ image. An ROI of the brain is used to improve the visualization of the results. The ROI was chosen at the edge of the corpus callosum (CC) in the left hemisphere due to the characteristic C-shape of the CC and the intersection with other well-known structures e.g. the cingulum. Furthermore, the ROI is in an area with crossing fiber tracts and is near the cortex. A $b=1000$ DWI volume is used with a ROI size $11\times 71\times 41$ at 1.25mm isotropic voxels with 90 directions. Only the central sagittal slice along with corresponding ODFs is visualized (Figure 14a), while the deformation is applied to the whole ROI. The deformation field for the central slice is shown in Figure 14b.

4.5.2 Parameters effect

The experimental setup provides a ground truth: the identity transformation. As the similarity measure, NMI, does not reflect a unique point-wise correspondence, we introduce three measures for evaluating the resulting registration: (1) mean squared error (MSE) between coordinates, (2) curl of the final deformation field, and (3) divergence of the final deformation field helmholz1858 .

Curl quantifies the amount of orientational change in the deformation field, also referred to as rotation or vorticity. We use the $L_{2}$-norm of the curl vector, i.e. the rotation velocity. A higher magnitude of the curl is expected for the results of a scalar-based registration compared to a registration over the full diffusion profile.

Divergence quantifies the density of the outward flux of a vector field which can be either positive or negative, indicating an expansion or a contraction at a given point. In combination, these three measures provide information about the magnitude of the voxel-wise distance and the state of the final deformation field. Curl and divergence have previously been used as regularization in a nonrigid registration framework riyahi2014regularization .

We report now results of experiments where we varied the intensity scale, the orientation scale, and the spatial scale: the intensity and orientational scales are the elements which distinguish LORD from the other registration frameworks. We investigated the following:

1. Isointensity curves.

The density-based formulationallows us to smooth the image inversely proportionally to image gradients, such as borders near the cerebrospinal fluid region, which can be strongly affected by partial volume effects. We iteratively decrease the control point-spacing in the FFD model to get an increasingly localized registration and increase the intensity range accordingly. This is realized with an initial histogram with relatively few bins and a fixed degree of smoothing and then successively increase the number of bins as the deformation resolution is increased.

2. Explicit reorientation.

The directional information is expected to provide a more stable and regularized transformation. To investigate this hypothesis, different levels of directional smoothing are examined from $\kappa=30$, a highly peaked ODF, to a scalar-based mean diffusivity registration, i.e. setting the concentration parameter $\kappa$ to [math]. The experiments on simulated data indicated that the directional information results in direct regularized solutions through the similarity, and these experiments seek to uncover if these observations hold for real data.

Performing multi-resolution registration (a continuation method) is a key element in most high resolution registration frameworks, e.g. in the FFD model rueckert1999nonrigid , ANTs avants2009advanced , Elastix klein2010elastix , FSL jenkinson2012fsl , etc. It provides stability while allowing for large and small deformations and reducing the computational complexity. While experiments regarding the spatial part will be performed, its impact on registration quality is fairly well-described in contrast to iso-intensity curves and explicit reorientation.

4.5.3 Parametric Setup

Hierarchical mesh resolution.

Similar to the simulated experiments, the spacing between the control points is successively decreased. To this end we use the composition $\delta_{1}=10$, $\delta_{2}=5$, $\delta_{3}=3.5$, $\delta_{4}=3$, (the control point spacing $\delta$ is scaled down according to the dimensions of the image). The low initial resolution of the deformation permits the generation of a large random and valid synthetic warp at around half of the control point spacing ($\sim$ $0.4\delta$) rueckert2006diffeomorphic . For the optimization, we use a quasi-Newton L-BFGS optimizer which runs for 50 iterations for all resolutions, unless the optimality condition of $\varepsilon=10^{-6}$ is fulfilled. Note that we refer to the result after each optimizer termination as a step, thus 4 steps in total.

Accumulated curl and divergence.

While the MSE iseasy to calculate at each step, the curl and divergence depend on the first-order derivatives of the spatial deformation, which are accumulated overeach step. Thus the final curl and divergence is defined from the product of Jacobians for spatial resolution Curl(\Phi_{\delta=3})\!=\!Curl(\Phi_{\delta=10})+Curl(\Phi_{\delta=5})+$$Curl(\Phi_{\delta=3.5})+Curl(\Phi_{\delta=3}) for all voxels.

Change in resolution.

Changing the spatial resolution is obtained by convolution.

Other fixed parameters.

The regularization is fixed to $\lambda=10^{-4}$, and we interpolate and optimize over 30 interpolated orientations out of the 90 ones present in the HCP data.

All the following experiments will show the results of the four accumulating steps. All error measures are reported for the entire ROI and not just the slice visualized.

4.5.4 The Quantitative Effect of Isoparametric Curves

This experiment illustrates the effect of changing the size of the smoothed joint histogram. Three experiments are performed with (i) a fixed histogram of $50\times 50$ bins, (ii) a histogram with $500\times 500$ bins, and (iii) gradually increasing the histogram size in 50, 100, 200, and 500 bins. Figure 15 shows the results in terms of MSE, curl, and divergence as mean value over all points along.

It is clear that a small histogram with wide isocurves results in an initially faster convergence Figure 15 (blue line). However, wide iso-curves result in flat regions with small gradients and little structure, which may cause the result to become sub-optimal or degenerate with increasing degrees of freedom in the transformation. In contrast, starting with a high bin resolution histogram with thin isocurves has the opposite effect in the initial steps (red line). However, by iteratively refining the joint histogram, we can refine image motion, from large to small, and this provides a better result (green line).

4.5.5 The Quantitative Effect of the Orientation Scale

In this experiment, we investigate whether directional information increases the stability and improves the registration. We employ the size progression of the histogram described in the previous experiment (i.e. $50\times 50$, $100\times 100$, $200\times 200$ and $500\times 500$ bins), and we vary the concentration parameter $\kappa$ from mean diffusivity at $\kappa=0$ to sharp angular features at $\kappa=30$. The results are shown in Figure 16.

The figure shows how the directional information results in better registration, with significantly less curl than the scalar registration at $\kappa=0$ and we observe that the best value $\kappa=30$ is stable in high directional resolution data such as the HCP. We suggest setting $\kappa=15$ as this should be enough for both high and low-resolution data, such as DTI.

4.5.6 The Quantitative Effect of the Spatial Resolution

In the last experiment, we use $\kappa=15$, set again the bins to $[50\times 50,100\times 100,200\times 200,500\times 500]$, and investigate the effects of the spatial scale.

The spatial resolution is set to $s=[4,3,2,1]$, which is equivalent to smoothly interpolating every 4th point, followed by every 3rd, etc. As with the control point spacing of the deformation field, we scale this to fit the image if the image is not cubical. For instance, if the image has the spatial dimension $100\times 150\times 50$ then space between spatial interpolations for $s=3$ will be $[2,3,1]$ with a bound on no less than 1 (i.e. full resolution). The results are shown in Figure 17. The results are very similar in the final step. However, the hierarchical resolution approach compared to the full resolution gave a speed-up factor of 2.4.

4.5.7 A Qualitative Example of the Results

Finally, we visualize the registration so as to perform a qualitative assessment of the warp, shape, and orientation of the individual ODFs. All registrations are based on the raw, noise-corrected HARDI data, while we show the tomographic inversion (FRT) indicating the direction of the diffusion and likely fiber tract orientations. Unlike the simulated experiments, we have fitted a B-spline to the image before deforming and visualizing the result, in contrast to the smoothing spline used in the artificial examples. The results where obtained using an increasing number of control point resolution with $\kappa=15$, same bin resolution refining and spatial resolution refining as above. Figure 18 shows the result for the central ROI slice, along with a zoomed-in version in Figure 19.

5 Conclusion

We have presented a scale-space formulation of density estimation that extends LOR to spatio-directional data, explicitly for the registration of DWI data with an explicit reorientation of the full diffusion profile. We have provided empirical evidence showing that the underlying structure of the data is preserved during registration while providing excellent registration results through many classical artificial examples known to be difficult for registration. We have shown that the formulation of the similarity itself provides regularization through the additional information given by the orientational dimension, illustrated clearly in the artificial examples. We have investigated the different scales provided by the framework and shown how the different parameters influence the registration results. In conclusion, LORD provides a smooth and well-matched deformation of DWI data, and the registration results are improved by the induced regularization obtained by integrating the orientational information in the objective function. Finally, we have shown that the deformation has very little impact on the shape of the deformed ODF.

Acknowledgements

This research was supported by Center for Stochastic Geometry and Advanced Bioimaging, funded by grant 8721 from the Villum Foundation.