Subject–Motion Correction in HARDI Acquisitions: Choices and Consequences

Diffusion-weighted imaging (DWI) is known to be prone to artifacts related to motion originating from subject movement, cardiac pulsation, and breathing, but also to mechanical issues such as table vibrations. Given the necessity for rigorous quality control and motion correction, users are often left to use simple heuristics to select correction schemes, which involves simple qualitative viewing of the set of DWI data, or the selection of transformation parameter thresholds for detection of motion outliers. The scientific community offers strong theoretical and experimental work on noise reduction and orientation distribution function (ODF) reconstruction techniques for HARDI data, where post-acquisition motion correction is widely performed, e.g., using the open-source DTIprep software (1), FSL (the FMRIB Software Library) (2), or TORTOISE (3). Nonetheless, effects and consequences of the selection of motion correction schemes on the final analysis, and the eventual risk of introducing confounding factors when comparing populations, are much less known and far beyond simple intuitive guessing. Hence, standard users lack clear guidelines and recommendations in practical settings. This paper reports a comprehensive evaluation framework to systematically assess the outcome of different motion correction choices commonly used by the scientific community on different DWI-derived measures. We make use of human brain HARDI data from a well-controlled motion experiment to simulate various degrees of motion corruption and noise contamination. Choices for correction include exclusion/scrubbing or registration of motion corrupted directions with different choices of interpolation, as well as the option of interpolation of all directions. The comparative evaluation is based on a study of the impact of motion correction using four metrics that quantify (1) similarity of fiber orientation distribution functions (fODFs), (2) deviation of local fiber orientations, (3) global brain connectivity via graph diffusion distance (GDD), and (4) the reproducibility of prominent and anatomically defined fiber tracts. Effects of various motion correction choices are systematically explored and illustrated, leading to a general conclusion of discouraging users from setting ad hoc thresholds on the estimated motion parameters beyond which volumes are claimed to be corrupted.

1. Introduction

Diffusion-weighted (DW)-MRI enables probing the fiber architecture of biological tissues – in vivo – by encoding the microscopic direction and speed of the diffusion of water molecules (4), while reflecting the amount of hindrance experienced by such molecules along the axis of the applied diffusion gradient due to barriers and obstacles imposed by micro-structures (5). Today, diffusion tensor imaging (DTI) is the method of choice for most neuroimaging studies, e.g., autism (6), schizophrenia (7), and Huntington’s disease (8). Nonetheless, DTI assumes a homogeneous axon population inside a single voxel (9) and fails at modeling more realistic heterogeneous populations. High angular resolution diffusion imaging (HARDI) (10), on the other hand, allows the diffusion acquisition to focus on the angular component of the DW signal using strong gradients and long diffusion times (5), while revealing the intra-voxel orientational heterogeneity, such as crossing and merging fiber bundles. The promising potential of HARDI-based DW-MRI in describing fiber tracts within the human brain comes with a price tag of a wide variety of artifacts related to the gradient system hardware, pulse sequence, acquisition strategy, and subject motion (11). Such artifacts render the quality of diffusion imaging questionable and reduce the accuracy of findings when left uncorrected (1).

1.1. Motion Artifacts

In today’s clinical DW-MRI acquisitions, the presence of the long and strong gradient pulses have made diffusion MRI more sensitive to the detrimental effects of subject motion than other MRI techniques (9, 12, 13). During a scanning session, the degree of a patient’s cooperation may vary: elderly people who may become uncomfortable during large scanning sessions, patients in pain who become restless and agitated during a scan, and unsedated pediatric subjects who will not cooperate long enough to be imaged without motion artifacts. Hence, it is safe to assume that there are always motion artifacts in any given DW-MRI acquisition due to the increased likelihood of involuntary subject motion; especially with HARDI acquisitions, which use a large number of gradient directions resulting in longer scan times. A proof-of-concept of this hypothesis is presented in section 1.

Motion artifacts range from physiological motion (e.g., cardiac pulsation and respiration) to physical (voluntary or involuntary) bulk movement by the patient (14). Physiological motion can be controlled by gating or in the sequence design (15), but the patient bulk movement during the diffusion-encoding gradient pulses leads to severe signal perturbation (16–18), which results in a significant signal phase shift or signal loss (19). The effects of bulk motion are twofold: slow bulk motion can cause misalignment of diffusion data between subsequent gradient applications (i.e., DWI-volumes), resulting in an underestimation of diffusion anisotropy (4), whereas fast bulk motion during the application of a single diffusion gradient causes inhomogeneous signal dropout/attenuation artifacts in the diffusion-weighted images. This dropout effect arises due to signal dephasing within the voxels (13, 14), which is the very phenomenon that gives rise to the DW-MRI contrast, leading to an overestimation of diffusion anisotropy (4). Although misalignment can be tackled by registration-based correction methods (20), the signal dropout due to intragradient motion will persist (4), where such images are identified and excluded from further processing and/or scheduled for reacquisition during the same scan (13, 14, 21–23). Left uncorrected, motion-corrupted datasets introduce bias in the subsequent findings due to the induced variability of diffusion MRI measurements, while affecting the statistical properties of diffusion derived measures in heterogeneous brain regions.

1.2. Motion Correction Choices

The identification and elimination of slow bulk motion artifacts in HARDI data, which are characterized by a high b-value and low signal-to-noise (SNR) ratio, still remains a challenge. In order to allow correction approaches to proceed with reasonable accuracy, motion occurring between diffusion gradients can be treated as if it occurred all at once (24).

Motion effects can be reduced by real-time motion control during the acquisition (a.k.a. prospective motion correction) (25–27), where the acquisition and the source of motion are synchronized, so that the data are never corrupted. In addition, the development of accelerated acquisition methods [e.g., Ref. (28)] can reduce the duration of a scan to minimize the susceptibility of subject motion. A comfortable padding can also be used to minimize head motion while urging the participant to remain without movement (11). Nonetheless padding is not always effective in studies involving infants [e.g., autism diagnosis (29)], where remaining still in the scanner may be more challenging. Nevertheless, prospective methods for motion correction might affect the acquisition time due to the reacquisition of motion-corrupted gradients (14). Such methods might also require external optical tracking systems (23), free-induction decay navigators (26), or volumetric navigators (30), which are not always available on current scanners (27), coupled with the need of time-consuming calibration steps prior to their use (14). Furthermore, rapid modification of diffusion gradients may induce eddy current artifacts (13), and there is no guarantee that the head will move back to the original position.

Motion compensation can also be performed as a post-processing step after acquisition, i.e., retrospective, to guarantee voxel-wise correspondence between different DWIs referring to the same anatomical structure. A common practice is to heuristically select transformation parameter thresholds for detection of motion outliers, where registration and interpolation are applied to gradient directions that are claimed to be corrupted. Software packages for image-based registration of DWIs are becoming readily available, e.g., FSL-MCFLIRT (2, 31), the Advanced Normalization Tools (ANTS) (32), TORTOISE (3), and BRAINSFit (33) employed in DTIPrep (1).

A typical retrospective motion correction algorithm involves two stages (20): first, finding the global transformation parameters that would transform all DWIs to the same coordinate frame, and then, applying the estimated transformations to the diffusion data. Solving for the transformation parameters usually involves rigidly registering the DWIs to a reference volume representing the same anatomical structure, but without being contaminated by motion artifacts. Examples of such a reference include a T2-weighted image (16), or a non-diffusion-weighted image (a.k.a baseline with b-value = 0) due to its high SNR and lesser vulnerability to eddy current distortion (34), where the difference in intensity profiles is compensated for using normalized mutual information similarity measure. Another alternative is a model-based reference volume computed for each diffusion-weighted image based on tensor fitting (35, 36). Model-based motion correction implicitly assumes that the original position defined by the baseline volume is the reference position to be aligned to Sakaie and Lowe (20). Recently, it has been shown that model-based motion correction becomes a more powerful choice for correcting higher b-value diffusion imaging, which does not contain enough anatomical features to be registered accurately (36).

Applying the estimated transformation parameters is performed using interpolation, which computes intensities at transformed voxel coordinates as a weighted sum of the scaled intensities at surrounding voxels. The diffusion gradient vectors are also reoriented to incorporate the rotational component of subject motion (37). Interpolation is usually carried out by an exact fit of a continuously defined model to discrete data samples. Nonetheless, this exact fit is less appropriate when data are noise-corrupted, since the model is forced to fit the noise too. Although using regularized interpolation can tackle noisy data, it is only preferable to applying denoising followed by standard interpolation under the assumption that the signal is a stationary Gaussian process (38); a situation that is not applicable for diffusion-weighted images, which are contaminated by Rician noise. Based on the central limit theorem, the (weighted) average of a large set of i.i.d. samples tends to follow a normal distribution. Thus, interpolation between Rician distributed samples might change the distribution toward a Gaussian PDF (39). We can, therefore, argue that the denoising process decreases the effect of standard interpolation on altering the underlying data distribution.

Another retrospective approach is to cast motion correction as an outlier rejection process, ranging from simply excluding one or more gradients bearing strong motion artifacts beyond acceptable levels of motion (11, 14, 40), to statistical methods for detecting and discarding voxel-wise diffusion measurements as outliers (17, 41, 42). Usually discarding entire scans (a.k.a motion scrubbing in functional MRI) either can be performed by visual inspection or based on predefined thresholds on estimated motion parameters (4). Nevertheless, removing gradients limits the ability to reconstruct crossing fibers, especially at small separation angles, due to the decreased number of distinct gradient directions needed for diffusion reconstruction. Moreover, scrubbing would introduce intersubject SNR and bias differences that would in turn affect subsequent statistical analysis (1). On the other hand, local exclusion of corrupted voxels for robust diffusion reconstruction in the presence of outliers is based on the deviation of the observed measurements (usually after motion correction) from the assumed diffusion model. Using these approaches for motion correction itself would mingle the effect of being an outlier to an assumed model with that of being corrupted due to motion. Further, local exclusion would lead to a different number of DWIs locally available for each voxel, complicating subsequent analysis to avoid bias due to different SNR values for different brain regions (1).

A common concern with retrospective methods in clinical studies, whether registration-based and/or outlier-based, is that data with different levels of motion will be subject to different schemes of motion correction. For instance, patients may show more motion than controls, or sedated subjects may be different from non-sedated. Applying different motion correction schemes could introduce a confounding factor for statistical analysis of populations that show different motion patterns. Nonetheless, eyeballing the acquired/preprocessed DWIs prior to proceeding to further analysis is highly recommended.

