- Email: [email protected]
Automated correction of improperly rotated diffusion gradient orientations in diffusion weighted MRI
Accepted Manuscript Automated correction of improperly rotated diffusion gradient orientations in diffusion weighted MRI Ben Jeurissen, Alexander Leemans, Jan Sijbers PII: DOI: Reference:
S1361-8415(14)00093-0 http://dx.doi.org/10.1016/j.media.2014.05.012 MEDIMA 897
To appear in:
Medical Image Analysis
Received Date: Revised Date: Accepted Date:
9 December 2013 20 May 2014 27 May 2014
Please cite this article as: Jeurissen, B., Leemans, A., Sijbers, J., Automated correction of improperly rotated diffusion gradient orientations in diffusion weighted MRI, Medical Image Analysis (2014), doi: http://dx.doi.org/ 10.1016/j.media.2014.05.012
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Automated correction of improperly rotated diffusion gradient orientations in diffusion weighted MRI Ben Jeurissena,∗, Alexander Leemansb , Jan Sijbersa a
b
iMinds-Vision Lab, Department of Physics, University of Antwerp, Belgium Image Sciences Institute, University Medical Center Utrecht, The Netherlands
Abstract Ensuring one is using the correct gradient orientations in a diffusion MRI study can be a challenging task. As different scanners, file formats and processing tools use different coordinate frame conventions, in practice, users can end up with improperly oriented gradient orientations. Using such wrongly oriented gradient orientations for subsequent diffusion parameter estimation will invalidate all rotationally variant parameters and fiber tractography results. While large misalignments can be detected by visual inspection, small rotations of the gradient table (e.g. due to angulation of the acquisition plane), are much more difficult to detect. In this work, we propose an automated method to align the coordinate frame of the gradient orientations with that of the corresponding diffusion weighted images, using a metric based on whole brain fiber tractography. By transforming the gradient table and measuring the average fiber trajectory length, we search for the transformation that results in the best global ‘connectivity’. To ensure a fast calculation of the metric we included a range of algorithmic optimizations in our tractography routine. To make the optimization routine robust to spurious local maxima, we use a stochastic optimization routine that selects a random set of seed points on each evaluation. Using simulations, we show Corresponding author Email: [email protected] Address: Vision Lab University of Antwerp (CDE) Universiteitsplein 1, N.1.18 B-2610 Wilrijk, Belgium Phone: +32 3 265 24 77 Fax: +32 3 265 22 45 ∗
Preprint submitted to Medical Image Analysis
May 31, 2014
that our method can recover the correct gradient orientations with high accuracy and precision. In addition, we demonstrate that our technique can successfully recover rotated gradient tables on a wide range of clinically realistic data sets. As such, our method provides a practical and robust solution to an often overlooked pitfall in the processing of diffusion MRI. Keywords: diffusion-weighted MRI, coordinate frame, rotationally variant parameter, fiber tractography, diffusion gradient reorientation, stochastic approximation 1. Introduction Diffusion-weighted (DW) MRI is a non-invasive imaging technique that is sensitive to the random microscopic motion of water molecules (Stejskal and Tanner, 1965). In the brain white matter (WM), which consists of tightly packed bundles of neuronal axons, the regular arrangement of fibers introduces a directional dependence of this motion (Moseley et al., 1990). The degree of this directionality provides an indirect measure for the organization of the underlying tissue, while the preferred diffusion orientation provides an estimate of the local fiber orientation of the bundles (Basser, 1995). These local fiber orientation estimates can be pieced together to infer the long-range pathways connecting distant regions of the brain, using a technique called WM fiber tractography (Conturo et al., 1999; Mori et al., 1999). DW-MRI is currently the only imaging technique that allows assessment of WM organization and connectivity in vivo and it is increasingly used in clinical studies to investigate brain WM, and its development and disorders (Johansen-Berg and Behrens, 2006; Mori and Zhang, 2006; Assaf and Pasternak, 2008). Diffusion metrics extracted from DW images can be categorized into two categories based on their rotational (in)variance. Rotationally invariant diffusion metrics embody intrinsic features of the local diffusion pattern and are independent of the orientation of the anisotropic structure, the orientation of the specimen with respect to the laboratory coordinate system, and ¨ even the choice of the laboratory coordinate system (Basser and Ozarslan, 2010). Notable examples are the scalar diffusion parameters obtained from diffusion tensor imaging (DTI) such as the widely used fractional anisotropy (FA) and mean diffusivity (MD) (Basser and Pierpaoli, 1996). Both parameters are extracted from the eigenvalues of the diffusion tensor, making them rotationally invariant by definition. 2
Rotationally variant diffusion parameters, on the other hand, embody orientational features of the local diffusion pattern and are dependent of the orientation of the anisotropic structure, the orientation of the specimen with respect to the laboratory coordinate system, and the choice of the laboratory coordinate system. Notable examples are the principal diffusion vector (PDV) estimated by DTI, diffusion/fiber orientation distribution functions obtained from high angular resolution diffusion imaging (HARDI) (Tournier et al., 2011) and fiber tractography results. For rotationally variant parameters to be meaningful, it is important that the coordinate system of the diffusion sensitizing gradient orientations and the spatial coordinate system of the voxels in the diffusion weighted data set coincide. In theory, the meta data generated by MRI scanner software should provide all the necessary information to properly align the gradient orientations with the subject coordinate frame. In practice, however, this procedure is error-prone, as the appropriate steps are dependent on the coordinate system conventions of the scanner software, the file format and even the post-processing tool. For example, while both MRtrix (Tournier et al., 2012) and FSL (Jenkinson et al., 2012) support DWI data in the NifTI format and while both represent the gradient vectors as simple text files, MRtrix assumes that the gradient vectors are provided with respect to real-world/scanner coordinates, whereas FSL assumes that the gradient vectors are provided with respect to the image axes. While large misalignments (e.g. flipping or permutations of the gradient orientation coordinate axes) can be detected by visual inspection of PDV maps or tractography results, small rotations of the gradient table (e.g. due to angulation of the acquisition plane), are much more difficult to detect. Nevertheless, even small rotations have a large impact on tractography results due to an accumulation of errors. In this work, we propose a data-driven method that does not rely on header information to automatically detect and correct for the rigid misalignment of the gradient table and the DW data, using a metric based on whole brain tractography.
3
2. Methods The use of incorrect gradient orientations will result in the premature termination of fiber trajectories, which can be measured using whole brain fiber tractography. By transforming the gradient table and measuring the average fiber trajectory length, one can search for the transformation that results in the best global ‘connectivity’. Finding the correct DW gradient settings then becomes akin to a registration problem (Figure 1). The gradient orientations (moving entity) are subjected to a 3D rotation (transformation), such that they are aligned with the DW images (reference entity). Reference entity: DW images
Metric: Average tract length
Moving entity: Gradient table
iterate
Transform: 3D rotation
fitness value
Optimizer: Exhaustive search / Stochastic gradient ascent
transform parameters
Figure 1: Schematic overview of the gradient table reorientation method.
2.1. Transform The 3D rotation R of the gradient table can be parameterized as successive rotations around the main orthogonal axes: R(θ) = Rx (θx )Ry (θy )Rz (θz )
(1)
resulting in a parameter vector θ = (θx , θy , θz ) with 3 transformation parameters. Note that, since the gradient table is rotated as a whole, this is a rigid body transformation. Also note that translations do not have to be taken into account, as they do not affect the gradient orientations.
4
2.2. Metric The metric, which provides a measure of how well the gradient orientations match the DW images, is defined as the average trajectory length f (θ, n), measured with whole-brain streamline DTI tractography, where θ parametrizes the 3D rotation of the gradient table and n is the number of seed points, evenly distributed over the whole brain. As fiber tractography is more reliable in structures with high FA, the length contribution of each step is weighted by the local FA value. To ensure robust fiber tractography even on data of low quality, the DWI images are smoothed with a 3D Gaussian filter prior to tractography. The metric exploits the fact that the brain is a complicated 3D network of fiber bundles. While rotating the gradient orientations away from the true orientations may result in a few slightly longer false trajectories in some locations, it will at the same time reduce the length of many true trajectories. Since finding the optimal gradient table requires many repeated evaluations of the metric, it is imperative that the computation of f (θ, n) is fast. However, as whole-brain tractography is relatively time-consuming, several optimizations were implemented to speed up the standard streamline DTI tractography algorithm. Traditionally, DTI tractography is performed as follows (Basser et al., 2000). First, the DW signal at the current position of the trajectory is obtained using trilinear interpolation. Next, the diffusion tensor is estimated from the interpolated signal. Then an eigenvalue decomposition is performed on the diffusion tensor and the eigenvector corresponding to the largest eigenvalue is assumed to be the current fiber orientation. Finally, the trajectory is advanced by a fixed step size along this orientation and the whole process is repeated. The most time consuming operations in this algorithm are the interpolation, the diffusion tensor fit and the eigenvalue decomposition. The computational complexity of the DTI tractography algorithm can be reduced significantly using the Log Euclidian framework (Arsigny et al., 2006), as follows. Prior to tractography, the log diffusion tensor field is calculated for the whole brain. Then, the log diffusion tensor at the current position of the trajectory can be obtained using trilinear interpolation directly on the log diffusion tensor field. Next, an eigenvalue decomposition is performed on the diffusion tensor and the eigenvector corresponding to the largest eigenvalue is assumed to be the current fiber orientation. Finally, the trajectory is advanced by a fixed step size along this orientation and the whole process is repeated. Following this approach, only a single diffusion 5
tensor fit is required and the number of interpolations at each tractography step is limited to 6 as compared to the number of diffusion weighted images for the traditional method, greatly reducing computation time and memory consumption. Instead of rotating the gradient table and recalculating the log diffusion tensor at each metric evaluation, our algorithm directly rotates the log tensor, requiring only a single diffusion tensor fit for repeated evaluations of the metric. To avoid unnecessary calculations outside the brain, a rough skull stripping was performed on the b = 0 s/mm2 images, separating brain from background. The only remaining bottleneck is the eigenvalue decomposition which needs to be performed at each tractography step. To reduce the computational load of this operation, the analytical eigenvalue decomposition proposed in (Hasan et al., 2001) was used, which exploits the fact that the diffusion tensor is a 3×3 symmetric and semi-positive-definite matrix, and which is approximately 40 times faster than standard iterative decomposition techniques. Finally, to reduce computation time even further, whole brain tractography is not performed from all N voxels of the whole brain, but from a set of n < N seed points distributed over the whole brain, allowing a trade-off between metric calculation speed and metric precision. For reasons outlined below, we also define a stochastic cost function (Robbins and Monro, 1951) f˜(θ, n) which is the approximate average tract length resulting from whole brain tractography with rotation angles θ and n the number of seed points randomly distributed over the whole brain according to a uniform distribution. As a consequence, each time f˜(θ, n) is called, it probes a different set of seed points, resulting in a different metric outcome. 2.3. Optimization We provide two optimization approaches, depending on what type of gradient corruptions are present in the data. The first approach assumes that the gradient table is only corrupted by means of column order permutations and/or column sign changes, a situation that often occurs in practice, due to different coordinate system conventions. The number of possible column order permutations is 3! = 6, while the number of possible sign changes is 4 (sign change in the first column, the second column, the third column and no sign change), resulting in 6 × 4 = 24 possible solutions in total. This optimization problem can be solved in a 6
brute force fashion by exhaustively evaluating the metric f (θ, n). From this list of evaluations fi (with i = 1, ..., 24) we select the permutation and/or sign change that results in the longest average fiber length. To avoid selecting false solutions in the case of insufficient signal-to-noise ratio (SNR), instead of just evaluating f (θ, n), we repeatedly evaluate f˜(θ, n) and the solution is only deemed optimal if it is consistently selected as having the longest average fiber orientation across the stochastic function evaluations. If the same solution is not found consistently across stochastic function evaluations, the method reports that the solution is undecided due to insufficient SNR. In that case, the method is rerun with twice the amount of seed points n, increasing the metric precision. This way, the method dynamically adjusts the number of seed points required for precise metric evaluation to the quality of the data and avoids false solutions. The second approach assumes that the gradient table is corrupted by means of a general 3D transformation. Note that, due to the antipodal symmetry of the diffusion signal, this approach also resolves the column order permutations and column sign changes resolved with the first approach. Since a reliable evaluation of the metric requires a large number of seed points n, repeated evaluation has a very high computational cost. In addition, even with high n, the metric still exhibits spurious maxima, which would cause deterministic optimization routines to report suboptimal solutions (Figure 2). In order to keep the computation time of the metric to a minimum and allow accurate estimation of the rotation parameters, we chose to use the stochastic gradient ascent (Robbins and Monro, 1951) which has previously been used successfully for the purpose of image registration (Klein et al., 2007). The stochastic gradient ascent method follows the same scheme as the deterministic gradient ascent, with the distinction that the derivative of the cost function is replaced by a fast stochastic approximation. By using a new, randomly selected subset of seed points every iteration of the optimization process, a bias in the approximation error is avoided. More specifically, in order to find the optimal set of rotation angles θ, a stochastic gradient ascent is defined as follows: θk+1 = θk + γk ∇f˜(θk , n)
with
k = 0, ..., K
(2)
with γk the scalar gain factor or “step size” and f˜(θk , n) the approximate average tract length resulting from whole brain tractography with rotation 7
f
30
n = 50 n = 500 n = 5000
20
10 −180
−120
0
−60
60
120
180
θx ( ) (a) f over the full domain of θx n = 5000
f
28
27.5
27 −5
0
−2.5
2.5
5
θx ( ) (b) f near the optimum Figure 2: Cost function values f for a phantom data set with correctly oriented gradient table and SNR = 10 (defined on the b = 0 s/mm2 images) as a function of θx for selected values of n (a). A zoomed-in plot near the optimum is provided in (b).
angles θk and n randomly distributed seed points. The approximated gradiˆ in contrast to the exact ent does not necessarily vanish close to the solution θ, ˆ = 0. Thus, convergence must be forced by enderivative that satisfies ∇f (θ) suring γk → 0 as k → K. In this work, K = 150 and γk ≡ γ(k) = 100e−0.068k . Using the stochastic gradient ascent makes the optimization routine robust to spurious maxima, and also allows the number of seed points n to be small, significantly speeding up the metric evaluation (Klein et al., 2007). Unless stated otherwise, n = 200 seed points randomly distributed over the whole brain, were used in all subsequent experiments. While a small number of 8
seed points results in f˜(θk , n) having low precision, the repeated evaluations of f˜(θk , n) will converge to the true optimum. Since the gradient ascent can only find the local minimum, we use a multi-start approach starting from 30 points distributed evenly in the space of 3D rotations. 2.4. Simulations To evaluate the proposed method quantitatively, a noiseless full-brain DWI data set consisting of 6 b = 0 s/mm2 and 60 b = 1200 s/mm2 images and with a voxel size of 2.4 × 2.4 × 2.4 mm3 was simulated (Leemans et al., 2005). From this data set, noisy Rician distributed data sets were obtained for different noise levels, with SNR in the b = 0 s/mm2 image ranging from 1 to 30 and 100 noisy realizations for each noise level (Gudbjartsson and Patz, 1995). 2.4.1. Gradient column order permutations and sign changes For each noisy realization, we checked whether the method could distinguish the correct column order permutation and sign changes from the wrong ones. We also investigated the amount of seed points required in the metric to make this distinction. 2.4.2. General 3D rotations For each noisy realization, the gradient table was corrupted by a random 3D rotation. After correction with the proposed method, quality of the corrected gradient table was assessed with respect to accuracy and precision as follows. For each gradient orientation, an average gradient orientation was obtained from the 100 corrected orientations by extracting the first eigenvector of the mean dyadic tensor (Basser and Pajevic, 2000) and the minimum angles subtended between each estimated gradient orientation and the average gradient orientation were recorded. A measure of accuracy was obtained by calculating the angular bias of the average gradient orientation compared to the ground truth gradient orientation (Basser and Pajevic, 2000). A measure of precision was obtained by calculating the 95% confidence interval (CI) of the one-tailed distribution of minimum angles between the average gradient orientation and the individual gradient orientations (Jones, 2003).
9
2.5. Experiments with in vivo human brain data The method was also tested on a set of real human brain DWI data sets with different characteristics: • Data set A is a DTI data set acquired on a 3T Siemens Trio scanner, using a 12-channel receiver head coil. Diffusion weighting was applied along 30 gradient directions with b = 1000 s/mm2 . Additionally, 2 nonDW images were acquired. Other imaging parameters were: TR/TE: 8100/116 ms; voxel size: 2.5 × 2.5 × 2.5 mm3 ; matrix: 96 × 96; slices: 54; NEX: 1. • Data set B is a HARDI data set acquired on a 3T GE HDx Signa scanner, using an 8-channel receiver head coil and cardiac gating (Jeurissen et al., 2013). Diffusion weighting was applied along 60 gradient directions with b = 1200 s/mm2 . Additionally, 6 non-DW images were acquired. Other imaging parameters were: TR: 20 R-R intervals; TE: 86 ms; voxel size: 2.4 × 2.4 × 2.4 mm3 ; matrix: 96 × 96; slices: 60; NEX: 1. • Data set C is a HARDI data set with high spatial resolution acquired on a 3T Philips Achieva scanner, using an 8-channel receiver head coil (Tax et al., in press). Diffusion weighting was applied along 60 gradient directions with b = 2500 s/mm2 . Additionally, 1 non-DW image was acquired. Other imaging parameters were: TR/TE: 10265/107 ms; voxel size: 2 × 2 × 2 mm3 ; matrix: 112 × 112; slices: 70; NEX: 1 . • Data set D is a multi-shell DW data set acquired with high spatial resolution on a 3T GE Discovery MR750 scanner, using an 8-channel receiver head coil. Diffusion weightings of b = 0, 700, 1200 and 2800 s/mm2 were applied in 5, 25, 45 and 75 directions, respectively. Other imaging parameters were: TR/TE: 9500/100 ms; voxel size: 2 × 2 × 2 mm3 ; matrix: 120 × 120; slices: 68; NEX: 1. We deliberately corrupted the gradient tables of these data sets using a random column permutation, a random column sign change, a small random 3D rotation and a large random 3D rotation. Subsequently, the corrupted data sets were processed using our correction method, to test whether it could correctly restore the gradient orientations.
10
3. Results 3.1. Simulations Image quality for the different simulated SNR levels can be appreciated from Figure 3. Note that SNR is defined on the b = 0 s/mm2 image and that for realistic clinical studies it is in the range of [10-30]. However, to ensure proper functioning of our method even on low quality data, we also simulated SNR values well below 10.
DWI
DEC
SNR
30
20
10
5
2
Figure 3: Illustration of the image quality at different simulated SNR levels, for a single DW image and for the resulting DEC FA map. Note that SNR is defined on the b = 0 s/mm2
11
3.1.1. Gradient column order permutations and sign changes Figure 4 shows the results of the exhaustive evaluation of the different combinations of gradient table column order permutations (rows) and sign changes (columns), for a single realization of different SNR levels. The correct answer is the top left cell, which corresponds to no permutation and no sign changes. Note that one can clearly distinguish the correct solution up till SNR = 5. However, at SNR = 2, the method is still able to repeatedly recognize the first solution as the optimum, over repeated stochastic evaluations of the metric f . At SNR = 1, no solution stands out from the rest and the method reports no answer. Figure 5 shows the percentage of correct, no, and wrong solutions as a function of SNR. Note that the method achieves a 100% success rate all the way up to SNR = 2. Starting from SNR = 1, the method reports that it can no longer reliably estimate the correct gradient orientations in all cases. Note that the method makes no false claims about the correct gradient orientations, even at low SNR. Figure 6 shows the number of seed points required by the method to restore the gradient table as a function of SNR. Note that 50 seed points were sufficient all the way down to SNR = 4. Dropping below SNR = 4, the method gradually requires more seed points. At SNR = 2, the method requires at most 200 seed points, which is still reasonable.
12
column permutations
column permutations
column permutations column sign changes
column sign changes
column permutations
(c) SNR = 10
column sign changes
column permutations
(b) SNR = 20
column permutations
(a) SNR = 30
column sign changes
column sign changes
(d) SNR = 5
column sign changes
(e) SNR = 2
5
10
15
20
(f) SNR = 1
25
30
Figure 4: Average metric values (over repeated stochastic metric evaluations) for the different combinations of gradient table column order permutations and sign changes, for different SNRs (single realization). The correct answer is the top left cell, which corresponds to no permutation and no sign changes.
13
100 correct solution no solution wrong solution
80
%
60 40 20 0
30
20
10
5
4
3
2
1
SNR Figure 5: Percentage of correct (green), no (blue) and wrong (red) solutions as a function of SNR, over 200 noise realizations.
100 80 50 100 200 400 800 1600 3200 6400
%
60 40 20 0
30
20
10
5
4
3
2
1
SNR Figure 6: The number of seed points required by the method to restore the gradient table as a function of SNR, over 200 noise realizations.
14
3.1.2. General 3D rotations Figure 7 shows the gradient orientations recovered with our method after random 3D rotations for 100 noisy realizations of the data set, together with the true gradient orientations. For SNR values all the way down to 5, the corrected gradient orientations are tightly clustered around the true gradient orientations, with no discernible bias. For SNR = 2, precision of the corrected gradient orientations is markedly decreased, but the corrected gradient orientations are still clustered around the true orientations, suggesting only a small bias. Note that an SNR as low as 2 in the b = 0 s/mm2 images, is extremely low and does no longer represent realistic clinical data.
(a) SNR = 30
(b) SNR = 20
(c) SNR = 10
(d) SNR = 5
(e) SNR = 2
(f) SNR = 1
Figure 7: The recovered gradient orientations (colored dots) vs the true gradient orientations (black crosses) for 100 noisy realizations of the data set. Corresponding corrected gradient orientations are given the same color to aid visual clustering.
Figure 8 quantifies both accuracy and precision as a function of SNR. Subdegree accuracy is obtained starting from SNR = 2. Sub-degree precision is obtained, starting from SNR = 4. This is well within the boundaries of clinically realistic SNR. 15
bias ( )
10 1
0.1
12345
10
20
30
20
30
SNR (a) Accuracy
95% CI ( )
100
10
1
0.1
12345
10 SNR (b) Precision
Figure 8: Accuracy (a) and precision (b) of the recovered gradient tables as a function of SNR. The dashed red lines indicate the points at which sub-degree accuracy (a) and sub-degree precision (b) are achieved.
16
3.2. Experiments with in vivo human brain data 3.2.1. Data set A Figure 9 shows an FA map, a DEC (directionally encoded color) FA map (Pajevic and Pierpaoli, 1999) and a whole brain tractogram for data set A, after corruption with a column order permutation (a), after correcting for column order permutations and column sign changes (b), after correction of a full 3D rotation (c), and without corruption (d). Note that column order permutation can easily be detected by visual inspection of the DEC FA map and causes a significant disruption of the tracts in the whole brain tractogram (a). Both our correction strategies correctly restore the DEC FA map and whole brain tractogram (b-d).
(a) Before correction
(b) After column permutation and sign correction
(c) After full correction
(d) Original
Figure 9: Example of the gradient correction on in-vivo data set A with a gradient table column permutation. FA, DEC and a whole brain tractogram are shown before correction (a), after gradient column permutation and sign correction (b), after the full correction (c) and for the original data set (d).
17
3.2.2. Data set B Figure 10 shows an FA map, a DEC FA map and a whole brain tractogram for data set B, after corruption with a column sign change (a), after correcting for column order permutations and column sign changes (b), after correction of a full 3D rotation (c), and without corruption (d). Note that column sign change can no longer be detected by visual inspection of the DEC FA map (a). However, the column sign change still causes a significant disruption of the tracts in the whole brain tractogram (a). Both our correction strategies correctly restore the DEC FA map and whole brain tractogram (b-d).
(a) Before correction
(b) After column permutation and sign correction
(c) After full correction
(d) Original
Figure 10: Example of the gradient correction on in-vivo data set B with a gradient table column sign change. FA, DEC and a whole brain tractogram are shown before correction (a), after gradient column permutation and sign correction (b), after the full correction (c) and for the original data set (d).
18
3.2.3. Data set C Figure 11 shows an FA map, a DEC FA map and a whole brain tractogram for data set C, after corruption with a large, random 3D rotation (a), after correcting for column order permutations and column sign changes (b), after correction of a full 3D rotation (c), and without corruption (d). Note that the large rotation can detected by visual inspection of both the DEC FA map and the whole brain tractogram (a). However, manual realignment is now no longer feasible. Correcting for column order permutations and column sign changes only marginally improves the DEC FA map and the whole brain tractogram, and cannot completely restore it (b). The full 3D correction on the other hand results in a correctly restored DEC FA map and whole brain tractogram (c-d).
(a) Before correction
(b) After column permutation and sign correction
(c) After full correction
(d) Original
Figure 11: Example of the gradient correction on in-vivo data set C with a gradient table subject to a large 3D rotation. FA, DEC and a whole brain tractogram are shown before correction (a), after gradient column permutation and sign correction (b), after the full correction (c) and for the original data set (d).
19
3.2.4. Data set D Figure 12 shows an FA map, a DEC FA map and a whole brain tractogram for data set D, after corruption with a subtle, random 3D rotation (a), after correcting for column order permutations and column sign changes (b), after correction of a full 3D rotation (c), and without corruption (d). Note that the small rotation is virtually impossible to detect by visual inspection on both the DEC FA map and the whole brain tractogram and will go unnoticed (a). Only very subtle differences exist between the corrupted and the true tractogram (a,d). Correcting for column order permutations and column sign changes comes to no avail (b). The full 3D correction on the other hand results in a correctly restored DEC FA map and whole brain tractogram (c-d).
(a) Before correction
(b) After column permutation and sign correction
(c) After full correction
(d) Original
Figure 12: Example of the gradient correction on in-vivo data set D with a gradient table subject to a subtle 3D rotation. FA, DEC and a whole brain tractogram are shown before correction (a), after gradient column permutation and sign correction (b), after the full rigid correction (c) and for the original data set (d).
20
4. Discussion & conclusion For rotationally variant diffusion parameters to be meaningful, the diffusion gradient orientations should use the same coordinate frame as the voxels in the DW data set. The meta data generated by MRI scanners should provide all the necessary information to properly align the gradient table with the subject coordinate frame. In practice, however, this procedure can be error-prone and small angulation errors can easily go unnoticed. As such, wrongly oriented gradient tables are an often overlooked and hard to detect pitfall in the processing of diffusion MRI (Jones and Cercignani, 2010; Leemans and Jones, 2009). In this work, we have introduced a method that allows rigid registration of the DW gradient table to the corresponding DW images, using a metric based on whole brain DTI tractography. Our tool enables researchers to check, and, if necessary, correct the gradient tables used in their study in a fully automated way. The proposed method uses whole brain fiber tractography as a measure of goodness-of-fit, which means it takes into account the global, long range, alignment patterns of the local diffusion tensors, something which cannot be achieved with methods that simply check how well the tensors are aligned with their neighboring tensors. Indeed, rotating a set of aligned tensors, will keep them aligned, only with respect to a different local orientation. To enable a fast metric calculation, we have combined several algorithmic optimizations from the literature, that help speeding up the process of whole brain tractography. Despite these optimizations, repeated metric calculations can still require a significant amount time when many seed points are used. Since a reliable evaluation of the metric requires a large number of seed points, repeated evaluation will still have a significant computational cost. In addition, even with a large number of seed points, the metric still exhibits spurious maxima, due to the discrete nature of tractography, which would cause deterministic optimization routines to report suboptimal solutions. In order to keep the computation time of the metric to a minimum and allow accurate estimation of the rotation parameters, we chose to use the stochastic gradient ascent. The stochastic gradient ascent method follows the same scheme as the deterministic gradient ascent, with the distinction that the derivative of the cost function is replaced by a fast stochastic approximation. By using a new, randomly selected subset of seed points every iteration of the optimization process, a bias in the approximation error is avoided. 21
Indeed, our simulations show that, by adopting a stochastic optimization approach, we are able to obtain corrected gradient orientations with sub degree accuracy and precision for SNR values down to 4, which is much lower than the SNR reported in the bulk of clinical studies. We have also shown that the proposed method works on a range of real clinical data sets without requiring the user to change any parameters of the metric function or the optimization routine. While changing the default number of random seed points of 200 to a lower value will reduce the metric computation time proportionally, it will potentially result in the metric becoming so unreliable that the stochastic gradient descent fails. Also, while reducing the default number of gradient descent iterations of 150 will reduce the computation time, it is not advised to do so as this could prevent the gradient descent from converging to the optimum. An additional benefit of the stochastic gradient descent is that it can potentially avoid getting converging to local minima due to the presence of local anomalies in the fiber bundles (e.g. injured or missing tracts). As the stochastic gradient descent selects a small, different set of random seed points at each iteration. This way, a different random part of the brain is sampled at each iteration, avoiding getting stuck in the local optima that could result from local anomalies, should the complete brain be sampled. Using our current Matlab routines, the column permutation and sign correction took on average 23 seconds, while the full correction of 3D rotations required on average 27 minutes (3.5 GHz quad core desktop processor). Note that the current algorithm was written in Matlab, which is not the ideal platform to perform an iterative algorithm such as fiber tracking, and further speed ups could be obtained by moving to C or C++. Note that our current method uses DTI tractography, which is sub optimal in the presence of crossing fibers (Tuch et al., 2002). While using a tractography method that resolves crossing fibers will clearly result in different metric values, we expect the optima to be found at the same locations. However, DTI tractography can generally be performed much faster than any of the HARDI alternatives and is compatible with a large set of acquisition schemes (simple DTI schemes (Jones et al., 1999), HARDI schemes (Tuch et al., 2002), multi-shell schemes such as those used for diffusion kurtosis imaging (Jensen et al., 2005), and even full blown Cartesian-sampling (Wedeen et al., 2005)). So despite the fact that our approach uses DTI tractography, the recovered gradient orientations can be applied in general DWI post-processing. 22
In this work, we did not directly investigate the One of the benefits of the stochastic gradient descent as an optimizer, is that is robust against such local anomalies. The stochastic gradient descent works by selecting a small, different set of random seed points at each iteration. This way, a different random part of the brain is sampled at each iteration, avoiding getting stuck in the local optima that could result from local anomalies, should the complete brain be sampled. Acknowledgements This work was supported by the Interuniversity Attraction Poles Program (P7/11) References Arsigny, V., Fillard, P., Pennec, X., Ayache, N., 2006. Log-Euclidean metrics for fast and simple calculus on diffusion tensors. Magnetic Resonance in Medicine 56, 411–421. Assaf, Y., Pasternak, O., 2008. Diffusion tensor imaging (DTI)-based white matter mapping in brain research: a review. Journal of Molecular Neuroscience 34, 51–61. Basser, P.J., 1995. Inferring microstructural features and the physiological state of tissues from diffusion-weighted images. NMR in Biomedicine 8, 333–344. ¨ Basser, P.J., Ozarslan, E., 2010. Anisotropic diffusion: from the apparent diffusion coefficient to the apparent diffusion tensor, in: Jones, D.K. (Ed.), Diffusion MRI: theory, methods, and applications. Oxford University Press. Basser, P.J., Pajevic, S., 2000. Statistical artifacts in diffusion tensor MRI (DT-MRI) caused by background noise. Magnetic Resonance in Medicine 44, 41–50. Basser, P.J., Pajevic, S., Pierpaoli, C., Duda, J., Aldroubi, A., 2000. In vivo fiber tractography using DT-MRI data. Magnetic Resonance in Medicine 44, 625–632. 23
Basser, P.J., Pierpaoli, C., 1996. Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI. Journal of Magnetic Resonance Series B 111, 209–219. Conturo, T.E., Lori, N.F., Cull, T.S., Akbudak, E., Snyder, A.Z., Shimony, J.S., McKinstry, R.C., Burton, H., Raichle, M.E., 1999. Tracking neuronal fiber pathways in the living human brain. Proceedings of the National Academy of Sciences of the United States of America 96, 10422–10427. Gudbjartsson, H., Patz, S., 1995. The Rician distribution of noisy MRI data. Magnetic Resonance in Medicine 34, 910–914. Hasan, K.M., Basser, P.J., Parker, D.L., Alexander, A.L., 2001. Analytical computation of the eigenvalues and eigenvectors in DT-MRI. Journal of Magnetic Resonance 152, 41–47. Jenkinson, M., Beckmann, C.F., Behrens, T.E., Woolrich, M.W., Smith, S.M., 2012. Fsl. NeuroImage 62, 782–790. Jensen, J.H., Helpern, J.A., Ramani, A., Lu, H., Kaczynski, K., 2005. Diffusional kurtosis imaging: the quantification of non-gaussian water diffusion by means of magnetic resonance imaging. Magn Reson Med 53, 1432–1440. Jeurissen, B., Leemans, A., Jones, D.K., Tournier, J.D., Sijbers, J., 2013. Investigating the prevalence of complex fiber configurations. Human Brain Mapping 34, 2747–2766. Johansen-Berg, H., Behrens, T.E.J., 2006. Just pretty pictures? What diffusion tractography can add in clinical neuroscience. Current Opinion in Neurology 19, 379–385. Jones, D.K., 2003. Determining and visualizing uncertainty in estimates of fiber orientation from diffusion tensor MRI. Magnetic Resonance in Medicine 49, 7–12. Jones, D.K., Cercignani, M., 2010. Twenty-five pitfalls in the analysis of diffusion MRI data. NMR in Biomedicine 23, 803–820. Jones, D.K., Horsfield, M.A., Simmons, A., 1999. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magnetic Resonance in Medicine 42, 515–525. 24
Klein, S., Staring, M., Pluim, J.P., 2007. Evaluation of optimization methods for nonrigid medical image registration using mutual information and Bsplines. IEEE Trans Image Process 16, 2879–2890. Leemans, A., Jones, D.K., 2009. The B-matrix must be rotated when correcting for subject motion in DTI data. Magnetic Resonance in Medicine 61, 1336–1349. Leemans, A., Sijbers, J., Verhoye, M., Van der Linden, A., Van Dyck, D., 2005. Mathematical framework for simulating diffusion tensor MR neural fiber bundles. Magnetic Resonance in Medicine 53, 944–953. Mori, S., Crain, B.J., Chacko, V.P., Van Zijl, P., 1999. Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Annals of Neurology 45, 265–269. Mori, S., Zhang, J., 2006. Principles of diffusion tensor imaging and its applications to basic neuroscience research. Neuron 51, 527–539. Moseley, M.E., Cohen, Y., Kucharczyk, J., Mintorovitch, J., Asgari, H.S., Wendland, M.F., Tsuruda, J., Norman, D., 1990. Diffusion-weighted MR imaging of anisotropic water diffusion in cat central nervous system. Radiology 176, 439. Pajevic, S., Pierpaoli, C., 1999. Color schemes to represent the orientation of anisotropic tissues from diffusion tensor data: application to white matter fiber tract mapping in the human brain. Magnetic Resonance in Medicine 42, 526–540. Robbins, H., Monro, S., 1951. A stochastic approximation method. The Annals of Mathematical Statistics , 400–407. Stejskal, E.O., Tanner, J.E., 1965. Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. The Journal of Chemical Physics 42, 288. Tax, C., Jeurissen, B., Vos, S.B., Viergever, M.A., Leemans, A., in press. Recursive calibration of the fiber response function for spherical deconvolution of diffusion mri data. NeuroImage .
25
Tournier, J.D., Calamante, F., Connelly, A., 2012. Mrtrix: Diffusion tractography in crossing fiber regions. International Journal of Imaging Systems and Technology 22, 53–66. Tournier, J.D., Mori, S., Leemans, A., 2011. Diffusion tensor imaging and beyond. Magnetic Resonance in Medicine 65, 1532–1556. Tuch, D.S., Reese, T.G., Wiegell, M.R., Makris, N., Belliveau, J.W., Wedeen, V.J., 2002. High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity. Magnetic Resonance in Medicine 48, 577–582. Wedeen, V.J., Hagmann, P., Tseng, W., Reese, T.G., Weisskoff, R.M., 2005. Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging. Magnetic Resonance in Medicine 54, 1377–1386.
26
Appendix
(a) data set A
(b) data set B
(c) data set C
(d) data set D original before correction column permutation and sign correction full correction
Figure 1: Gradient correction on in-vivo data sets with a gradient table subject to: (a) a column permutation (data set A); (b) a column sign change (data set B); (c) a large, random 3D rotation (data set C); and (d) a subtle 3D rotation (data set D). True gradient orientations are displayed as black circles; corrupted gradients as red squares; gradient orientations corrected by column permutation and sign correction as blue plus signs; and gradient orientations corrected by the full rigid correction as green crosses.
27
(IGHLIGHTS
· · · · ·
Our method aligns the diffusion gradient orientations to the corresponding data. We optimize a metric based on whole brain fiber tractography length. Stochastic optimization is used to ensure fast and robust realignment. The method is fully automated, data driven and independent of scanner meta-data. It can be applied to data sets with different acquisition settings and properties.
'RAPHICAL !BSTRACT
S1361-8415(14)00093-0 http://dx.doi.org/10.1016/j.media.2014.05.012 MEDIMA 897
To appear in:
Medical Image Analysis
Received Date: Revised Date: Accepted Date:
9 December 2013 20 May 2014 27 May 2014
Please cite this article as: Jeurissen, B., Leemans, A., Sijbers, J., Automated correction of improperly rotated diffusion gradient orientations in diffusion weighted MRI, Medical Image Analysis (2014), doi: http://dx.doi.org/ 10.1016/j.media.2014.05.012
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Automated correction of improperly rotated diffusion gradient orientations in diffusion weighted MRI Ben Jeurissena,∗, Alexander Leemansb , Jan Sijbersa a
b
iMinds-Vision Lab, Department of Physics, University of Antwerp, Belgium Image Sciences Institute, University Medical Center Utrecht, The Netherlands
Abstract Ensuring one is using the correct gradient orientations in a diffusion MRI study can be a challenging task. As different scanners, file formats and processing tools use different coordinate frame conventions, in practice, users can end up with improperly oriented gradient orientations. Using such wrongly oriented gradient orientations for subsequent diffusion parameter estimation will invalidate all rotationally variant parameters and fiber tractography results. While large misalignments can be detected by visual inspection, small rotations of the gradient table (e.g. due to angulation of the acquisition plane), are much more difficult to detect. In this work, we propose an automated method to align the coordinate frame of the gradient orientations with that of the corresponding diffusion weighted images, using a metric based on whole brain fiber tractography. By transforming the gradient table and measuring the average fiber trajectory length, we search for the transformation that results in the best global ‘connectivity’. To ensure a fast calculation of the metric we included a range of algorithmic optimizations in our tractography routine. To make the optimization routine robust to spurious local maxima, we use a stochastic optimization routine that selects a random set of seed points on each evaluation. Using simulations, we show Corresponding author Email: [email protected] Address: Vision Lab University of Antwerp (CDE) Universiteitsplein 1, N.1.18 B-2610 Wilrijk, Belgium Phone: +32 3 265 24 77 Fax: +32 3 265 22 45 ∗
Preprint submitted to Medical Image Analysis
May 31, 2014
that our method can recover the correct gradient orientations with high accuracy and precision. In addition, we demonstrate that our technique can successfully recover rotated gradient tables on a wide range of clinically realistic data sets. As such, our method provides a practical and robust solution to an often overlooked pitfall in the processing of diffusion MRI. Keywords: diffusion-weighted MRI, coordinate frame, rotationally variant parameter, fiber tractography, diffusion gradient reorientation, stochastic approximation 1. Introduction Diffusion-weighted (DW) MRI is a non-invasive imaging technique that is sensitive to the random microscopic motion of water molecules (Stejskal and Tanner, 1965). In the brain white matter (WM), which consists of tightly packed bundles of neuronal axons, the regular arrangement of fibers introduces a directional dependence of this motion (Moseley et al., 1990). The degree of this directionality provides an indirect measure for the organization of the underlying tissue, while the preferred diffusion orientation provides an estimate of the local fiber orientation of the bundles (Basser, 1995). These local fiber orientation estimates can be pieced together to infer the long-range pathways connecting distant regions of the brain, using a technique called WM fiber tractography (Conturo et al., 1999; Mori et al., 1999). DW-MRI is currently the only imaging technique that allows assessment of WM organization and connectivity in vivo and it is increasingly used in clinical studies to investigate brain WM, and its development and disorders (Johansen-Berg and Behrens, 2006; Mori and Zhang, 2006; Assaf and Pasternak, 2008). Diffusion metrics extracted from DW images can be categorized into two categories based on their rotational (in)variance. Rotationally invariant diffusion metrics embody intrinsic features of the local diffusion pattern and are independent of the orientation of the anisotropic structure, the orientation of the specimen with respect to the laboratory coordinate system, and ¨ even the choice of the laboratory coordinate system (Basser and Ozarslan, 2010). Notable examples are the scalar diffusion parameters obtained from diffusion tensor imaging (DTI) such as the widely used fractional anisotropy (FA) and mean diffusivity (MD) (Basser and Pierpaoli, 1996). Both parameters are extracted from the eigenvalues of the diffusion tensor, making them rotationally invariant by definition. 2
Rotationally variant diffusion parameters, on the other hand, embody orientational features of the local diffusion pattern and are dependent of the orientation of the anisotropic structure, the orientation of the specimen with respect to the laboratory coordinate system, and the choice of the laboratory coordinate system. Notable examples are the principal diffusion vector (PDV) estimated by DTI, diffusion/fiber orientation distribution functions obtained from high angular resolution diffusion imaging (HARDI) (Tournier et al., 2011) and fiber tractography results. For rotationally variant parameters to be meaningful, it is important that the coordinate system of the diffusion sensitizing gradient orientations and the spatial coordinate system of the voxels in the diffusion weighted data set coincide. In theory, the meta data generated by MRI scanner software should provide all the necessary information to properly align the gradient orientations with the subject coordinate frame. In practice, however, this procedure is error-prone, as the appropriate steps are dependent on the coordinate system conventions of the scanner software, the file format and even the post-processing tool. For example, while both MRtrix (Tournier et al., 2012) and FSL (Jenkinson et al., 2012) support DWI data in the NifTI format and while both represent the gradient vectors as simple text files, MRtrix assumes that the gradient vectors are provided with respect to real-world/scanner coordinates, whereas FSL assumes that the gradient vectors are provided with respect to the image axes. While large misalignments (e.g. flipping or permutations of the gradient orientation coordinate axes) can be detected by visual inspection of PDV maps or tractography results, small rotations of the gradient table (e.g. due to angulation of the acquisition plane), are much more difficult to detect. Nevertheless, even small rotations have a large impact on tractography results due to an accumulation of errors. In this work, we propose a data-driven method that does not rely on header information to automatically detect and correct for the rigid misalignment of the gradient table and the DW data, using a metric based on whole brain tractography.
3
2. Methods The use of incorrect gradient orientations will result in the premature termination of fiber trajectories, which can be measured using whole brain fiber tractography. By transforming the gradient table and measuring the average fiber trajectory length, one can search for the transformation that results in the best global ‘connectivity’. Finding the correct DW gradient settings then becomes akin to a registration problem (Figure 1). The gradient orientations (moving entity) are subjected to a 3D rotation (transformation), such that they are aligned with the DW images (reference entity). Reference entity: DW images
Metric: Average tract length
Moving entity: Gradient table
iterate
Transform: 3D rotation
fitness value
Optimizer: Exhaustive search / Stochastic gradient ascent
transform parameters
Figure 1: Schematic overview of the gradient table reorientation method.
2.1. Transform The 3D rotation R of the gradient table can be parameterized as successive rotations around the main orthogonal axes: R(θ) = Rx (θx )Ry (θy )Rz (θz )
(1)
resulting in a parameter vector θ = (θx , θy , θz ) with 3 transformation parameters. Note that, since the gradient table is rotated as a whole, this is a rigid body transformation. Also note that translations do not have to be taken into account, as they do not affect the gradient orientations.
4
2.2. Metric The metric, which provides a measure of how well the gradient orientations match the DW images, is defined as the average trajectory length f (θ, n), measured with whole-brain streamline DTI tractography, where θ parametrizes the 3D rotation of the gradient table and n is the number of seed points, evenly distributed over the whole brain. As fiber tractography is more reliable in structures with high FA, the length contribution of each step is weighted by the local FA value. To ensure robust fiber tractography even on data of low quality, the DWI images are smoothed with a 3D Gaussian filter prior to tractography. The metric exploits the fact that the brain is a complicated 3D network of fiber bundles. While rotating the gradient orientations away from the true orientations may result in a few slightly longer false trajectories in some locations, it will at the same time reduce the length of many true trajectories. Since finding the optimal gradient table requires many repeated evaluations of the metric, it is imperative that the computation of f (θ, n) is fast. However, as whole-brain tractography is relatively time-consuming, several optimizations were implemented to speed up the standard streamline DTI tractography algorithm. Traditionally, DTI tractography is performed as follows (Basser et al., 2000). First, the DW signal at the current position of the trajectory is obtained using trilinear interpolation. Next, the diffusion tensor is estimated from the interpolated signal. Then an eigenvalue decomposition is performed on the diffusion tensor and the eigenvector corresponding to the largest eigenvalue is assumed to be the current fiber orientation. Finally, the trajectory is advanced by a fixed step size along this orientation and the whole process is repeated. The most time consuming operations in this algorithm are the interpolation, the diffusion tensor fit and the eigenvalue decomposition. The computational complexity of the DTI tractography algorithm can be reduced significantly using the Log Euclidian framework (Arsigny et al., 2006), as follows. Prior to tractography, the log diffusion tensor field is calculated for the whole brain. Then, the log diffusion tensor at the current position of the trajectory can be obtained using trilinear interpolation directly on the log diffusion tensor field. Next, an eigenvalue decomposition is performed on the diffusion tensor and the eigenvector corresponding to the largest eigenvalue is assumed to be the current fiber orientation. Finally, the trajectory is advanced by a fixed step size along this orientation and the whole process is repeated. Following this approach, only a single diffusion 5
tensor fit is required and the number of interpolations at each tractography step is limited to 6 as compared to the number of diffusion weighted images for the traditional method, greatly reducing computation time and memory consumption. Instead of rotating the gradient table and recalculating the log diffusion tensor at each metric evaluation, our algorithm directly rotates the log tensor, requiring only a single diffusion tensor fit for repeated evaluations of the metric. To avoid unnecessary calculations outside the brain, a rough skull stripping was performed on the b = 0 s/mm2 images, separating brain from background. The only remaining bottleneck is the eigenvalue decomposition which needs to be performed at each tractography step. To reduce the computational load of this operation, the analytical eigenvalue decomposition proposed in (Hasan et al., 2001) was used, which exploits the fact that the diffusion tensor is a 3×3 symmetric and semi-positive-definite matrix, and which is approximately 40 times faster than standard iterative decomposition techniques. Finally, to reduce computation time even further, whole brain tractography is not performed from all N voxels of the whole brain, but from a set of n < N seed points distributed over the whole brain, allowing a trade-off between metric calculation speed and metric precision. For reasons outlined below, we also define a stochastic cost function (Robbins and Monro, 1951) f˜(θ, n) which is the approximate average tract length resulting from whole brain tractography with rotation angles θ and n the number of seed points randomly distributed over the whole brain according to a uniform distribution. As a consequence, each time f˜(θ, n) is called, it probes a different set of seed points, resulting in a different metric outcome. 2.3. Optimization We provide two optimization approaches, depending on what type of gradient corruptions are present in the data. The first approach assumes that the gradient table is only corrupted by means of column order permutations and/or column sign changes, a situation that often occurs in practice, due to different coordinate system conventions. The number of possible column order permutations is 3! = 6, while the number of possible sign changes is 4 (sign change in the first column, the second column, the third column and no sign change), resulting in 6 × 4 = 24 possible solutions in total. This optimization problem can be solved in a 6
brute force fashion by exhaustively evaluating the metric f (θ, n). From this list of evaluations fi (with i = 1, ..., 24) we select the permutation and/or sign change that results in the longest average fiber length. To avoid selecting false solutions in the case of insufficient signal-to-noise ratio (SNR), instead of just evaluating f (θ, n), we repeatedly evaluate f˜(θ, n) and the solution is only deemed optimal if it is consistently selected as having the longest average fiber orientation across the stochastic function evaluations. If the same solution is not found consistently across stochastic function evaluations, the method reports that the solution is undecided due to insufficient SNR. In that case, the method is rerun with twice the amount of seed points n, increasing the metric precision. This way, the method dynamically adjusts the number of seed points required for precise metric evaluation to the quality of the data and avoids false solutions. The second approach assumes that the gradient table is corrupted by means of a general 3D transformation. Note that, due to the antipodal symmetry of the diffusion signal, this approach also resolves the column order permutations and column sign changes resolved with the first approach. Since a reliable evaluation of the metric requires a large number of seed points n, repeated evaluation has a very high computational cost. In addition, even with high n, the metric still exhibits spurious maxima, which would cause deterministic optimization routines to report suboptimal solutions (Figure 2). In order to keep the computation time of the metric to a minimum and allow accurate estimation of the rotation parameters, we chose to use the stochastic gradient ascent (Robbins and Monro, 1951) which has previously been used successfully for the purpose of image registration (Klein et al., 2007). The stochastic gradient ascent method follows the same scheme as the deterministic gradient ascent, with the distinction that the derivative of the cost function is replaced by a fast stochastic approximation. By using a new, randomly selected subset of seed points every iteration of the optimization process, a bias in the approximation error is avoided. More specifically, in order to find the optimal set of rotation angles θ, a stochastic gradient ascent is defined as follows: θk+1 = θk + γk ∇f˜(θk , n)
with
k = 0, ..., K
(2)
with γk the scalar gain factor or “step size” and f˜(θk , n) the approximate average tract length resulting from whole brain tractography with rotation 7
f
30
n = 50 n = 500 n = 5000
20
10 −180
−120
0
−60
60
120
180
θx ( ) (a) f over the full domain of θx n = 5000
f
28
27.5
27 −5
0
−2.5
2.5
5
θx ( ) (b) f near the optimum Figure 2: Cost function values f for a phantom data set with correctly oriented gradient table and SNR = 10 (defined on the b = 0 s/mm2 images) as a function of θx for selected values of n (a). A zoomed-in plot near the optimum is provided in (b).
angles θk and n randomly distributed seed points. The approximated gradiˆ in contrast to the exact ent does not necessarily vanish close to the solution θ, ˆ = 0. Thus, convergence must be forced by enderivative that satisfies ∇f (θ) suring γk → 0 as k → K. In this work, K = 150 and γk ≡ γ(k) = 100e−0.068k . Using the stochastic gradient ascent makes the optimization routine robust to spurious maxima, and also allows the number of seed points n to be small, significantly speeding up the metric evaluation (Klein et al., 2007). Unless stated otherwise, n = 200 seed points randomly distributed over the whole brain, were used in all subsequent experiments. While a small number of 8
seed points results in f˜(θk , n) having low precision, the repeated evaluations of f˜(θk , n) will converge to the true optimum. Since the gradient ascent can only find the local minimum, we use a multi-start approach starting from 30 points distributed evenly in the space of 3D rotations. 2.4. Simulations To evaluate the proposed method quantitatively, a noiseless full-brain DWI data set consisting of 6 b = 0 s/mm2 and 60 b = 1200 s/mm2 images and with a voxel size of 2.4 × 2.4 × 2.4 mm3 was simulated (Leemans et al., 2005). From this data set, noisy Rician distributed data sets were obtained for different noise levels, with SNR in the b = 0 s/mm2 image ranging from 1 to 30 and 100 noisy realizations for each noise level (Gudbjartsson and Patz, 1995). 2.4.1. Gradient column order permutations and sign changes For each noisy realization, we checked whether the method could distinguish the correct column order permutation and sign changes from the wrong ones. We also investigated the amount of seed points required in the metric to make this distinction. 2.4.2. General 3D rotations For each noisy realization, the gradient table was corrupted by a random 3D rotation. After correction with the proposed method, quality of the corrected gradient table was assessed with respect to accuracy and precision as follows. For each gradient orientation, an average gradient orientation was obtained from the 100 corrected orientations by extracting the first eigenvector of the mean dyadic tensor (Basser and Pajevic, 2000) and the minimum angles subtended between each estimated gradient orientation and the average gradient orientation were recorded. A measure of accuracy was obtained by calculating the angular bias of the average gradient orientation compared to the ground truth gradient orientation (Basser and Pajevic, 2000). A measure of precision was obtained by calculating the 95% confidence interval (CI) of the one-tailed distribution of minimum angles between the average gradient orientation and the individual gradient orientations (Jones, 2003).
9
2.5. Experiments with in vivo human brain data The method was also tested on a set of real human brain DWI data sets with different characteristics: • Data set A is a DTI data set acquired on a 3T Siemens Trio scanner, using a 12-channel receiver head coil. Diffusion weighting was applied along 30 gradient directions with b = 1000 s/mm2 . Additionally, 2 nonDW images were acquired. Other imaging parameters were: TR/TE: 8100/116 ms; voxel size: 2.5 × 2.5 × 2.5 mm3 ; matrix: 96 × 96; slices: 54; NEX: 1. • Data set B is a HARDI data set acquired on a 3T GE HDx Signa scanner, using an 8-channel receiver head coil and cardiac gating (Jeurissen et al., 2013). Diffusion weighting was applied along 60 gradient directions with b = 1200 s/mm2 . Additionally, 6 non-DW images were acquired. Other imaging parameters were: TR: 20 R-R intervals; TE: 86 ms; voxel size: 2.4 × 2.4 × 2.4 mm3 ; matrix: 96 × 96; slices: 60; NEX: 1. • Data set C is a HARDI data set with high spatial resolution acquired on a 3T Philips Achieva scanner, using an 8-channel receiver head coil (Tax et al., in press). Diffusion weighting was applied along 60 gradient directions with b = 2500 s/mm2 . Additionally, 1 non-DW image was acquired. Other imaging parameters were: TR/TE: 10265/107 ms; voxel size: 2 × 2 × 2 mm3 ; matrix: 112 × 112; slices: 70; NEX: 1 . • Data set D is a multi-shell DW data set acquired with high spatial resolution on a 3T GE Discovery MR750 scanner, using an 8-channel receiver head coil. Diffusion weightings of b = 0, 700, 1200 and 2800 s/mm2 were applied in 5, 25, 45 and 75 directions, respectively. Other imaging parameters were: TR/TE: 9500/100 ms; voxel size: 2 × 2 × 2 mm3 ; matrix: 120 × 120; slices: 68; NEX: 1. We deliberately corrupted the gradient tables of these data sets using a random column permutation, a random column sign change, a small random 3D rotation and a large random 3D rotation. Subsequently, the corrupted data sets were processed using our correction method, to test whether it could correctly restore the gradient orientations.
10
3. Results 3.1. Simulations Image quality for the different simulated SNR levels can be appreciated from Figure 3. Note that SNR is defined on the b = 0 s/mm2 image and that for realistic clinical studies it is in the range of [10-30]. However, to ensure proper functioning of our method even on low quality data, we also simulated SNR values well below 10.
DWI
DEC
SNR
30
20
10
5
2
Figure 3: Illustration of the image quality at different simulated SNR levels, for a single DW image and for the resulting DEC FA map. Note that SNR is defined on the b = 0 s/mm2
11
3.1.1. Gradient column order permutations and sign changes Figure 4 shows the results of the exhaustive evaluation of the different combinations of gradient table column order permutations (rows) and sign changes (columns), for a single realization of different SNR levels. The correct answer is the top left cell, which corresponds to no permutation and no sign changes. Note that one can clearly distinguish the correct solution up till SNR = 5. However, at SNR = 2, the method is still able to repeatedly recognize the first solution as the optimum, over repeated stochastic evaluations of the metric f . At SNR = 1, no solution stands out from the rest and the method reports no answer. Figure 5 shows the percentage of correct, no, and wrong solutions as a function of SNR. Note that the method achieves a 100% success rate all the way up to SNR = 2. Starting from SNR = 1, the method reports that it can no longer reliably estimate the correct gradient orientations in all cases. Note that the method makes no false claims about the correct gradient orientations, even at low SNR. Figure 6 shows the number of seed points required by the method to restore the gradient table as a function of SNR. Note that 50 seed points were sufficient all the way down to SNR = 4. Dropping below SNR = 4, the method gradually requires more seed points. At SNR = 2, the method requires at most 200 seed points, which is still reasonable.
12
column permutations
column permutations
column permutations column sign changes
column sign changes
column permutations
(c) SNR = 10
column sign changes
column permutations
(b) SNR = 20
column permutations
(a) SNR = 30
column sign changes
column sign changes
(d) SNR = 5
column sign changes
(e) SNR = 2
5
10
15
20
(f) SNR = 1
25
30
Figure 4: Average metric values (over repeated stochastic metric evaluations) for the different combinations of gradient table column order permutations and sign changes, for different SNRs (single realization). The correct answer is the top left cell, which corresponds to no permutation and no sign changes.
13
100 correct solution no solution wrong solution
80
%
60 40 20 0
30
20
10
5
4
3
2
1
SNR Figure 5: Percentage of correct (green), no (blue) and wrong (red) solutions as a function of SNR, over 200 noise realizations.
100 80 50 100 200 400 800 1600 3200 6400
%
60 40 20 0
30
20
10
5
4
3
2
1
SNR Figure 6: The number of seed points required by the method to restore the gradient table as a function of SNR, over 200 noise realizations.
14
3.1.2. General 3D rotations Figure 7 shows the gradient orientations recovered with our method after random 3D rotations for 100 noisy realizations of the data set, together with the true gradient orientations. For SNR values all the way down to 5, the corrected gradient orientations are tightly clustered around the true gradient orientations, with no discernible bias. For SNR = 2, precision of the corrected gradient orientations is markedly decreased, but the corrected gradient orientations are still clustered around the true orientations, suggesting only a small bias. Note that an SNR as low as 2 in the b = 0 s/mm2 images, is extremely low and does no longer represent realistic clinical data.
(a) SNR = 30
(b) SNR = 20
(c) SNR = 10
(d) SNR = 5
(e) SNR = 2
(f) SNR = 1
Figure 7: The recovered gradient orientations (colored dots) vs the true gradient orientations (black crosses) for 100 noisy realizations of the data set. Corresponding corrected gradient orientations are given the same color to aid visual clustering.
Figure 8 quantifies both accuracy and precision as a function of SNR. Subdegree accuracy is obtained starting from SNR = 2. Sub-degree precision is obtained, starting from SNR = 4. This is well within the boundaries of clinically realistic SNR. 15
bias ( )
10 1
0.1
12345
10
20
30
20
30
SNR (a) Accuracy
95% CI ( )
100
10
1
0.1
12345
10 SNR (b) Precision
Figure 8: Accuracy (a) and precision (b) of the recovered gradient tables as a function of SNR. The dashed red lines indicate the points at which sub-degree accuracy (a) and sub-degree precision (b) are achieved.
16
3.2. Experiments with in vivo human brain data 3.2.1. Data set A Figure 9 shows an FA map, a DEC (directionally encoded color) FA map (Pajevic and Pierpaoli, 1999) and a whole brain tractogram for data set A, after corruption with a column order permutation (a), after correcting for column order permutations and column sign changes (b), after correction of a full 3D rotation (c), and without corruption (d). Note that column order permutation can easily be detected by visual inspection of the DEC FA map and causes a significant disruption of the tracts in the whole brain tractogram (a). Both our correction strategies correctly restore the DEC FA map and whole brain tractogram (b-d).
(a) Before correction
(b) After column permutation and sign correction
(c) After full correction
(d) Original
Figure 9: Example of the gradient correction on in-vivo data set A with a gradient table column permutation. FA, DEC and a whole brain tractogram are shown before correction (a), after gradient column permutation and sign correction (b), after the full correction (c) and for the original data set (d).
17
3.2.2. Data set B Figure 10 shows an FA map, a DEC FA map and a whole brain tractogram for data set B, after corruption with a column sign change (a), after correcting for column order permutations and column sign changes (b), after correction of a full 3D rotation (c), and without corruption (d). Note that column sign change can no longer be detected by visual inspection of the DEC FA map (a). However, the column sign change still causes a significant disruption of the tracts in the whole brain tractogram (a). Both our correction strategies correctly restore the DEC FA map and whole brain tractogram (b-d).
(a) Before correction
(b) After column permutation and sign correction
(c) After full correction
(d) Original
Figure 10: Example of the gradient correction on in-vivo data set B with a gradient table column sign change. FA, DEC and a whole brain tractogram are shown before correction (a), after gradient column permutation and sign correction (b), after the full correction (c) and for the original data set (d).
18
3.2.3. Data set C Figure 11 shows an FA map, a DEC FA map and a whole brain tractogram for data set C, after corruption with a large, random 3D rotation (a), after correcting for column order permutations and column sign changes (b), after correction of a full 3D rotation (c), and without corruption (d). Note that the large rotation can detected by visual inspection of both the DEC FA map and the whole brain tractogram (a). However, manual realignment is now no longer feasible. Correcting for column order permutations and column sign changes only marginally improves the DEC FA map and the whole brain tractogram, and cannot completely restore it (b). The full 3D correction on the other hand results in a correctly restored DEC FA map and whole brain tractogram (c-d).
(a) Before correction
(b) After column permutation and sign correction
(c) After full correction
(d) Original
Figure 11: Example of the gradient correction on in-vivo data set C with a gradient table subject to a large 3D rotation. FA, DEC and a whole brain tractogram are shown before correction (a), after gradient column permutation and sign correction (b), after the full correction (c) and for the original data set (d).
19
3.2.4. Data set D Figure 12 shows an FA map, a DEC FA map and a whole brain tractogram for data set D, after corruption with a subtle, random 3D rotation (a), after correcting for column order permutations and column sign changes (b), after correction of a full 3D rotation (c), and without corruption (d). Note that the small rotation is virtually impossible to detect by visual inspection on both the DEC FA map and the whole brain tractogram and will go unnoticed (a). Only very subtle differences exist between the corrupted and the true tractogram (a,d). Correcting for column order permutations and column sign changes comes to no avail (b). The full 3D correction on the other hand results in a correctly restored DEC FA map and whole brain tractogram (c-d).
(a) Before correction
(b) After column permutation and sign correction
(c) After full correction
(d) Original
Figure 12: Example of the gradient correction on in-vivo data set D with a gradient table subject to a subtle 3D rotation. FA, DEC and a whole brain tractogram are shown before correction (a), after gradient column permutation and sign correction (b), after the full rigid correction (c) and for the original data set (d).
20
4. Discussion & conclusion For rotationally variant diffusion parameters to be meaningful, the diffusion gradient orientations should use the same coordinate frame as the voxels in the DW data set. The meta data generated by MRI scanners should provide all the necessary information to properly align the gradient table with the subject coordinate frame. In practice, however, this procedure can be error-prone and small angulation errors can easily go unnoticed. As such, wrongly oriented gradient tables are an often overlooked and hard to detect pitfall in the processing of diffusion MRI (Jones and Cercignani, 2010; Leemans and Jones, 2009). In this work, we have introduced a method that allows rigid registration of the DW gradient table to the corresponding DW images, using a metric based on whole brain DTI tractography. Our tool enables researchers to check, and, if necessary, correct the gradient tables used in their study in a fully automated way. The proposed method uses whole brain fiber tractography as a measure of goodness-of-fit, which means it takes into account the global, long range, alignment patterns of the local diffusion tensors, something which cannot be achieved with methods that simply check how well the tensors are aligned with their neighboring tensors. Indeed, rotating a set of aligned tensors, will keep them aligned, only with respect to a different local orientation. To enable a fast metric calculation, we have combined several algorithmic optimizations from the literature, that help speeding up the process of whole brain tractography. Despite these optimizations, repeated metric calculations can still require a significant amount time when many seed points are used. Since a reliable evaluation of the metric requires a large number of seed points, repeated evaluation will still have a significant computational cost. In addition, even with a large number of seed points, the metric still exhibits spurious maxima, due to the discrete nature of tractography, which would cause deterministic optimization routines to report suboptimal solutions. In order to keep the computation time of the metric to a minimum and allow accurate estimation of the rotation parameters, we chose to use the stochastic gradient ascent. The stochastic gradient ascent method follows the same scheme as the deterministic gradient ascent, with the distinction that the derivative of the cost function is replaced by a fast stochastic approximation. By using a new, randomly selected subset of seed points every iteration of the optimization process, a bias in the approximation error is avoided. 21
Indeed, our simulations show that, by adopting a stochastic optimization approach, we are able to obtain corrected gradient orientations with sub degree accuracy and precision for SNR values down to 4, which is much lower than the SNR reported in the bulk of clinical studies. We have also shown that the proposed method works on a range of real clinical data sets without requiring the user to change any parameters of the metric function or the optimization routine. While changing the default number of random seed points of 200 to a lower value will reduce the metric computation time proportionally, it will potentially result in the metric becoming so unreliable that the stochastic gradient descent fails. Also, while reducing the default number of gradient descent iterations of 150 will reduce the computation time, it is not advised to do so as this could prevent the gradient descent from converging to the optimum. An additional benefit of the stochastic gradient descent is that it can potentially avoid getting converging to local minima due to the presence of local anomalies in the fiber bundles (e.g. injured or missing tracts). As the stochastic gradient descent selects a small, different set of random seed points at each iteration. This way, a different random part of the brain is sampled at each iteration, avoiding getting stuck in the local optima that could result from local anomalies, should the complete brain be sampled. Using our current Matlab routines, the column permutation and sign correction took on average 23 seconds, while the full correction of 3D rotations required on average 27 minutes (3.5 GHz quad core desktop processor). Note that the current algorithm was written in Matlab, which is not the ideal platform to perform an iterative algorithm such as fiber tracking, and further speed ups could be obtained by moving to C or C++. Note that our current method uses DTI tractography, which is sub optimal in the presence of crossing fibers (Tuch et al., 2002). While using a tractography method that resolves crossing fibers will clearly result in different metric values, we expect the optima to be found at the same locations. However, DTI tractography can generally be performed much faster than any of the HARDI alternatives and is compatible with a large set of acquisition schemes (simple DTI schemes (Jones et al., 1999), HARDI schemes (Tuch et al., 2002), multi-shell schemes such as those used for diffusion kurtosis imaging (Jensen et al., 2005), and even full blown Cartesian-sampling (Wedeen et al., 2005)). So despite the fact that our approach uses DTI tractography, the recovered gradient orientations can be applied in general DWI post-processing. 22
In this work, we did not directly investigate the One of the benefits of the stochastic gradient descent as an optimizer, is that is robust against such local anomalies. The stochastic gradient descent works by selecting a small, different set of random seed points at each iteration. This way, a different random part of the brain is sampled at each iteration, avoiding getting stuck in the local optima that could result from local anomalies, should the complete brain be sampled. Acknowledgements This work was supported by the Interuniversity Attraction Poles Program (P7/11) References Arsigny, V., Fillard, P., Pennec, X., Ayache, N., 2006. Log-Euclidean metrics for fast and simple calculus on diffusion tensors. Magnetic Resonance in Medicine 56, 411–421. Assaf, Y., Pasternak, O., 2008. Diffusion tensor imaging (DTI)-based white matter mapping in brain research: a review. Journal of Molecular Neuroscience 34, 51–61. Basser, P.J., 1995. Inferring microstructural features and the physiological state of tissues from diffusion-weighted images. NMR in Biomedicine 8, 333–344. ¨ Basser, P.J., Ozarslan, E., 2010. Anisotropic diffusion: from the apparent diffusion coefficient to the apparent diffusion tensor, in: Jones, D.K. (Ed.), Diffusion MRI: theory, methods, and applications. Oxford University Press. Basser, P.J., Pajevic, S., 2000. Statistical artifacts in diffusion tensor MRI (DT-MRI) caused by background noise. Magnetic Resonance in Medicine 44, 41–50. Basser, P.J., Pajevic, S., Pierpaoli, C., Duda, J., Aldroubi, A., 2000. In vivo fiber tractography using DT-MRI data. Magnetic Resonance in Medicine 44, 625–632. 23
Basser, P.J., Pierpaoli, C., 1996. Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI. Journal of Magnetic Resonance Series B 111, 209–219. Conturo, T.E., Lori, N.F., Cull, T.S., Akbudak, E., Snyder, A.Z., Shimony, J.S., McKinstry, R.C., Burton, H., Raichle, M.E., 1999. Tracking neuronal fiber pathways in the living human brain. Proceedings of the National Academy of Sciences of the United States of America 96, 10422–10427. Gudbjartsson, H., Patz, S., 1995. The Rician distribution of noisy MRI data. Magnetic Resonance in Medicine 34, 910–914. Hasan, K.M., Basser, P.J., Parker, D.L., Alexander, A.L., 2001. Analytical computation of the eigenvalues and eigenvectors in DT-MRI. Journal of Magnetic Resonance 152, 41–47. Jenkinson, M., Beckmann, C.F., Behrens, T.E., Woolrich, M.W., Smith, S.M., 2012. Fsl. NeuroImage 62, 782–790. Jensen, J.H., Helpern, J.A., Ramani, A., Lu, H., Kaczynski, K., 2005. Diffusional kurtosis imaging: the quantification of non-gaussian water diffusion by means of magnetic resonance imaging. Magn Reson Med 53, 1432–1440. Jeurissen, B., Leemans, A., Jones, D.K., Tournier, J.D., Sijbers, J., 2013. Investigating the prevalence of complex fiber configurations. Human Brain Mapping 34, 2747–2766. Johansen-Berg, H., Behrens, T.E.J., 2006. Just pretty pictures? What diffusion tractography can add in clinical neuroscience. Current Opinion in Neurology 19, 379–385. Jones, D.K., 2003. Determining and visualizing uncertainty in estimates of fiber orientation from diffusion tensor MRI. Magnetic Resonance in Medicine 49, 7–12. Jones, D.K., Cercignani, M., 2010. Twenty-five pitfalls in the analysis of diffusion MRI data. NMR in Biomedicine 23, 803–820. Jones, D.K., Horsfield, M.A., Simmons, A., 1999. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magnetic Resonance in Medicine 42, 515–525. 24
Klein, S., Staring, M., Pluim, J.P., 2007. Evaluation of optimization methods for nonrigid medical image registration using mutual information and Bsplines. IEEE Trans Image Process 16, 2879–2890. Leemans, A., Jones, D.K., 2009. The B-matrix must be rotated when correcting for subject motion in DTI data. Magnetic Resonance in Medicine 61, 1336–1349. Leemans, A., Sijbers, J., Verhoye, M., Van der Linden, A., Van Dyck, D., 2005. Mathematical framework for simulating diffusion tensor MR neural fiber bundles. Magnetic Resonance in Medicine 53, 944–953. Mori, S., Crain, B.J., Chacko, V.P., Van Zijl, P., 1999. Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Annals of Neurology 45, 265–269. Mori, S., Zhang, J., 2006. Principles of diffusion tensor imaging and its applications to basic neuroscience research. Neuron 51, 527–539. Moseley, M.E., Cohen, Y., Kucharczyk, J., Mintorovitch, J., Asgari, H.S., Wendland, M.F., Tsuruda, J., Norman, D., 1990. Diffusion-weighted MR imaging of anisotropic water diffusion in cat central nervous system. Radiology 176, 439. Pajevic, S., Pierpaoli, C., 1999. Color schemes to represent the orientation of anisotropic tissues from diffusion tensor data: application to white matter fiber tract mapping in the human brain. Magnetic Resonance in Medicine 42, 526–540. Robbins, H., Monro, S., 1951. A stochastic approximation method. The Annals of Mathematical Statistics , 400–407. Stejskal, E.O., Tanner, J.E., 1965. Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. The Journal of Chemical Physics 42, 288. Tax, C., Jeurissen, B., Vos, S.B., Viergever, M.A., Leemans, A., in press. Recursive calibration of the fiber response function for spherical deconvolution of diffusion mri data. NeuroImage .
25
Tournier, J.D., Calamante, F., Connelly, A., 2012. Mrtrix: Diffusion tractography in crossing fiber regions. International Journal of Imaging Systems and Technology 22, 53–66. Tournier, J.D., Mori, S., Leemans, A., 2011. Diffusion tensor imaging and beyond. Magnetic Resonance in Medicine 65, 1532–1556. Tuch, D.S., Reese, T.G., Wiegell, M.R., Makris, N., Belliveau, J.W., Wedeen, V.J., 2002. High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity. Magnetic Resonance in Medicine 48, 577–582. Wedeen, V.J., Hagmann, P., Tseng, W., Reese, T.G., Weisskoff, R.M., 2005. Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging. Magnetic Resonance in Medicine 54, 1377–1386.
26
Appendix
(a) data set A
(b) data set B
(c) data set C
(d) data set D original before correction column permutation and sign correction full correction
Figure 1: Gradient correction on in-vivo data sets with a gradient table subject to: (a) a column permutation (data set A); (b) a column sign change (data set B); (c) a large, random 3D rotation (data set C); and (d) a subtle 3D rotation (data set D). True gradient orientations are displayed as black circles; corrupted gradients as red squares; gradient orientations corrected by column permutation and sign correction as blue plus signs; and gradient orientations corrected by the full rigid correction as green crosses.
27
(IGHLIGHTS
· · · · ·
Our method aligns the diffusion gradient orientations to the corresponding data. We optimize a metric based on whole brain fiber tractography length. Stochastic optimization is used to ensure fast and robust realignment. The method is fully automated, data driven and independent of scanner meta-data. It can be applied to data sets with different acquisition settings and properties.
'RAPHICAL !BSTRACT