Keywords

1 Introduction

The combined analysis of diffusion and structural MRI is extensively used in adult and neonatal [20] brain studies. Structural MRI has the highest contrast for the cortex region, while dMRI primarily provides information about white matter (WM) structures.

The uncertainty of deformation fields in the regions characterised by low contrast or homogeneous intensities (e.g., low WM fibre density regions in dMRI) is one the primary challenges associated with both longitudinal and inter-subject registration. Multi-channel registration that includes both anatomical and diffusion channels has been shown to improve registration and label-propagation results [2, 7, 19]. The reported MC registration solutions generally employ fractional anisotropy (FA) [7, 8, 14, 19] or DTI [2, 9, 13] as an additional channel. However, DTI-extracted metrics are characterised by inconsistencies in fibre-crossing regions. On the other hand, higher-order techniques such as constrained spherical deconvolution (CSD) [21] alleviate some of the limitations of the DTI model and allow extracting orientation-resolved microstructural information as so-called orientation distribution functions (ODFs) from HARDI data.

The classical approach for the fusion of information from different channels is based on simple averaging of individual channel updates [2]. More recently proposed solutions include scalar weighs for ROIs defined by thresholded FA maps [13] or local certainty maps based on normalised gradients correlated to structural content [7]. While the detailed overview of the choice of registration metrics is out-of-scope of this work, it can be summarised that the published works on intensity-based multi-channel registration primarily use the sum of squared differences (SSD) [2, 7, 8, 14] or local normalised cross-correlation (LNCC) [4] metrics. There is also a reported approach for T1-DTI atlas generation where datasets are spatially normalized only according to the structural channel [9].

Contributions. In this work we present a framework for multi-channel brain registration that allows local certainty-based fusion of dMRI-derived ODFs and structural MRI and is based on a novel similarity metric for dMRI. The solution is an extension of the multi-contrast ODF registration framework [15, 17]. The novel elements include implementation of local angular correlation (AC) as a metric for ODF channels, LNCC for structural channels and weighted fusion based on local certainty maps. The pipeline was implemented in MRtrix3 [22].

The method is evaluated on 20 longitudinal (intra-patient) neonatal MRI datasets from the developing Human Connectome Project (dHCP)Footnote 1 which constitutes a particularly challenging task for registration due to the rapid changes that occur in volume, structure and intensities during brain development. In addition, we demonstrate an example of a MC template of neonatal brain generated from 10 datasets (40–43 weeks PMA) using the proposed registration approach.

2 Method

2.1 Datasets, Acquisition and Pre-processing

The data used for evaluation of the proposed method include 20 longitudinal datasets of neonates scanned as a part of the dHCP project at St. Thomas Hospital, London. The gap between the scans is in the range of 7–11 weeks which is associated with significant changes in volume, myelination and cortical folding [16]. The postmenstrual age (PMA) at the first scan is within 30–35 weeks.

Each dataset includes two scans with diffusion and structural MRI volumes acquired on a Philips 3T scanner. The multi-shell HARDI volumes were acquired with four phase-encode directions on four shells with b-values of 0, 400, 1000 and 2600 s/mm\(^{2}\) with TE 90 ms, TR 3800 ms [10] with \(1.5\times 1.5\times 3\) mm resolution and 1.5 mm slice overlap and reconstructed to 1.5 mm isotropic resolution using the SHARD pipeline [5]. The structural T2-weighted volumes were acquired using a TSE sequence with TR = 12 s, TE = 156 ms. The T1-weighted volumes were acquired using an IR TSE sequence with TR 4.8 s, TE 8.7 ms. The isotropic T2 and T1 volumes with 0.5 mm resolution were reconstructed using a combination of motion correction [6] and super-resolution reconstruction [11]. All volumes of the same modality were normalised to the same global intensity ranges. The tissue segmentations were generated by the Draw-EM pipeline [12].

The preprocessing of the datasets was performed in MRtrix3 including: (i) decomposition of WM ODF from HARDI data via constrained spherical deconvolution (CSD) [21] followed by intensity normalisation; (ii) extraction of FA and mean diffusivity (MD) DTI-metrics; (iii) alignment of the structural to dMRI volumes based on affine registration of T2 to MD volumes using global NCC metric; (iv) resampling of all channels to 1 mm isotropic resolution with B-Spline interpolation. In addition, we manually segmented internal capsules (IC) in FA volumes for all datasets.