1.3. Objective and Contributions

The lack of a comprehensive/rigorous quality control (QC) for HARDI datasets can result in considerable error and bias in subsequent analyses, which may affect research studies using these datasets. Most current software packages such as DTIPrep (1), TORTOISE (3), and FSL (2), which offer various tools for processing and analysis of diffusion-weighted images, are mostly limited to DTI datasets, which are characterized by low b-values (i.e., higher SNR) and fewer gradients (i.e., shorter acquisition times). Nonetheless, special care is needed for HARDI datasets due to their low SNR and longer acquisition times, which increase the likelihood of subject motion. As a part of a thorough pipeline for HARDI-QC, this paper addresses the motion correction aspect for slow bulk motion where users often do not fully understand the consequences of different types of correction schemes on the final analysis, and whether those choices may introduce confounding factors when comparing populations. Therefore, the presented work is directed toward clear guidelines and recommendations to the standard users in practical settings.

The optimal preprocessing pipeline for HARDI sequences remains an open question and a challenge for real data. Questions that might arise include: is there a threshold that would identify a motion-corrupted volume? How sensitive are HARDI reconstructions to such a predefined threshold? What is the impact of various motion correction schemes on subsequent HARDI-based reconstructions and tractography? So far, these questions have received, surprisingly, little attention in various DW-MRI studies of clinical populations. This study, then, focuses on the effect of preprocessing schemes, in particular motion correction, commonly deployed as a post-acquisition step, on succeeding steps. We propose a comprehensive experimental framework (see Figure 1) that enables making use of human brain HARDI data from a well-controlled motion experiment to simulate various degrees of motion/noise corruption. The comprehensiveness is related to the systematic evaluation of the outcome of different motion correction choices commonly used by the scientific community on different DWI-derived measures. To our knowledge, this evaluation does not exist in the literature and has not been discussed in detail.

FIGURE 1

Figure 1. A comprehensive experimental framework for subject motion simulation to systematically evaluate the outcome of different motion correction choices commonly used by the scientific community on HARDI-based reconstructions and tractography. (A) A human brain HARDI data were acquired from a well-controlled motion experiment. (a.1) Acquired DWIs were preprocessed to obtain nearly noise-free and motion-free datasets. (a.2) For automated tractography selection and the quantification of whole brain connectivity, a subject-specific unbiased atlas was constructed via DTI-derived data from HARDI sequences resulting in a tensor atlas, where we can define a detailed parcelation of neuroanatomical structures, and map it back to each raw scan. (B) Noticeable motion was then simulated by randomly mixing gradients from the acquired datasets. (C) Motion correction involves four main decision variables where each distinct combination of choices defines a correction scheme. (D) Reconstruction of a corrected or motion-free dataset entails reconstructing the voxel-wise fiber orientation distribution functions, detecting local (voxel-wise) fiber orientation, preforming whole brain tractography, and automatically selecting anatomical pathways. (E) The evaluation of the effect of a motion correct scheme has been investigated based on voxel-wise metrics, global brain connectivity metric, and tract-based metric.

Choices for correction include exclusion or registration of motion corrupted directions, with different choices of interpolation, as well as the option of registration/interpolation of all directions versus corrupted directions only. The effect of denoising as a preprocessing step applied prior to motion correction is also investigated. Further, the choice of the reference volume used in the registration framework is also discussed. The comparative evaluation covers four metrics: (1) the similarity of fiber orientation distribution functions (fODFs) via Jensen–Shannon divergence (JSD), (2) the deviation of multiple fiber orientations at each voxel, (3) the global brain connectivity via graph diffusion distance (GDD), and (4) the reproducibility of seven anatomically defined fiber pathways via Cohen’s Kappa statistics. On the basis of our findings, we recommend assuming that motion is inevitable, even subtle, in the acquired scans. Motion correction, therefore, needs to be applied to all gradient directions without relying heuristically on a threshold that determines a gradient direction to be claimed as motion corrupted.

2. Materials and Methods

2.1. Motion is Inevitable: Proof-of-Concept

To back up our assumption that motion is omnipresent, we analyzed data from three healthy human phantoms (males 30–40 years old) visiting each of the four clinical sites (Chapel Hill, Philadelphia, St. Louis, and Seattle) as a part of the ACE-IBIS study [Autism Centers for Excellence, Infant Brain Imaging study (6)], using a total of six MRI systems (two sites using both research and hospital scanners). All study procedures were approved by the institutional review board at each clinical site, and informed, written consent was obtained for all participants. In addition, the traveling phantoms sign consent forms at each of the sites, as per their own institutional IRBs. The sites include the University of Washington, Seattle, the Washington University in St. Louis, the Childrens hospital of Philadelphia, and the University of North Carolina at Chapel Hill. Each subject was scanned twice on a 3 T Siemens Tim Trio scanner¹ with a strict calibration of image acquisition parameters. Test–retest reliability at each site was established with two scans within 24 h. The scans were acquired within 1 week to guarantee that there were no major brain changes over time. The scanning environment was well controlled. Comfortable padding was used to minimize head motion, and patients were urged to remain without movement. Eddy current was compensated for using a Twice Refocused Spin Echo (TRSE) protocol², with FoV = 209 mm, 76 transversal slices, thickness = 2 mm (2 mm)³ voxel resolution, matrix size = 106 × 106, TR = 11100 ms, TE = 103 ms, one baseline image with zero b-value and 64 DWI with b-value at 2000 s/mm², with a total scan time of 12.5 min.

Initially, we ran automated Quality Control on the DWIs via DTIPrep (1), which includes among other steps interlaced correlation analysis for detection and removal of fast bulk motion within a single DWI volume, where no quantitative within-gradient motion was detected. Inspired by Sakaie and Lowe (20), FSL-MCFLIRT (31) was then used to provide the rigid transformation matrix (i.e., six degrees of freedom) for each volume having the baseline image as the reference for motion correction and normalized mutual information as the cost function. It is worth noting that MCFLIRT employs a global-local hybrid optimization method for robust affine registration that is specifically tailored to brain images. Within a multiresolution framework, four scales were used (8, 4, 2, and 1 mm, i.e., supervoxel vs. subvoxel). At each scale, volumes were resampled after initial filtering to reduce the effect of noise. Further, we tested motion correction based on denoised HARDI sequences using the Joint Rician LMMSE filter (43) implemented as part of 3D Slicer (www.slicer.org), and found that the quantified motion with and without noise reduction was very similar.

To quantify motion, we used the magnitude of the translation vector (in millimeters) as well as the axis–angle rotation representation (in degrees) (4). The boxplots in Figure 2 show the rotational and translational components of the motion being detected from a total of 24 DWI datasets, showing an average of 0.39° rotation and 0.61 mm translation. The graphs in Figure 2 illustrate the arbitrariness of a common calculation of percentage of motion correction to determine the number of affected scans, here shown as a function of thresholding on the estimated motion parameters. While this experiment attributes the estimated rotation and translation parameters to actual subject motion, a part of the experimentally obtained parameters may be due to some imaging/image-processing uncertainty and also to image differences due to anatomical properties of the object (e.g., tissue orientation) that make the images “look” different even if they were perfectly aligned. To backup our analysis, we conducted another experiment where we contaminated a single DWI dataset with two independent realizations of Rician noise such that the two generated DWI images were perfectly aligned because they were the exact same image. Then, we ran motion correction where all DWI images were aligned to the same baseline, we obtain similar motion parameters although we are registering two independent acquisitions of the same subject. We therefore conclude that the transformation parameter estimates from FSL-MCFLIRT (31) are resilient to noise and may primarily caused by subject motion during a DWI scan, or eventually also by relative motion between subject and scans if considering artifacts due to pulse sequence and scanner technology.

FIGURE 2

Figure 2. Average and standard deviation of the percentage of motion-corrupted gradient directions as a function of thresholding on the estimated rotation angle in degrees (left) and the estimated translation magnitude in millimeter (right) for three human phantoms scanned twice at four clinical sites. The boxplots show the overall statistics of estimated motion parameters.

2.2. Living Phantom: Acquisition and Gold Standard Generation

Unlike conventional MRI, where realistic phantoms exist (44), there is no widely acceptable realistic DWI phantom for the assessment of different processing tasks (43, 45). Existing phantoms simulate crossing sections in two and three dimensions, but they are not representative of white matter complex architecture with multiple fiber crossing, bending, and branching. The lack of realistic phantoms motivates us to base our analysis on living (human) phantoms being scanned under well-controlled environments and propose a HARDI-based QC to yield motion- and noise-free datasets. Acquired DWIs were preprocessed (refer to Figure 1A) to obtain nearly noise-free and motion-free datasets according to the following pipeline, and therefore, to be used as a gold standard for reconstruction and tractography.

2.2.1. HARDI-based quality control

The QC process starts with identifying individual volumes having fast bulk (intra/within-gradient) motion using the signal dropout score proposed in Ref. (14). The score was computed for each slice in each volume, where slices with a score >1 were considered to have suspect signal dropout. Based on a zero-tolerance strategy, any volume having at least one slice with signal dropout was excluded from further analysis. It is worth noting that no within-gradient motion was detected in our phantom acquisitions. Each gradient was then independently denoised to reduce noise using the Rician LMMSE estimator with an 11 × 11 neighborhood (46) implemented in 3D Slicer (www.slicer.org) where the noise parameter is automatically estimated. Using DTIPrep (1), interslice brightness artifacts were detected via normalized correlation analysis between successive slices within a single DWI volume, where corrupted gradients were excluded before being streamed into the next steps. Further, interlaced correlation analysis (1) was used for detection and removal of Venetian blind artifacts (seen when motion occurs between the interleaved parts of an individual gradient volume) and fast bulk motion within a single DWI volume, where no quantitative within-gradient motion was detected.

For each DW-MRI scan, iterative FSL-MCFLIRT (20) was used to correct for intergradient subtle motion (<1° rotation and < 0.8 mm translation), with the baseline volume as the reference for rigid alignment (i.e., six degrees of freedom with normalized mutual information as the cost function). The corresponding diffusion-weighting gradient vectors were reoriented accordingly (37). To palliate the effect of spatial intensity inhomogeneities, N4 correction (47) was performed where the bias field was computed from the baseline volume and subsequently applied to all diffusion-weighted images. For further noise reduction, the Joint LMMSE (43) (www.slicer.org) was used to exploit the joint information from neighboring gradients from motion-corrected sequences. To avoid over-blurring, we used a 2 × 2 × 2 neighborhood with six neighboring gradients.

2.2.2. Atlas-guided parcelation

For automated tractography selection and the quantification of whole brain connectivity, we defined a subject-specific unbiased atlas via DTI-derived data from HARDI sequences belonging to the same subject/phantom. This results in a tensor atlas, where we can define a detailed parcelation of neuroanatomical structures, and map it back to each raw scan. This reduces registration variability between each phantom data when defining the parcelation in subject spaces. The full process entails atlas creation and parcelation definition, as detailed in the following.

2.2.2.1. Co-registration and atlas building. To define a common reference space, our framework is centered around the creation of a DTI atlas, generated as an unbiased average atlas from the study dataset via a deformable atlas building strategy. Unbiased atlas building is used to provide one-to-one mapping between the image data and the template atlas, wherein the atlas is built from the population of data as the centered image with the smallest deformation distances. The overall registration framework, similar to what has been presented in Ref. (48), proceeds in four steps: (1) image preprocessing via skull-stripping and tensor estimation, (2) affine alignment, (3) unbiased diffeomorphic atlas computation via GreedyAtlas module in AtlasWerks³ software (49), and (4) a refinement step via symmetric diffeomorphic registration using the advanced normalization tools – ANTS (32).

Image preprocessing. A brain masking is first performed on the baseline images using FSL-BET2 (brain extraction tool) (50) to remove all non-brain parts of the image. BET2 uses a surface model approach to robustly and accurately carry out the segmentation. We then model tensors using the brain masks from the initial DWI datasets by using weighted least squares estimation, and then extract related scalar maps such as fractional anisotropy (FA) images.

Affine alignment. The second step applies affine registration of baseline images to a previously defined baseline template. A multithreaded, coarse-to-fine registration scheme using mattes mutual information metric is employed in that regard (33). The transformations are applied to curvature FA maps. The use of curvature FA as feature to derive registration has initially been presented by Goodlett et al. (51). It is defined as the maximum eigenvalue of the Hessian of the FA image, therefore measuring image intensity curvature (second derivative) in the direction of largest curvature, which acts like a 3D ridge detector. It is computed by convolution of the FA image with a set of Gaussian second derivatives with a fixed aperture, proportional to the size of the white matter structures. The curvature feature image proved to be an efficient detector of the 3D manifold skeleton of major fiber bundles, which occur as tubular or sheet-like thin structures (similarly to the TBSS software), with the strongest response at their center. It is thus commonly used by our group when building population atlases to optimize correspondence of fiber tract geometries, and integrated into our freely distributed software package (48). The curvature FA maps are thus mapped to this template space, and then the intensity is rescaled via histogram matching.

Atlas building. We then use an unbiased deformable atlas-building procedure (52) that applies large deformation diffeomorphic metric mapping transformations to these intensity rescaled mapped curvature FA images. The procedure relates individual datasets to the subject-specific atlas template space by means of non-linear, invertible transformation. Tensor maps are transformed into the atlas space with tensor reorientation by the finite strain approach (53), taking into account both affine transformation and non-linear deformation. The transformed tensor images are finally averaged using the Riemannian framework proposed in Fletcher and Joshi (54), resulting in a three dimensional average tensor atlas.

Atlas refinement. An additional step is performed by direct symmetric diffeomorphic registration of initial FA images to the previously created DTI-FA atlas via the Advanced Normalization Tools – ANTS (32). In our experience, this dual stage procedure has been shown to produce a sharper atlas with improved registration accuracy, most likely attributable to the use of local normalized cross-correlation as the image similarity metric. Final affine transformation and deformation fields are then available from subject space to atlas space.

2.2.2.2. White matter parcelation. We used the publicly available JHU-DTI-SS (a.k.a. “Eve”) atlas described in Oishi et al. (55). Defined as a single subject template, it includes both structural (T1w, T2w) and DTI images with white matter map parcelations, defining 176 hand-segmented core and peripheral regions of interest (ROIs). A multithreaded, coarse-to-fine diffeomorphic registration scheme using the cross-correlation metric via ANTS is employed on FA images between the Eve atlas and the subject-specific atlas. The computed deformation field is then applied to the Eve white matter label map. We can then map the parcelation, now defined in our subject atlas space, back to raw data in the initial image space, via the use of previously computed displacement fields. On a specific note, we concatenated the transformations from Eve atlas space to our initial images in order to directly map the parcelation and avoid the use of multiple interpolations. The white matter parcelation map is then defined both in the subject-specific atlas space and in each individual subject space.

2.3. Subject Motion: Between Simulation and Correction

2.3.1. Human motion simulation

As a pilot study, one human phantom was asked to be rescanned with his head tilted to simulate noticeable motion. The two datasets, after being QCed (see 1), were then used to construct motion-corrupted sequences (see Figure 1B). Based on the alignment of the baseline images of the two scans (original and tilted) using FSL-MCFLIRT, about 12° of rotation and 7 mm of translation were detected, whereas <1° of rotation and 0.8 mm of translation were found when aligning individual DWIs to their corresponding baseline image. It is worth noting that the quantified motion between the acquired datasets (i.e., untilted versus tilted brains) can be classified as severe subject motion (36). We then arbitrarily considered the first out of the two scans as the “motion-free” sequence and used it as a reference for performance evaluation of different motion correction schemes. A random percentage of DW images (10, 30, 50, 70, and 90%, each with five distinct random sets of gradient directions) drawn from the second scan (tilted brain) were mixed with the first scan to construct 25 motion-corrupted datasets. Noisy sequences were generated by simulating Rician noise based on seven levels of SNRs from 4 to 20 (56), yielding 175 (5 experiments × 5 corruption percentages × 7 SNR levels) sequences.

2.3.2. Motion correction schemes

Correction for subject motion involves four main decision variables (see Figure 1C), where each distinct combination of choices defines a motion correction scheme. The first variable is which reference volume is to be used in the alignment process. Two options are available (20): baseline-based [e.g., Ref. (16)] and model-based [e.g., Ref. (35, 36)]. In this context, we use the FMAM (Fit Model to All Measurements) method (35) where target images for registration were generated by first fitting the diffusion tensor to the DWIs, followed by diffusion simulation to provide target images of similar contrast to the DWIs. Notice that with >50% motion corrupted, model-based reconstruction infers the spatial position/orientation from the gradients corresponding to the tilted brain due to its majority (i.e., gradients of the untilted brain are considered the motion-corrupted directions). Therefore, with model-based correction for sequences having more than 50% corrupted directions, the tilted brain was used as a reference for performance evaluation.

The second variable denotes whether the correction is performed based on raw or denoised DWIs, where the denoising process should not take into account joint information between diffusion gradients due to motion corruption. In our experiments, we denoised motion-corrupted sequences using the Rician LMMSE estimator (46), where each gradient was independently denoised.

The third variable entails the mode of correction, i.e., registration-based versus outlier-based. The first choice explores two options: (1) only aligning and interpolating the corrupted gradient directions to mimic the situation where a predefined motion parameter threshold is used to claim whether a DWI volume is motion-corrupted, (2) assuming there is always motion, forcing the alignment and interpolation of all DWI volumes. Note that both options involve the reorientation of the diffusion gradient vectors corresponding to the corrupted volumes (37) to incorporate the rotational component of subject motion. In the second choice, i.e., outlier-based, we mimic the motion scrubbing approach, where we exclude the affected gradient directions from subsequent computations (i.e., diffusion profile reconstruction and tractography). Eventually, the interpolation step in the registration-based choices introduces the fourth variable where we study the impact of using trilinear and sinc interpolants.

It is important to stress that, in our motion simulation paradigm (i.e., randomly mixing DW volumes from a tilted-brain dataset), the identity of the motion-corrupted directions is known a priori without any use of parameters. This prior information is used via the outlier-based correction, as well as the interpolate corrupted directions choices. Nonetheless, in practice, this a priori information corresponds to heuristically set thresholds on the estimated motion parameters beyond which volumes are claimed to be corrupted/outliers. For example, a rotation threshold of 0.5° and a translation threshold of about one voxel spacing are set by default in DTIprep (1).

2.4. Reconstruction and Tractography

The reconstruction and whole brain tractography were computed for the motion corrected sequences as well as the motion-free sequences [gold standard generated in Section 2, followed by automatic tractography selection for seven major fiber bundles (see Figure 1D)].

We employed the constrained spherical deconvolution (CSD) technique (57) to reconstruct fiber orientation distributions functions (fODFs) from the DWI data using the diffusion imaging Python (DiPy) library (58). The fiber response function was estimated from the corpus callosum region, defined by the white matter parcelation (see 2), where it is known to have single fibers. In particular, we used an ROI at the center of the corpus callosum and of a radius that would include all its voxels. The response function was estimated in that region from the voxels with FA higher than 0.7.

Part of our analysis is based on comparing brain connectivity graphs, which are represented as weighted graphs and computed from fiber tractography results. Whole brain tractography was performed using the EuDX deterministic tracking technique (59), which is implemented in the DiPy library (58), using random seeding inside the brain region and a turning-angle threshold of 30° between two connected voxels [as suggested by Parizel et al. (60) to provide sufficient fiber density while minimizing the number of spurious tracts].

To extract brain connectivity graphs from the fiber tractography results, we used the 176 core and peripheral ROIs defined in the white matter parcelation (see 2). Let N_ij denote the total number of streamlines connecting the i-th and j-th ROIs, each with length $l_k^{\mathit{ij}}\forall k\in \left[1,N_{\mathit{ij}}\right]$, and the edge weights w_ij computed as follows (61): $w_{\mathit{ij}}=\frac{1}{N_{\mathit{ij}}}{\sum }_{k=1}^{N_{\mathit{ij}}}\frac{1}{l_k^{\mathit{ij}}}$. The normalization by the tracts length gives a higher connection strength to short tracts to compensate for the signal attenuation as a function of tract length. It is worth noting that the concept of using the connection strength or other measures to weight the graph edges was previously discussed in several papers [e.g., Ref. (62, 63)].