2.2 Multi-channel Registration Pipeline

The proposed registration pipeline is an extension of the multi-contrast ODF registration framework [15, 17]. The original method is based on SyN Demons [2] with an SSD metric and reorientation of ODF using apodized point spread functions [18]. In order to decrease the sensitivity to acquisition or physiology related changes in signal intensities, we replace SSD with a novel similarity metric based on angular correlation [1] and add certainty-maps weighting for fusion of structural and diffusion channels.

The input channels for each of the cases include: WM ODFs, structural MRI (T2-weighted and T1-weighted) volumes and FA maps. At first, the cases are globally aligned using affine registration of structural volumes using the global NCC metric. Next, we employ symmetric diffeomorphic LNCC demons [3] for structural and FA channels and local angular correlation metric [1] for ODF channels. In comparison to the classical ODF registration approach based on SSD metric in [15], using AC provides a more robust solution since it is less susceptible to the local changes in signal intensities while preserving directional information. However, unlike SSD, AC might be affected by the noise in the directional information.

Angular correlation \(r_{A}\) between two ODFs \(F^{ODF}\) and \(G^{ODF}\) represented with real valued spherical harmonic (SH) orthonormal basis functions \(Y_{lm}(\theta ,\phi )\)

$$\begin{aligned} F^{ODF}(\theta ,\phi ) = \sum _{l=0}^{\propto }{\sum _{m=-l}^{l}}{f_{lm}}Y_{lm}(\theta ,\phi ),\quad G^{ODF}(\theta ,\phi ) = \sum _{l=0}^{\propto }{\sum _{m=-l}^{l}}{g_{lm}}Y_{lm}(\theta ,\phi ) \end{aligned}$$
(1)

is computed as [1]:

$$\begin{aligned} r_{A} = \frac{{\sum _{l=2}^{L}}{\sum _{m=-l}^{l}}{f_{lm}g_{lm}}}{[{\sum _{l=2}^{L}}{\sum _{m=-l}^{l}}{|f_{lm}|^{2}}]^{\frac{1}{2}}[[{\sum _{l=2}^{L}}{\sum _{m=-l}^{l}}{|g_{lm}|^{2}}]^{\frac{1}{2}}}, \end{aligned}$$
(2)

where \(g_{lm}\) and \(f_{lm}\) are the SH coefficients of \(G^{ODF}(\theta ,\phi )\) and \(F^{ODF}(\theta ,\phi )\) of order L with even \(l=\{2,4,...,L\}\) harmonic degree terms, correspondingly. The \(l=0\) term does not contribute to AC values.

Since this is a correlation metric, the corresponding symmetric updates to the displacement fields \(\varLambda ^{F}\) and \(\varLambda ^{G}\) can be computed in a similar manner to LNCC demons [3] but with respect to the 4D ODFs rather than only the 3D local neighbourhood (Eq. 3).

$$\begin{aligned} \varLambda ^{F} = \frac{2FG}{F^{2}G^{2}}\left( G - \frac{FG}{F^{2}}F\right) \nabla F, \quad \varLambda ^{G} = \frac{2FG}{F^{2}G^{2}}\left( F - \frac{FG}{G^{2}}G\right) \nabla G, \end{aligned}$$
(3)

where \(G=\{g_{lm}^{n}\}_{{l=2,...,l_{max},m=-l,...,l}}\) and \(F=\{f_{lm}^{n}\}_{{l=2,...,l_{max}, m=-l,..., l}}\) are the vectors of SH coefficients at a given location in the 3D volume space with local neighbourhood \({n=1,...,N}\). We refer to this registration metric as local angular cross-correlation (LAC).

The updates from the structural channels are computed similarly to [3]. We also consider \(Y_{00}(\theta ,\phi )\) as a separate channel and use the LNCC metric for its contributions since it is excluded from the AC metric formalisation (Eq. 2).

The contributions from each of the channels i to the global symmetric displacement field update \(\varLambda ^{global}\) are locally weighted with respect to the 3D certainty maps based on the approach proposed in [7]. At first, the certainty maps \(\alpha _{i}^{F}\) and \(\alpha _{i}^{G}\) are computed from the original volumes F and G for each of the channels (including structural and ODF volumes) and normalised as:

$$\begin{aligned} \alpha _{i}^{F}\,=\,\parallel \nabla F_{i}^{T} \nabla F_{i} \parallel , \quad \widehat{\alpha _{i}}^{F} = \frac{\alpha _{i}^{F}}{max (\alpha _{i}^{F})} \end{aligned}$$
(4)

Then, the global symmetric updates to the displacement fields are computed by weighted averaging of the channel-specific update fields with respect to the certainty maps:

$$\begin{aligned} \varLambda _{global}^{F} = \frac{\sum _{i}\widehat{\alpha _{i}}^{F}\varLambda _{i}^{F}}{\sum _{i}\widehat{\alpha _{i}}^{F}}, \quad \varLambda _{global}^{G} = \frac{\sum _{i}\widehat{\alpha _{i}}^{G}\varLambda _{i}^{G}}{\sum _{i}\widehat{\alpha _{i}}^{G}} \end{aligned}$$
(5)

Figure 1 shows an example of the certainty gradient maps \(\widehat{\alpha _{i}}\) for structural, FA and one of the ODF component channels and the \(\sum _{i}\widehat{\alpha _{i}}\) of all channels. This approach ensures that the output deformation fields are defined by the contribution of the local channel regions with the highest structural content. This is relevant for the ROIs where one of the channels has low intensity contrast. In comparison, the multi-variate SyN (MVSyN) approach [2] is based on averaging of the individual channel updates.

Fig. 1.
figure 1

An example of the certainty maps \(\widehat{\alpha _{i}}\) for T2, one of the ODF component channels (\(l_{max}=2\)), FA, and the sum of all channels (\(\sum _{i}\widehat{\alpha _{i}}\)).

Full size image

2.3 Implementation Details

The method was implemented in MRtrix3 [22]. The new elements include: LAC metric for ODF registration and certainty-based weighting of the channels. In addition, we transferred ANTsFootnote 2 implementation of LNCC Demons metric [3] to MRtrix3 for registration of structural, \(Y_{00}(\theta ,\phi )\) ODF and FA channels.

It was experimentally identified that multi-resolution \(\{0.5; 0.75; 1.0\}\) and SH order \(l_{max}=\{0; 2; 4\}\) schemes and 3 voxel radius for the local neighbourhood for both structural and ODF channels are optimal for deformable registration of the investigated datasets. We used the standard MRtrix3 regularisation of gradient update and displacement fields based on Gaussian smoothing with 1 voxel standard deviation. The MRtrix3 parameter settings employed for generation of the multi-channel template are based on the pipeline formalised in [16].

3 Experiments and Results

3.1 Longitudinal Registration Study

For each of the investigated 20 cases, we performed a set of longitudinal (intra-patient) registrations with different settings including different combinations of channels \(\{\)T2; FA; T2+FA; STR (structural: T1+T2); ODF; ODF+STR; ODF+STR+FA\(\}\) and similarity metrics \(\{\)SSD; LNCC; LAC\(\}\). The channel weighting options include average and weighted: \(\{\)A-; W-\(\}\). The employed parameter settings are given in Sect. 2.3.

The MRtrix3-based implementation of LNCC Demons [3] is based on the ANTs toolbox and provides similar performance for structural registration. Therefore, we compare the proposed method directly to the existing MRtrix3 registration module. The main aim is an improvement of the combined quality of label propagation and image similarity for the structural and diffusion channels.

Table 1 presents the results of the comparison study. The quantitative evaluation is performed with respect to the quality of label propagation (Dice score) and similarity of the registered ODF volumes. The labels include: cortical grey matter (C-GM), hippocampus (HIP) and internal capsule (IC). The intensity-based similarity is assessed in terms of ODF AC (Eq. 2) for \(l_{max}=4\). The best performance results are highlighted in blue.

Firstly we can observe that locally weighted fusion [7] of T2 and FA improves the combined results while slightly lowering the dice score for cortex compared to single channel T2.

Table 1. Quantitative evaluation of the proposed multi-channel registration approach on longitudinal dMRI datasets of 20 neonates: Dice coefficient for brain tissue labels and AC between ODF volumes.
Full size table
Fig. 2.
figure 2

An example of longitudinal intra-patient registration for 31 \(\Longrightarrow \) 42 weeks PMA datasets. Difference between the original and transformed WM ODF SH coefficients (\(l=2,m=0\)) for the classical ODF registration with SSD metric, weighted MC registration of T2+FA channels with LNCC metric and weighted MC registration of ODF+T1+T2 channels with LAC and LNCC metrics.

Full size image

In general, myelination and cortical folding occurring during 7–11 weeks period significantly change local intensities in both structural and diffusion MRI data [16]. Therefore, even though all input volumes were normalised, using SSD metric for ODF or structural MC registration leads to the lower quality results in comparison to the proposed LAC metric which produced statistically significant (p < 0.05) improvement for C-GM and HIP Dice scores and ODF AC, while there was not a significant difference in Dice scores of IC.

Figure 2 demonstrates an example of the original and transformed WM ODF SH coefficients (\(l=2,m=0\)) for longitudinal registration of 31 to 42 weeks PMA datasets. There is a clear difference in the magnitude of SH coefficients between the original scans. Using MC registration with LAC for ODF and LNCC for structural channels produces visually sharper results for the IC region in comparison to both classical SSD ODF registration [15, 17] or fusion of T2 and FA [7]. This is in agreement with the higher AC values reported in Table 1.

Fig. 3.
figure 3

An example of longitudinal symmetric MC registration for 31 \(\Longleftrightarrow \) 42 weeks PMA at scan case including: original and transformed T2 volumes, original and transformed ODF over T2, original and transformed ODF over masked original T1 volume (used as a template), original and transformed labels.

Full size image

Furthermore, there is a clear indication that additional structural channels (in this case T1+T2) and certainty-based weighting increase the quality of label propagation and AC similarity of ODF volumes. Adding the FA channel did not significantly affect the results since ODFs contain the WM structure information.

All ODF-based options resulted in approximately the same range for the IC Dice score values due to its high contrast and showed significant improvement (p < 0.05) in comparison to using the FA and T2 channels only. Apart from the IC values for ODF registration, the improvement in performance of the proposed method (ODF+STR: W-LAC/LNCC) in comparison to the baseline methods (structural LNCC Demons, ODF MRtrix registration as well as fused T2+FA) is statistically significant with p < 0.05.

An example of symmetric LAC+LNCC MC registration for 31 \(\longleftrightarrow \) 42 weeks PMA at scan case is presented in Fig. 3. The registration of the structural and ODF channels was successful even though there are significant differences in contrast of both structural and ODF volumes, cortex folding surface and the global shape. Visualisation of the original and transformed normalised ODFs over the same padded T1 volume (third row) confirms that the global shape and features of the volumes are sufficiently well aligned. Label propagation for tissue and IC segmentations also resulted in relatively similar results. This, however, might also be affected by the quality of the original segmentations produced by the automated Draw-EM method [12].

3.2 Multi-channel Template Example

Figure 4 shows and example of a multi-channel template generated using MRtrix3 \(population\_template\) tool with the proposed MC registration pipeline and LAC+LNCC metrics. The template with 1 mm resolution was generated from 10 neonatal MRI datasets from 40–43 weeks PMA. It includes T2, T1, normalised WM ODF and FA channels. The resulting volumes are characterised by well defined features of both cortex and WM structures.

Fig. 4.
figure 4

An example of multi-channel template of neonatal brain generated from 40–43 PMA datasets including: WM ODF, FA, T1 and T2 channels.

Full size image

4 Discussion and Conclusions

This paper presents a solution for multi-channel registration combining multi-shell HARDI and structural MRI data. It is based on a novel similarity metric for diffusion MRI and certainty-based weighting of the channels. The method was implemented in MRtrix3 and can be integrated into neuroimaging pipelines.

The quantitative evaluation was performed on 20 longitudinal neonatal datasets from the dHCP project. The results showed that fusion of structural and diffusion ODF channels improves overall results, compared to single-channel registrations. The weighting of channels based on certainty maps also improves the results thus potentially minimising the uncertainty of deformation fields. Furthermore, the proposed LAC metric outperforms the state-of-the-art ODF registration method for challenging cases.

An example of the generated multi-modal template shows that this tool has a potential application for generation of spatio-temporal multi-modal brain MRI templates that require robust similarity metrics. Simultaneous segmentation of WM and cortex structures could also potentially improve the accuracy of morphometry in structural MRI processing pipelines.

However, these results also emphasise the fact that accurate alignment of diffusion and structural volumes is a critical step for multi-channel registration since affine registration might not fully solve this due to distortions in dMRI data. Future work will focus on further optimisation of the MC ODF registration pipeline and extensive evaluation on adult and multi-site datasets.