For tract-based analysis, an automatic tractography selection method was performed to select a subset of detected tracts from the whole brain tractography result corresponding to a specific white matter structure. Starting from the Eve-atlas-based white matter parcelation map defined in the subject space (see 2), the pass-through and not-pass-through volumes of seven fundamental fiber bundles (left and right hemispheres) were defined. To remove fibers that do not belong to the pathway of interest, we used the geometrical constraints specific for different fiber bundles as defined in Ref.(64), where the anatomical characteristics of these fiber bundles are defined in Ref. (65). We report the matching results from seven major fiber bundles: corpus callosum (CC), cingulum of the cingulate gyrus (CG), corticospinal tract (CST), fornix (FX), inferior fronto-occipital tract (IFO), inferior longitudinal fasciculus (ILF), and uncinate fasciculus (UNC).

2.5. Motion Correction Consequences: Evaluation Metrics

The influence of various motion correction choices on subsequent reconstruction and tractography is evaluated according to voxel-based, global connectivity-based as well as tract-based metrics (see Figure 1E), detailed as follows.

2.5.1. Voxel-based metrics

In order to measure similarities between the original motion-free fODFs and the fODFs corresponding to the motion corrected images, we use the Jensen–Shannon divergence (JSD), which has been used to quantify differences between ODFs in various studies, e.g., Ref. (66, 67). Given two probability distributions P and Q, the JSD metric is defined as follows:

$\mathrm{JSD}\left(\mathrm{P}\parallel \mathrm{Q}\right)=\frac{1}{2}\left[D_{\mathrm{KL}}\left(\mathit{P}\parallel M\right)+D_{\mathrm{KL}}\left(\mathit{Q}\parallel M\right)\right],$ (1)

where M = (P + Q)/2 and D_KL is the Kullback–Leibler divergence. In our case, P and Q are represented as discrete distributions; therefore, the KL divergence takes the following form: $D_{\mathrm{KL}}\left(\mathit{P}\parallel Q\right)={\sum }_i{P}_i\log \frac{P_i}{Q_i}$, where i is the discrete sample index. The JSD is for PDFs, but we compute it for normalized fODFs. We believe that it is a good measure since it reveals subtle changes in PDFs so we can also keep track of changes in fiber volumes as well as orientations.

In addition to comparing fODFs, we are interested in quantifying local deviations in fiber orientations due to motion correction. Since brain connectivity maps are inferred by tracking local fiber orientations extracted from fODFs, distortions in those directions may lead to unreliable brain connectivity maps. Therefore, it is important to study the impact of motion correction on fiber orientations by directly comparing the local fiber orientations before and after correction. To that end, we use the mean angular deviation measure θ defined as follows:

${\mathrm{\theta }}_{i,j}^k=\frac{180}{\mathrm{\pi }}\left|{\cos }^{-1}\left(v_i^k.v_j^k\right)\right|,\mathrm{\theta }=\frac{1}{N}{\displaystyle \sum _{k=1}^N}{\mathrm{\theta }}_{i,j}^k,$ (2)

where N is the number of fibers compared, and $v_i^k$ and $v_j^k$ correspond to the orientations of fiber k, with and without motion correction. Before averaging the deviations, we match the fibers, such that fiber j has the closest direction to fiber i. If the number of fibers is different, we compare the fibers that are present in both voxels. For example, if we have three fibers after motion correction, whereas before correction there were only two, we compare the two closest fiber directions. The fiber orientations were computed using the DiPy peak extraction tool (with 0.4 relative peak threshold and 20° minimum separation angle). We allowed up to five orientations in each voxel (N = 5). Since general image transformation does not necessarily preserve the original ordering of the fiber orientations, we first match the fibers based on the angular distance between each pair before computing the mean deviation.

2.5.2. Global connectivity-based metric

Once the brain connectivity graphs were generated for the different sequences, we compared them by means of the graph diffusion distance (GDD) metric, which has been proposed in Ref. (68). The GDD is a novel distance measure for comparing weighted graphs, which takes into account the graph structure in addition to the edge weights, compared to the Frobenius norm, which is sensitive only to the edge weights. For an explanation of the differences between the GDD and the Frobenius norm, see the Barbell graph example in Ref. (68).

The GDD is based on the diffusion maps framework (69). Let A₁ and A₂ be weighted adjacency matrices for N vertices, that is, A₁ and A₂ are symmetric, non-negative, N × N real matrices with zeros along the principle diagonal. The (unnormalized) graph Laplacian operator is defined by L_n = D_n − A_n (for n = 1, 2), where D_n is a diagonal degree matrix for the adjacency A_n, i.e., ${\left(D_n\right)}_{i,i}={\sum }_{j=1}^N{\left(A_n\right)}_{i,j}$.

Given adjacency matrices A₁ and A₂, the columns of the Laplacian exponential kernels, exp(-tL₁) and exp(-tL₂), describe the different diffusion patterns centered at each vertex generated by diffusion up to time t under the two sets of weighted edges. Measuring the sum of squared differences between these patterns, summed over all the vertices, yields

$\begin{array}{llll}\hfill {\mathrm{\xi }}_{\mathit{gdd}}^2\left(A_1,A_2;t\right)& ={\displaystyle \sum _{i,j}}{\left({\left(\exp \left(-t{L}_1\right)\right)}_{i,j}-{\left(\exp \left(-t{L}_2\right)\right)}_{i,j}\right)}^2& \hfill & \\ \hfill & =\left|\right|\exp \left(-t{L}_1\right)-\exp \left(-t{L}_2\right)|{|}_F^2\end{array}$ (3)

where ||⋅||_F is the matrix Frobenius norm. This defines a family of distance measures ξ, depending on the information propagation time t. The graph diffusion distance is given by ξ at the time of maximal difference, i.e., d_gdd(A₁, A₂) = max_t ξ_gdd(A₁, A₂; t). Here, we compute d_gdd(A₁, A₂) by first diagonalizing L₁ and L₂ and using the exponential mapping. Then, Equation (3) allows the computation of ξ(A₁, A₂; t) for any fixed t. Finally, we optimize over t by a line search to give d_gdd(A₁,A₂).

2.5.3. Tract-based metric

The spatial matching between motion-free and motion-corrected tracts was examined using Cohen’s Kappa statistic (70). The streamlines for a specific fiber tract (e.g., CST, IFO.) are first converted to a binary volume with the same dimension and spacing of the raw DWI, where voxels that were occupied by at least one streamline were assigned a value 1. The two tracking results to be matched were then superimposed to identify: (1) voxels that did not contain streamlines in either result (NN), (2) voxels that contain streamlines in both results (PP) and (3) voxel that contain streamlines in one of the results (PN or NP)⁴. The Kappa statistic measures the level of agreement of the tracking results and corrects for agreement expected by chance. Hence, Kappa is computed based on the probability of agreement P(a) and the probability of expected agreement due to chance P(e) as (71),

$\mathrm{\kappa }=\frac{P\left(a\right)-P\left(e\right)}{1-P\left(e\right)},$ (4)

where,

$\begin{array}{llll}\hfill P\left(a\right)& =\frac{\mathit{NN}+\mathit{PP}}{\mathit{PP}+\mathit{PN}+\mathit{NP}+\mathit{NN}},& \hfill & \\ \hfill P\left(e\right)& =\frac{\left(\mathit{NP}+\mathit{PP}\right)\left(\mathit{PN}+\mathit{PP}\right)+\left(\mathit{NP}+\mathit{NN}\right)\left(\mathit{PN}+\mathit{NN}\right)}{{\left(\mathit{PP}+\mathit{PN}+\mathit{NP}+\mathit{NN}\right)}^2}.& \hfill & \end{array}$

3. Results

The fODFs and the whole brain tractography were computed for the 3,150 motion corrected sequences (175 datasets × 18 correction schemes), as well as the motion-free sequences, followed by automatic tractography selection for seven major fiberbundles.

3.1. Voxel-Based Metrics

The average JSD metric was computed using the fODF reconstruction from the “motion-free” dataset, not corrupted by mixing DWI directions from the tilted-brain scan, as a reference (i.e., presenting only subtle motion inherent to a scan). We differentiated between regions where multiple fibers were detected versus single fiber regions. Figure 3 shows the average JSD values for single and multiple fiber regions as a function of motion corrupted percentage for different SNR levels and as a function of SNR levels for different motion corrupted percentages. Figure 4 illustrates sample reconstructions from motion-free versus motion-corrected datasets for different corrupted percentages and different motion correction choices. Table 1 shows the effect of the denoising process prior to applying motion correction on the average JSD values for single and multiple fiber regions as a function of SNR levels for different motion corrupted percentages. Figure 5 shows the average deviation of local fiber orientations (for the first two dominant detected fibers per voxel) as a function of motion corrupted percentage, as well as SNR levels.

FIGURE 3

Figure 3. The average Jensen–Shannon divergence (JSD) values (lower is better) for reconstructions based on raw datasets (denoised ones share similar performance) as (A) a function of motion corrupted percentage for different SNR levels and (B) a function of SNR levels for different motion corrupted percentage. The first and third columns show JSDs single fiber regions while the second and fourth columns show such values for reconstructions based on multiple fiber regions. Notice the impact of motion scrubbing (removing gradient directions), which becomes more significant with more motion-corrupted directions when compared to registration-based correction. Further the impact of motion scrubbing is rendered evident for 10% corrupted gradients.

FIGURE 4

Figure 4. Sample fODFs reconstruction from untilted and tilted motion-free datasets as well as reconstruction from motion-corrected datasets with 10, 30, and 70% corrupted gradient directions. Correction choices shown include outlier-based (i.e., motion scrubbing) and registration-based (using baseline and model-based reference volumes).

TABLE 1

Table 1. The effect of denoising on the average ± standard deviation of Jensen–Shannon divergence (JSD) values for single fiber regions and multiple fiber regions as a function of SNR levels for different motion corrupted percentages.

FIGURE 5

Figure 5. The average fiber orientation deviation (lower is better) for reconstructions based on raw datasets (denoised ones share similar performance) as (A) a function of motion corrupted percentage for different SNR levels and (B) a function of SNR levels for different motion corrupted percentage. The first and third columns show orientation deviation for the first detected fiber having the largest volume fraction while the second and fourth columns show such values for the second detected fiber having the second largest volume fraction. Notice that local fiber orientations are more affected by motion scrubbing as SNR decreases and/or corrupted directions increase.

3.2. Global Connectivity Metric

Figure 6 shows the average graph diffusion distance (GDD) metric as function of both the corrupted directions percentage and the SNR levels. The metric compares the weighted connectivity graphs from the whole brain tractography of the “motion-free” dataset to that of the motion-corrected datasets. It is worth noting that the tractography of the tilted brain dataset is used as a reference for model-based corrections when the corrupted percentage exceeds 50%. Figure 7 visualizes the brain connectivity being represented circularly using the Circos software (72) where the parcelated structures (refer to Table S1 in Supplementary Material for their full names) are displayed on a connectogram representing left and right hemispheres symmetrically positioned along the vertical axis. The weighted connectivity matrix computed as described in Section 4 was normalized to attain a unit maximum. Each entry in the normalized connectivity matrix corresponds to an interregion link with thickness proportional to the entry weight. To avoid dense visualization, all entries with weight <0.15 were discarded.

FIGURE 6

Figure 6. The average graph diffusion distance (GDD) (lower is better) for the whole brain tractography derived from the raw datasets (denoised ones share similar performance) as (A) a function of the corrupted directions percentage for different SNR levels and (B) a function of SNR levels for different motion corrupted percentages. Notice the different behavior displayed by motion scrubbing for ≥50% corrupted directions, which due to having more short tracts connecting nearby region of interests while being assigned to larger weights in the graph construction step.

FIGURE 7

Figure 7. Sample reconstructed connectomic profile (i.e., connectogram) from untilted and tilted motion-free datasets as well as connectograms from motion-corrected datasets with 10, 30, and 70% corrupted gradient directions. Correction choices shown include outlier-based (i.e., motion scrubbing) and registration-based (using baseline and model-based reference volumes). Notice the tendency of motion scrubbing to add more links between nearby ROIs at corruption percentages, implying the detection of more short tracts.

3.3. Tract-Based Metric

Table 2 shows the average Cohen’s Kappa statistic computed for corpus callosum (CC), corticospinal tract (CST), and inferior fronto-occipital tract (IFO) (where other pathways showed similar trend) based on automatic tractography selection using whole brain tractography of raw datasets (denoised datasets showed similar trends due to the robust fODF estimation, yet their graphs were omitted due to space limitation). Figures 8–12 show sample tractography selections for the aforementioned tracts from the untilted motion-free dataset as well as selections from motion-corrected datasets with different corrupted gradient directions. Correction choices shown include outlier-based (i.e., motion scrubbing) and registration-based (using baseline and model-based reference volumes). Also pass-through (in green) and not-pass-through (in red) volumes (i.e., isosurfaces) are shown. Their definitions along with the geometric constraints employed to remove fibers, which do not belong to the pathway of interest, can be found in Ref. (64).

TABLE 2

Table 2. The average Cohen’s Kappa statistic (higher is better) of different anatomically defined fiber pathways (other pathways show similar trend) based on automatic tractography selection based on whole brain tractography of raw datasets (denoised ones share similar performance) for different corrupted directions percentages.

FIGURE 8

Figure 8. Sample tractography selection for the corpus callosum (CC) from the untilted motion-free dataset as well as selections from motion-corrected datasets with 10, 30, and 70% corrupted gradient directions. Correction choices shown include outlier-based (i.e., motion scrubbing) and registration-based (using baseline and model-based reference volumes). One can observe the short tracts being detected by motion scrubbing at high corruption percentages due to the exclusion of too many gradient directions.

FIGURE 9

Figure 9. Sample tractography selection for the cingulum of the cingulate gyrus (CG) from the untilted motion-free dataset as well as selections from motion-corrected datasets with 10, 30, and 70% corrupted gradient directions. Correction choices shown include outlier-based (i.e., motion scrubbing) and registration-based (using baseline and model-based reference volumes). Notice the inability of motion scrubbing to detect an anatomically realized CG at high corrupted percentages.

FIGURE 10

Figure 10. Sample tractography selection for the corticospinal tract (CST) from the untilted motion-free dataset as well as selections from motion-corrected datasets with 10, 30, and 70% corrupted gradient directions. Correction choices shown include outlier-based (i.e., motion scrubbing) and registration-based (using baseline and model-based reference volumes). Note that motion scrubbing cannot recover long tracts such as CST beyond 10% motion corruption.

FIGURE 11

Figure 11. Sample tractography selection for the inferior fronto-occipital tract (IFO) from the untilted motion-free dataset as well as selections from motion-corrected datasets with 10, 30, and 70% corrupted gradient directions. Correction choices shown include outlier-based (i.e., motion scrubbing) and registration-based (using baseline and model-based reference volumes). Note that motion scrubbing cannot recover long tracts such as IFO beyond 10% motion corruption. Further, motion-based motion correction tends to recover longer tracts at high motion corruption compared to baseline-based correction.

FIGURE 12

Figure 12. Sample tractography selection for the uncinate fasciculus (UNC) from the untilted motion-free dataset as well as selections from motion-corrected datasets with 10, 30, and 70% corrupted gradient directions. Correction choices shown include outlier-based (i.e., motion scrubbing) and registration-based (using baseline and model-based reference volumes). Notice the inaccurate UNC tract being detected from the motion scrubbing choice at high percentages of motion corruption.

4. Discussion

In this section, we discuss the impact of different motion correction choices using local as well as global metrics.

4.1. Voxel-Based Results

Heterogeneous regions are more affected by motion correction, showing larger average JSD in general when compared to the single fiber regions, regardless of the correction mode, interpolation scheme, or reference volume employed (see Figure 3).

The impact of motion scrubbing (removing gradient directions) becomes more pronounced with more motion-corrupted directions when compared to registration-based correction (see Figure 3A). Meanwhile, the JSD values indicate minimal deformations in fODFs reconstructed for baseline-based correction at high SNR levels compared to model-based correction, whereas both choices show comparable average JSD values at low SNR levels. This complies with the conclusions presented in Ref. (20).

Forcing the correction and interpolation of all gradient directions shows comparable performance compared to the correction and interpolation of only the corrupted directions (see Figure 3A). This observation discourages the choice of heuristic parameters on motion parameters beyond which directions are claimed to be corrupted and interpolated. Further, interpolation of all directions causes less impact on the reconstructed fODFs at low corrupted percentages (<50%). We can assume, therefore, that motion is omnipresent and can be corrected for by the alignment and interpolation of all gradient directions.

On the interpolation aspect of correction, the sampling theory suggests the sinc kernel as the ideal interpolation kernel; nonetheless, this gives rise to the Gibbs phenomena (i.e., ringing) due to kernel truncation. This explains the smaller fODF deformation when using trilinear interpolation compared to sinc interpolation. Trilinear interpolation, which is much faster, is probably sufficient for motion correction.

In Figure 3B, one can observe the comparable impact of different motion correction choices at low motion corruption percentages (<30%). Whereas with higher motion corruption, a situation that is encountered in studies including infants, for example, motion scrubbing shows a significant impact on the reconstructed fODFs even at high SNR levels. This effect is more pronounced in regions with crossing fibers where the ability to resolve fiber crossings is deteriorated especially as the separation angle of the fibers decreases.

Further, baseline-based motion corrections show minimal JSD values with higher corruption levels (> 50%) when compared to model-based corrections, regardless of the interpolation scheme employed. The difference in performance between baseline-based and model-based becomes more significant as the SNR level increases.

The denoising process yields smaller JSD values for low SNR levels (<12) (see Table 1), while providing comparable performance for baseline-based and model-based motion correction choices. The slight decrease of JSD values for denoised datasets compared to the raw ones is due to the fODF reconstruction processes where we use the constrained spherical deconvolution (CSD) technique (57). In an iterative manner, the deconvolution process in CSD applies a non-negativity constraint on the estimated fODFs as negative fiber orientation densities are physically impossible. This process provides fODFs estimates that preserve the angular resolution while being robust to noise. Yet, as a word of caution, the denoising process, when applied to motion-corrupted datasets, should not take into consideration the joint information from diffusion gradients since voxel-wise correspondence between different diffusion volumes is not guaranteed.

In Figure 4, one can observe the significant impact of motion scrubbing (i.e., outlier-based correction) on the reconstructed fODFs for mildly corrupted datasets (e.g., 30% corrupted directions). Further, it can be noticed that with >50% motion corruption, model-based reconstruction infers the spatial position from the gradients corresponding to the tilted brain due to its majority (i.e., gradients of the untilted brain are considered the motion-corrupted directions).

Due to the insufficient number of gradients and unbalanced sampling of the q-space, the impact of motion scrubbing on the estimated fiber orientations becomes evident as SNR decreases and/or corrupted directions increase (see Figure 5).

Although interpolating all directions versus corrupted directions reports comparable orientation deviation with lower impact on fractionally corrupted datasets (<50%), we still favor forcing such a process to all directions to avoid the ad hoc process of thresholding motion parameters.

Nonetheless, one can notice the peaked performance of the orientation deviation at 50% corrupted directions for model-based motion correction choices. The explanation of this phenomenon is based on the fact that, with >50% of the gradients being corrupted (i.e., corresponding to the tilted brain), the formed reference volumes would instead infer its anatomical structure from the tilted brain. For highly corrupted datasets, the gradients corresponding to the untilted brains become the corrupted directions (i.e., a 70% corruption will have a performance similar to the 30% case).

Model-based corrections display higher impact on the JSD of the reconstructed fODFs at higher levels of motion corruption, but such corrections have a smaller impact on the fiber orientation deviations especially when interpolating all directions (trilinear interpolant). This change of JSD metric implies an increase in the overall fODF volume when compared to the reconstructions from the motion-free dataset, yet the fODFs maintain the voxel-wise fiber crossing structure. This observation is more pronounced for fibers with the largest fiber volume fraction.

4.2. Global Connectivity-Based Results

Whereas there is a slight performance difference between GDD values computed based on raw datasets versus those from denoised dataset, thanks to the fODF reconstruction that is robust to noise contamination, one may observe consistent findings when GDD is compared to the JSD metric. In particular, the global brain connectivity is least affected by the motion correction step when forcing the alignment and interpolation of all gradient directions without setting a predefined threshold to claim corrupted volumes. There is a significant difference between GDD values obtained from trilinear interpolation compared to sinc interpolation. This implies that the impact of sinc interpolation on the fODFs, being encoded by the JSD metric, yields global brain connectivity that is different from the “motion-free”-based brain connectivity.

Whereas the effect of motion correction is evident at higher corrupted percentages (except for motion scrubbing), one can notice the effect of noise where the impact of motion correction becomes more significant at low SNRs (<12), while different correction choices (except motion scrubbing) render slight performance difference at high SNRs (>12). Moreover, being consistent with different SNR levels, the baseline-based correction choices yield connectivity graphs with minimal deviations (smaller GDD) compared to their corresponding model-based choices.

On the contrary, motion scrubbing displays a different behavior. The GDD values from the scrubbed datasets, though maximal compared to the other correction choices, are decreasing with higher SNR levels for <50% corrupted directions, but this behavior is soon changed to the opposition direction for ≥50% corrupted directions (see Figure 6B). This change of behavior is perceivable in Figure 6A where the GDD values are maximal at 50% corruption percentage for high SNR levels (>12), whereas such a peak occurs even at low corrupted percentages (e.g., 30%) for low SNRs (<12). This phenomenon can be explained as follows: with high percentage of motion-contaminated gradients, the scrubbing (outlier-based) option tends to produce an inadequate set of gradients for accurate fODF estimation due to the exclusion of too many gradients. This unbalanced sampling of the q-space, henceforth, biases the CSD process to converge to an incorrect solution, producing inaccurate fiber orientation, and in turn imprecise brain connectivity. Hence, the increase of the GDD values with higher SNRs beyond 30% corrupted directions is due to having more short tracts connecting nearby region of interests while being assigned to larger weights in the graph construction step (see 4).

In Figure 7, one can observe the motion scrubbing behavior where the links become denser with higher corrupted percentages, implying the detection of more short tracts connecting nearby ROIs. On the other hand, the baseline-based choice reveals comparable connectograms to the motion-free ones while model-based counterpart tends to add more shorter tracts.

4.3. Tract-Based Results

Being consistent with the results from the other metrics, motion scrubbing shows a significant decrease in the degree of tract agreement when increasing the percentage of motion corruption, which in turn leads to discarding more gradient directions. With <50% corrupted directions, the tract agreement degree increases with higher SNR levels, yet such a trend changes with ≥50% where shorter or no tracts being detected, which deviates from being anatomically realistic; see, for example, the top row of Figures 8–12 where tracts can be even missed even at 70% corruption. The CST and IFO tracts are good examples of long tracts that are not recovered by motion scrubbing beyond 10% motion corruption, see Figures 10 and 11. Nonetheless, the maximal agreement is achieved when aligning and interpolating all gradient directions to correct for motion regardless of the reference volume used in the registration process (i.e., baseline versus model-based). It can be observed in Figures 8–12 that model-based motion correction is able to recover longer tracts at high corruption percentages compared to the baseline-based motion correction.

5. Conclusion: Guidelines for Motion Correction in HARDI Acquisitions

Although there is excellent theoretical work on DWI acquisition parameters and ODF reconstruction schemes, as well as their effects on the quality and crossing fiber resolution, standard users lack clear guidelines and recommendations on the best ways to approach and correct for motion in practical settings. This work investigated motion correction using transformation and interpolation of affected DWI directions versus the exclusion of subsets of DWIs, and its impact on the reconstructed fODFs, local fiber orientations, brain connectivity, and detection of fiber tracts. The various effects were systematically explored and illustrated via living phantom data, leading to the general conclusion that motion, even subtle, exists in every acquired DW scan and special care is needed to correct for motion. In the following, we summarize the findings of our analysis, which might serve as guidelines for users in practice:

– Although least recommended, motion scrubbing (removing corrupted gradient directions) can be used in studies with well-controlled environments and involving not-in-pain adults or sedated subjects, where minimal subject motion is anticipated (i.e., <10%motion corruption). Yet, this gradient removal should not result in unbalanced sampling of the q-space since the gradient distribution should be as uniform as possible on the sphere.

– Voxel-wise reconstructions, tractography, and global brain connectivity are least affected by the motion correction step when forcing the alignment and interpolation of all gradient directions without setting predefined thresholds to claim corrupted volumes.

– Using voxel-wise reconstructions that are robust to noise, the denoising process can be considered unnecessary prior to applying motion correction. Nonetheless, if applied, the denoising algorithms should not take into account joint information from different diffusion gradients since voxel-wise correspondence is not guaranteed.

– Baseline-based correction choices can be used in studies involving voxel-wise scalars, which depend on the volume of the reconstructed ODFs, especially with highly motion-corrupted datasets.

– Model-based correction choices, on the other hand, are recommended for studies requiring the recovery and analysis of long tracts, e.g., CST and IFO, especially with highly motion-corrupted datasets.

– Trilinear interpolation, although much faster compared to sinc, is probably sufficient for motion correction, where the global brain connectivity is least affected.

One may wonder that using a gold standard, which was obtained by motion correction (among other QC steps) using some of the methods under investigation could raise questions on reliability of the conclusions presented. Hence, in order to support the validity of the conclusions drawn from this study, we conducted the same set of experiments using the raw acquired data without performing any quality control. Figure 13 shows a sample result of the average JSD and local fiber orientation deviation metric for reconstructions based on gold standards generated from the QCed phantom datasets as well as the raw phantom datasets. Being consistent with the conclusions drawn from the reconstructions based on the QCed datasets, regions with crossing fibers are more affected by motion correction, showing larger average JSD in general when compared to the single fiber regions. The impact of motion scrubbing becomes more evident with more motion-corrupted directions when compared to the registration-based correction. Moreover, the peaked performance of the orientation deviation at 50% corrupted directions for model-based motion correction is also maintained. Further, forcing the interpolation of all gradients directions would have minimal impact on the reconstructions when compared to the choice of interpolating motion corrupted directions via setting a predefined threshold beyond which a direction is claimed to be corrupted.

FIGURE 13

Figure 13. The average Jensen–Shannon divergence (JSD) values (first row) and the average fiber orientation deviation (second and third row), a function of motion corrupted percentage for reconstructions based on gold standards generated from (A) the QCed phantom dataset and (B) the raw phantom dataset. Note the agreement between (A) and (B) where the impact of motion scrubbing becomes more significant with more motion-corrupted directions when compared to registration-based correction. This effect is also rendered evident for local fiber orientations.

5.1. Limitations and Future Work

The primary message of this paper is that care should be taken in deciding the processing pipeline for any DW-MRI (esp. HARDI) at hand, this involves, for example, the acquisition protocol (i.e., less redundant gradients would discourage the choice of motion scrubbing) and the participating subjects (i.e., elderly in pain, infants, unsedated subjects versus healthy adults where variable motion severity levels are anticipated). Nonetheless, the presented analysis attains some limitations, which can be outlined as follows:

• One-subject analysis: as a controlled motion experiment, we could use a scan session of subjects with repeated scans where the second shows bulk motion relative to the first one. The existing phantom data contain repeated scans taken in different sessions within 24 h and hence they have to be seen as independent scans for the same subject. As a pilot study, we therefore asked one healthy volunteer to be scanned twice in a single scan session while tilting the head between the two scans. This enables us to mix gradients between the two scans from the same subject; this cannot be done with the existing repeated independent scans. We understand that reporting our results with more than a pair of datasets (tilted and untilted brain) would support our analysis, and we will collect more scans with this experimental design in our future annual phantom scan sessions. Nonetheless, we think that this experiment, even with its limitations, contributes to establish an experimental framework that would guide the scientific community in systematically evaluating the outcomes of different preprocessing steps. In the future, we will prospectively plan to obtain more of such datasets, also including navigator shots for estimation of rotation, to extend this analysis.

• Anatomical geometric correction: echo-planner imaging (EPI) distortion, in contrast to Eddy current that affects only diffusion-weighted images, would affect all images in the acquired sequence regardless of their level of diffusion sensitization. Hence, EPI distortion correction would involve acquiring additional data for either B0 mapping or a dedicated T1 or T2-weighted structural target. That’s a primary reason behind ignoring EPI correction in most MRI processing pipelines (73). With the availability of such additional data, EPI correction would involve non-linear spatial warping that employ interpolation, a decision variable under investigation of the presented work. Hence, we favored to bypass this step in order not to inter-mingle interpolation due to motion correction and that of EPI correction. However, we think that the analysis/correction of inter-gradient spatial distortions, and its effect on ODF reconstruction, is an important issue, which we together with the scientific community need to address.

• Better gold standard generation: the living phantoms were healthy volunteers who were aware of the whole process and were keen to remain without motion. Nonetheless, the investigation of prospective navigators is a promising idea for future work to provide different types of ground truth data and to get motion estimates directly from the scanner rather than only via post-processing.

Conflict of Interest Statement

The editor and reviewer Maxime Descoteaux and the reviewer Emmanuel Caruyer declare that, despite having collaborated with author Dr. Yaniv Gur, the review process was handled objectively and no conflict of interest exists. The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Funding: this work is supported by NIH grants ACE RO1 HD 055741 and NA-MIC Roadmap U54 EB005149 and the cocaine infant project (CAMID NIDA DA022446-01). The NIH funded Autism Centers of Excellence Infant Brain Imaging Study (ACE-IBIS) Network: clinical sites: University of North Carolina: Joseph Piven (IBIS Network PI), H. C. Hazlett, C. Chappell; University of Washington: S. Dager, A. Estes, D. Shaw; Washington University: K. Botteron, R. McKinstry, J. Constantino, J. Pruett; Children’s Hospital of Philadelphia: R. Schultz, S. Paterson; University of Alberta: L. Zwaigenbaum; Data Coordinating Center: Montreal Neurological Institute: A. C. Evans, D. L. Collins, G. Bruce Pike, V. Fonov, P. Kostopoulos; Samir Das; Image Processing Core: University of Utah: G. Gerig; University of North Carolina: M. Styner; Statistical Analysis Core: University of North Carolina: H. Gu. This work was also funded in part by the NIH/NCRR Center for Integrative Biomedical Computing, P41RR12553.

Supplementary Material

The Supplementary Material for this article can be found online at http://www.frontiersin.org/Journal/10.3389/fneur.2014.00240/abstract

Footnotes

1. ^The protocol used a GRAPPA parallel imaging factor of 2 and a partial Fourier factor of 3/4, which does indeed result in non-Rician noise distributions. However, the effect of the noise distribution is expected to be relatively small at a b values of 2000 s/mm2 [e.g., Ref. (74)] and we do not expect the difference in noise profile to affect our conclusions in terms of the motion correction schemes.
2. ^In our analysis, we opt to using a prospective approach [a Twice-Refocused Spin Echo (TRSE) sequence] for eddy current compensation in order not to introduce any alignment-based preprocessing before running motion correction that could already have a confounding effect (otherwise we would work on resampled and interpolated images before motion detection). Further, we have eye-balled the FA map of the acquired sequences (prior to motion correction) where the prominent edge artifact (regions of anomalous diffusion contrast resulted from misregistration of dissimilar materials) that should be visible in case of Eddy current distortion was almost entirely absent in the TRSE images. To further support our decision, we used FSL-MCFLIRT (31) to provide the affine transformation matrix (i.e., 12 degrees of freedom) to detect the scaling and skewing parameters, which might occur due to induced eddy current where the affine transformation matrix can be written as:
   $M=R\left(\begin{array}{ccc}\hfill s_x\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill s_y\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill s_z\hfill \\ \hfill \hfill \end{array}\right)\left(\begin{array}{ccc}\hfill 1\hfill & \hfill a\hfill & \hfill b\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill c\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill \hfill \end{array}\right)+\left(\begin{array}{c}\hfill t_x\hfill \\ \hfill t_y\hfill \\ \hfill t_z\hfill \\ \hfill \hfill \end{array}\right)$
   where R is the rotation matrix, s_x, s_y, s_z ∈ ℝ are the scaling parameters, a, b, c ∈ ℝ are the skewing parameters and t_x, t_y, t_z ∈ ℝ are the translation parameters. Table S2 in Supplementary Material reports the average and standard deviation (along all gradient directions per dataset) of the estimated transformation parameters where all datasets tend to have a unit scale with minimal skewing values. These values confirm the decision of bypassing eddy current compensation in our analysis.
3. ^http://www.sci.utah.edu/software/atlaswerks.html
4. ^P denotes positive and N denotes negative

References

1. Oguz I, Farzinfar M, Matsui J, Budin F, Liu Z, Gerig G, et al. DTIPrep: quality control of diffusion-weighted images. Front Neuroinformatics (2014) 8:4. doi: 10.3389/fninf.2014.00004

CrossRef Full Text | Google Scholar

2. Jenkinson M, Beckmann CF, Behrens TE, Woolrich MW, Smith SM. FSL. Neuroimage (2012) 62(2):782–90. doi:10.1016/j.neuroimage.2011.09.015

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

3. Pierpaoli C, Walker L, Irfanoglu M, Barnett A, Basser P, Chang L, et al. Tortoise: an integrated software package for processing of diffusion MRI data. Book TORTOISE: An Integrated Software Package for Processing of Diffusion MRI Data ISMRM 18th Annual Meeting. Stockholm: (2010). 1597 p.

Google Scholar

4. Yendiki A, Koldewyn K, Kakunoori S, Kanwisher N, Fischl B. Spurious group differences due to head motion in a diffusion MRI study. Neuroimage (2014) 88:79–90. doi:10.1016/j.neuroimage.2013.11.027

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

5. Jones DK, Knösche TR, Turner R. White matter integrity, fiber count, and other fallacies: the do’s and don’ts of diffusion MRI. Neuroimage (2013) 73:239–54. doi:10.1016/j.neuroimage.2012.06.081

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

6. Wolff JJ, Gu H, Gerig G, Elison JT, Styner M, Gouttard S, et al. Differences in white matter fiber tract development present from 6 to 24 months in infants with autism. Am J Psychiatry (2012) 169(6):589–600. doi:10.1176/appi.ajp.2011.11091447

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

7. Gilmore JH, Kang C, Evans DD, Wolfe HM, Smith JK, Lieberman JA, et al. Prenatal and neonatal brain structure and white matter maturation in children at high risk for schizophrenia. Am J Psychiatry (2010) 167(9):1083–91. doi:10.1176/appi.ajp.2010.09101492

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

8. Dumas EM, van den Bogaard SJ, Ruber ME, Reilmann R, Stout JC, Craufurd D, et al. Early changes in white matter pathways of the sensorimotor cortex in premanifest Huntington’s disease. Hum Brain Mapp (2012) 33(1):203–12. doi:10.1002/hbm.21205

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

9. Le Bihan D, Poupon C, Amadon A, Lethimonnier F. Artifacts and pitfalls in diffusion MRI. J Magn Reson Imaging (2006) 24(3):478–88. doi:10.1002/jmri.20683

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

10. Tuch DS, Reese TG, Wiegell MR, Makris N, Belliveau JW, Wedeen VJ. High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity. Magn Reson Med (2002) 48(4):577–82. doi:10.1002/mrm.10268

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

11. Soares JM, Marques P, Alves V, Sousa N. A hitchhiker’s guide to diffusion tensor imaging. Front Neurosci (2013):7. doi:10.3389/fnins.2013.00031

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

12. Pierpaoli C. Artifacts in diffusion MRI. Diffusion MRI: Theory, Methods and Applications. New York, NY: Oxford University Press (2010). p. 303–18.

Google Scholar

13. Gumus K, Keating B, Poser BA, Armstrong B, Chang L, Maclaren J, et al. Prevention of motion-induced signal loss in diffusion-weighted echo-planar imaging by dynamic restoration of gradient moments. Magnetic Resonance in Medicine (2014) 71:2006–13. doi:10.1002/mrm.24857

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

14. Benner T, van der Kouwe AJ, Sorensen AG. Diffusion imaging with prospective motion correction and reacquisition. Magn Reson Med (2011) 66(1):154–67. doi:10.1002/mrm.22837

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

15. Nunes RG, Jezzard P, Clare S. Investigations on the efficiency of cardiac-gated methods for the acquisition of diffusion-weighted images. J Magn Reson (2005) 177(1):102–10. doi:10.1016/j.jmr.2005.07.005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

16. Rohde G, Barnett A, Basser P, Marenco S, Pierpaoli C. Comprehensive approach for correction of motion and distortion in diffusion-weighted MRI. Magn Reson Med (2004) 51(1):103–14. doi:10.1002/mrm.10677

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

17. Chang L-C, Walker L, Pierpaoli C. Informed restore: a method for robust estimation of diffusion tensor from low redundancy datasets in the presence of physiological noise artifacts. Magn Reson Med (2012) 68(5):1654–63. doi:10.1002/mrm.24173

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

18. Mohammadi S, Möller HE, Kugel H, Müller DK, Deppe M. Correcting eddy current and motion effects by affine whole-brain registrations: evaluation of three-dimensional distortions and comparison with slicewise correction. Magn Reson Med (2010) 64(4):1047–56. doi:10.1002/mrm.22501

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

19. Tournier J-D, Mori S, Leemans A. Diffusion tensor imaging and beyond. Magn Reson Med (2011) 65(6):1532–56. doi:10.1002/mrm.22924

CrossRef Full Text | Google Scholar

20. Sakaie KE, Lowe MJ. Quantitative assessment of motion correction for high angular resolution diffusion imaging. Magn Reson Imaging (2010) 28(2):290–6. doi:10.1016/j.mri.2009.07.004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

21. Shi X, Kholmovski EG, Kim S-E, Parker DL, Jeong E-K. Improvement of accuracy of diffusion MRI using real-time self-gated data acquisition. NMR Biomed (2009) 22(5):545–50. doi:10.1002/nbm.1368

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

22. Porter DA, Heidemann RM. High resolution diffusion-weighted imaging using readout-segmented echo-planar imaging, parallel imaging and a two-dimensional navigator-based reacquisition. Magn Reson Med (2009) 62(2):468–75. doi:10.1002/mrm.22024

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

23. Aksoy M, Forman C, Straka M, Skare S, Holdsworth S, Hornegger J, et al. Real-time optical motion correction for diffusion tensor imaging. Magn Reson Med (2011) 66(2):366–78. doi:10.1002/mrm.22787

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

24. Oakes T, Johnstone T, Ores Walsh K, Greischar L, Alexander A, Fox A, et al. Comparison of fMRI motion correction software tools. Neuroimage (2005) 28(3):529–43. doi:10.1016/j.neuroimage.2005.05.058

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

25. Herbst M, Maclaren J, Weigel M, Korvink J, Hennig J, Zaitsev M. Prospective motion correction with continuous gradient updates in diffusion weighted imaging. Magn Reson Med (2012) 67(2):326–38. doi:10.1002/mrm.23230

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

26. Kober T, Gruetter R, Krueger G. Prospective and retrospective motion correction in diffusion magnetic resonance imaging of the human brain. Neuroimage (2012) 59(1):389–98. doi:10.1016/j.neuroimage.2011.07.004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

27. Caruyer E, Aganj I, Lenglet C, Sapiro G, Deriche R. Motion detection in diffusion MRI via online ODF estimation. Int J Biomed Imaging (2013) 2013:849363. doi:10.1155/2013/849363

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

28. Feinberg DA, Setsompop K. Ultra-fast MRI of the human brain with simultaneous multi-slice imaging. J Magn Reson (2013) 229:90–100. doi:10.1016/j.jmr.2013.02.002

CrossRef Full Text | Google Scholar

29. Alexander AL, Lee JE, Lazar M, Boudos R, DuBray MB, Oakes TR, et al. Diffusion tensor imaging of the corpus callosum in autism. Neuroimage (2007) 34(1):61–73. doi:10.1016/j.neuroimage.2006.08.032

CrossRef Full Text | Google Scholar

30. Alhamud A, Tisdall MD, Hess AT, Hasan KM, Meintjes EM, van der Kouwe AJ. Volumetric navigators for real-time motion correction in diffusion tensor imaging. Magn Reson Med (2012) 68(4):1097–108. doi:10.1002/mrm.23314

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

31. Jenkinson M, Bannister P, Brady M, Smith S. Improved optimization for the robust and accurate linear registration and motion correction of brain images. Neuroimage (2002) 17(2):825–41. doi:10.1006/nimg.2002.1132

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

32. Avants BB, Epstein CL, Grossman M, Gee JC. Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain. Med Image Anal (2008) 12(1):26–41. doi:10.1016/j.media.2007.06.004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

33. Johnson H, Harris G, Williams K. Brainsfit: mutual information rigid registrations of whole-brain 3D images, using the insight toolkit. Insight J (2007). Available from: http://hdl.handle.net/1926/1291

Google Scholar

34. Netsch T, van Muiswinkel A. Quantitative evaluation of image-based distortion correction in diffusion tensor imaging. IEEE Trans Med Imaging (2004) 23(7):789–98. doi:10.1109/TMI.2004.827479

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

35. Bai Y, Alexander DC. Model-based registration to correct for motion between acquisitions in diffusion MR imaging. Biomedical Imaging: From Nano to Macro, 2008. ISBI 2008. 5th IEEE International Symposium. Paris: IEEE (2008). p. 947–50.

Google Scholar

36. Ben-Amitay S, Jones DK, Assaf Y. Motion correction and registration of high b-value diffusion weighted images. Magn Reson Med (2012) 67(6):1694–702. doi:10.1002/mrm.23186

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

37. Leemans A, Jones DK. The B-matrix must be rotated when correcting for subject motion in DTI data. Magn Reson Med (2009) 61(6):1336–49. doi:10.1002/mrm.21890

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

38. Ramani S, Thevenaz P, Unser M. Regularized interpolation for noisy images. IEEE Trans Med Imaging (2010) 29(2):543–58. doi:10.1109/TMI.2009.2038576

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

39. Veraart J, Rajan J, Peeters RR, Leemans A, Sunaert S, Sijbers J. Comprehensive framework for accurate diffusion MRI parameter estimation. Magn Reson Med (2013) 70(4):972–84. doi:10.1002/mrm.24529

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

40. Liu Z, Wang Y, Gerig G, Gouttard S, Tao R, Fletcher T, et al. Quality control of diffusion weighted images. SPIE Medical Imaging. San Jose, CA: International Society for Optics and Photonics (2010). p. 76280J–76280J.

Google Scholar

41. Chang L-C, Jones DK, Pierpaoli C. Restore: robust estimation of tensors by outlier rejection. Magn Reson Med (2005) 53(5):1088–95. doi:10.1002/mrm.20426

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

42. Pannek K, Raffelt D, Bell C, Mathias JL, Rose SE. Homor: higher order model outlier rejection for high b-value MR diffusion data. Neuroimage (2012) 63(2):835–42. doi:10.1016/j.neuroimage.2012.07.022

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

43. Tristán-Vega A, Aja-Fernández S. DWI filtering using joint information for DTI and HARDI. Med Image Anal (2010) 14(2):205–18. doi:10.1016/j.media.2009.11.001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

44. Collins DL, Zijdenbos AP, Kollokian V, Sled JG, Kabani NJ, Holmes CJ, et al. Design and construction of a realistic digital brain phantom. IEEE Trans Med Imaging (1998) 17(3):463–8. doi:10.1109/42.712135

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

45. Tristán-Vega A, Aja-Fernández S. Design and construction of a realistic DWI phantom for filtering performance assessment. Medical Image Computing and Computer-Assisted Intervention – MICCAI 2009. London: Springer (2009). p. 951–8.

Google Scholar

46. Aja-Fernández S, Niethammer M, Kubicki M, Shenton ME, Westin C-F. Restoration of DWI data using a Rician LMMSE estimator. IEEE Trans Med Imaging (2008) 27(10):1389–403. doi:10.1109/TMI.2008.920609

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

47. Tustison NJ, Avants BB, Cook PA, Zheng Y, Egan A, Yushkevich PA, et al. N4itk: improved N3 bias correction. IEEE Trans Med Imaging (2010) 29(6):1310–20. doi:10.1109/TMI.2010.2046908

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

48. Verde AR, Budin F, Berger J-B, Gupta A, Farzinfar M, Kaiser A, et al. Unc-Utah NA-MIC framework for DTI fiber tract analysis. Front Neuroinform (2014) 7:51. doi:10.3389/fninf.2013.00051

CrossRef Full Text | Google Scholar

49. SCI Institute. AtlasWerks: An Open-Source (BSD License) Software Package for Medical Image Atlas Generation. Scientific Computing and Imaging Institute (SCI) (2014). Available from: http://www.sci.utah.edu/software/atlaswerks.html

Google Scholar

50. Smith SM. Fast robust automated brain extraction. Hum Brain Mapp (2002) 17(3):143–55. doi:10.1002/hbm.10062

CrossRef Full Text | Google Scholar

51. Goodlett CB, Fletcher PT, Gilmore JH, Gerig G. Group analysis of DTI fiber tract statistics with application to neurodevelopment. Neuroimage (2009) 45(1):133–42. doi:10.1016/j.neuroimage.2008.10.060

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

52. Joshi S, Davis B, Jomier M, Gerig G. Unbiased diffeomorphic atlas construction for computational anatomy. Neuroimage (2004) 23:S151–60. doi:10.1016/j.neuroimage.2004.07.068

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

53. Alexander DC, Pierpaoli C, Basser PJ, Gee JC. Spatial transformations of diffusion tensor magnetic resonance images. IEEE Trans Med Imaging (2001) 20(11):1131–9. doi:10.1109/42.963816

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

54. Fletcher PT, Joshi S. Riemannian geometry for the statistical analysis of diffusion tensor data. Signal Processing (2007) 87(2):250–62. doi:10.1016/j.sigpro.2005.12.018

CrossRef Full Text | Google Scholar

55. Oishi K, Faria A, Jiang H, Li X, Akhter K, Zhang J, et al. Atlas-based whole brain white matter analysis using large deformation diffeomorphic metric mapping: application to normal elderly and Alzheimer’s disease participants. Neuroimage (2009) 46(2):486–99. doi:10.1016/j.neuroimage.2009.01.002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

56. Burdette JH, Durden DD, Elster AD, Yen Y-F. High b-value diffusion-weighted MRI of normal brain. J Comput Assist Tomogr (2001) 25(4):515–9. doi:10.1097/00004728-200107000-00002

CrossRef Full Text | Google Scholar

57. Tournier JD, Calamante F, Connelly A. Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution. Neuroimage (2007) 35(4):1459–72. doi:10.1016/j.neuroimage.2007.02.016

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

58. Garyfallidis E, Brett M, Amirbekian B, Rokem A, Van Der Walt S, Descoteaux M, et al. Dipy, a library for the analysis of diffusion MRI data. Front Neuroinform (2014) 8:8. doi:10.3389/fninf.2014.00008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

59. Garyfallidis E. Towards an accurate brain tractography, PhD thesis. Cambridge, UK: University of Cambridge (2012).

Google Scholar

60. Parizel P, Van Rompaey V, Van Loock R, Van Hecke W, Van Goethem J, Leemans A, et al. Influence of user-defined parameters on diffusion tensor tractography of the corticospinal tract. Neuroradiol J (2007) 20(2):139–47.

Pubmed Abstract | Pubmed Full Text | Google Scholar

61. Hammond DK, Scherrer B, Warfield SK. Cortical graph smoothing: a novel method for exploiting DWI-derived anatomical brain connectivity to improve eeg source estimation. IEEE Trans Med Imaging (2013) 32(10):1952–63. doi:10.1109/TMI.2013.2271486

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

62. Kaiser M. A tutorial in connectome analysis: topological and spatial features of brain networks. Neuroimage (2011) 57(3):892–907. doi:10.1016/j.neuroimage.2011.05.025

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

63. Rubinov M, Sporns O. Complex network measures of brain connectivity: uses and interpretations. Neuroimage (2010) 52(3):1059–69. doi:10.1016/j.neuroimage.2009.10.003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

64. de Luis-García R, Westin C-F, Alberola-López C. Geometrical constraints for robust tractography selection. Neuroimage (2013) 81:26–48. doi:10.1016/j.neuroimage.2013.04.096

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

65. Jellison BJ, Field AS, Medow J, Lazar M, Salamat MS, Alexander AL. Diffusion tensor imaging of cerebral white matter: a pictorial review of physics, fiber tract anatomy, and tumor imaging patterns. Am J Neuroradiol (2004) 25(3):356–69.

Google Scholar

66. Chiang M-C, Barysheva M, Lee AD, Madsen S, Klunder AD, Toga AW, et al. Brain fiber architecture, genetics, and intelligence: a high angular resolution diffusion imaging (hardi) study. Medical Image Computing and Computer-Assisted Intervention–MICCAI 2008. New York: Springer (2008). p. 1060–7.

Google Scholar

67. Cohen-Adad J, Descoteaux M, Wald LL. Quality assessment of high angular resolution diffusion imaging data using bootstrap on Q-ball reconstruction. J Magn Reson Imaging (2011) 33(5):1194–208. doi:10.1002/jmri.22535

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

68. Hammond DK, Gur Y, Johnson CR. Graph diffusion distance: A difference measure for weighted graphs based on the graph Laplacian exponential kernel. Global Conference on Signal and Information Processing (GlobalSIP), 2013. IEEE (2013). p. 419–22. doi:10.1109/GlobalSIP.2013.6736904

CrossRef Full Text | Google Scholar

69. Nadler B, Lafon S, Coifman R, Kevrekidis I. Diffusion maps, spectral clustering and eigenfunctions of Fokker-Planck operators. In: Weiss Y, Schölkopf B, Platt J, editors. Advances in Neural Information Processing Systems 18. Cambridge, MA: MIT Press (2005). p. 955–62.

Google Scholar

70. Landis JR, Koch GG. The measurement of observer agreement for categorical data. Biometrics (1977) 33(1):159–74. doi:10.2307/2529310

CrossRef Full Text | Google Scholar

71. Hallgren KA. Computing inter-rater reliability for observational data: an overview and tutorial. Tutor Quant Methods Psychol (2012) 8(1):23–34.

Pubmed Abstract | Pubmed Full Text | Google Scholar

72. Krzywinski M, Schein J, Birol I, Connors J, Gascoyne R, Horsman D, et al. Circos: an information aesthetic for comparative genomics. Genome Res (2009) 19(9):1639–45. doi:10.1101/gr.092759.109

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

73. Irfanoglu MO, Walker L, Sarlls J, Marenco S, Pierpaoli C. Effects of image distortions originating from susceptibility variations and concomitant fields on diffusion MRI tractography results. Neuroimage (2012) 61(1):275–88. doi:10.1016/j.neuroimage.2012.02.054

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

74. Jones DK, Basser PJ. Squashing peanuts and smashing pumpkins: how noise distorts diffusion-weighted MR data. Magn Reson Med (2004) 52(5):979–93. doi:10.1002/mrm.20283

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar