lop-DWI: A Novel Scheme for Pre-Processing of Diffusion-Weighted Images in the Gradient Direction Domain

We describe and evaluate a pre-processing method based on a periodic spiral sampling of diffusion-gradient directions for high angular resolution diffusion magnetic resonance imaging. Our pre-processing method incorporates prior knowledge about the acquired diffusion-weighted signal, facilitating noise reduction. Periodic spiral sampling of gradient direction encodings results in an acquired signal in each voxel that is pseudo-periodic with characteristics that allow separation of low-frequency signal from high frequency noise. Consequently, it enhances local reconstruction of the orientation distribution function used to define fiber tracks in the brain. Denoising with periodic spiral sampling was tested using synthetic data and in vivo human brain images. The level of improvement in signal-to-noise ratio and in the accuracy of local reconstruction of fiber tracks was significantly improved using our method.

Introduction

Diffusion-weighted imaging (DWI) is a powerful tool for inferring tissue structure and has been used extensively to map white matter pathways in healthy and diseased brains (1, 2). Errors due to thermal and physiological noise and to eddy currents affect individual diffusion-weighted images and influence the accuracy of measures obtained from DWI data such as fractional anisotropy, fiber orientation, and separation between fibers (3–5).

Averaging over a number of data acquisitions improves the signal-to-noise ratio (SNR) but is time consuming, limiting practicality. Hence, several post-processing techniques have been proposed to enhance the quality of diffusion-weighted images. Parametric and non-parametric statistical approaches have been used to describe the noise distribution and to derive the best fit for tensor parameters (6–8). Parametric methods rely on the distributional model of noise being correct. Non-parametric statistical approaches are model independent but are less powerful in hypothesis testing and more computationally demanding (7).

A number of techniques that do not require statistical analysis have been employed to enhance information contained in diffusion-weighted images (7). Image SNR can be enhanced by filtering individual images using the wavelet shrinkage technique (9). This approach is highly efficient due to the sparse representation; however, early implementations assumed spatially homogeneous variance of noise. Nowak and Pajevic described a method of wavelet-based denoising with spatially inhomogeneous noise but this approach suffered from the disadvantage of blurring edges in the image (7, 10). SNR can also be improved by exploiting spatial variations and similarities across an image, for example, using denoising techniques such as non-linear anisotropic diffusion, non-local means, or their variants (11–19). Alternatively, regularization functions can be used in the model fitting stage to increase the accuracy of parameter extraction (20, 21).

Incorporation of prior knowledge may also increase the effectiveness of denoising. In this paper, we propose a sampling and pre-processing scheme in which gradient direction encodings vary in continuous steps on a unit sphere, resulting in a periodic diffusion-weighted signal in each voxel. In this domain, which we call the “gradient direction domain,” the diffusion signal in each voxel is concentrated in the low-frequency portion of the Fourier spectrum of the acquired signal. This prior knowledge can be used for effective filtering of signal from noise using simple low-pass filtering in the Fourier domain.

To assess the impact of our method on the accuracy of diffusion-weighted data analysis, we first evaluated it on a digital phantom, where ground truth signal and fiber orientations are known. We compared mean squared error (MSE) before and after pre-processing datasets across a wide range of SNR levels. Local reconstruction accuracy was evaluated according to fiber population estimation (number of fiber bundles in each voxel) and angular estimation. For human brain data, we concatenated eight DWI acquisitions in the same subject to provide “gold standard” data similar to (22). The accuracy of the new method was assessed in terms of SNR improvement and accuracy of fiber population estimation.

Method

Pre-Processing with Low-Pass DWI

A periodic and continuous arrangement of gradient directions according to orientation results in periodicity of the diffusion-weighted signal in each voxel, which can be exploited to denoise DWI data using frequency domain techniques in a voxel-wise manner. Gradient directions used in DWI can be mapped to points on a unit sphere, which are parameterized by the longitude angle and latitude angle. Here, we term this parameterization the “gradient direction domain.” A spherical spiral curve can be defined to traverse gradient direction domain, forming a one-dimensional continuous indexing of gradient directions (Figure 1C). This one-dimensional indexing is discretized to generate a sample set of gradient directions d_n. The acquired diffusion-weighted signals in each voxel form a one-dimensional signal, s, in the gradient direction domain, which is defined as follows:

$\begin{array}{lll}\hfill \mathbf{\text{s}}& =& \left[s_1,\dots ,s_{\mathrm{N}}\right],\\ s_{\mathrm{n}}& =& E\left(b,{\mathbf{\text{d}}}_n\right)∕E_0,\\ b& =& \left(\mathrm{\Delta }-\mathrm{\delta }∕3\right){\left(\mathrm{\delta \gamma }\left|G\right|\right)}^2,\end{array}$ (1)

where E(b,d_n) is the diffusion-weighted signal when the gradient is applied in direction d_n, b (the b-value) is an example for a Stejskal–Tanner pair of diffusion gradients (23), E₀ is the signal without diffusion weighting (b = 0), Δ is the diffusion time, δ is the pulse length, γ is the gyromagnetic ratio, and G is the gradient strength. The periodic and continuous samples yield a signal s that is pseudo-periodic, continuous, and “smooth” (i.e., s_n ≈s_n+1) and that concentrates at low frequencies (Figure 1D).

FIGURE 1

Figure 1. lop-DWI: sampling and filtering. (A) Evenly distributed spiral samples over the unit sphere, (B) periodic spiral sampling over the unit sphere, sampled with 20° azimuthal angular velocity, (C) noise-free and noisy diffusion-weighted signals of a given voxel of the phantom in the gradient direction domain, (D) a representative view of diffusion-weighted signal in frequency domains, where blue area represents low-frequency part of the diffusion-weighted signal, (E) noise-free and filtered diffusion-weighted signals of a given voxel of the phantom in the gradient direction domain.

In our previous work, we outlined a spiral sampling scheme (24) as shown in Figure 1B. Such a spiral sampling scheme can be used to improve reconstruction of local extrema using techniques such as low-pass filtering but the samples are unevenly distributed over the sphere, potentially leading to orientation-dependent bias (25, 26). Uniform distribution of periodic spiral samples over the unit sphere avoids orientation-dependent bias (25). Thus, we used the sampling scheme proposed by Koay (27), which yields uniformly distributed periodic spiral samples based on an analytically exact spiral over the sphere (Figure 1A). Since the acquired DWI data using this sampling scheme retain intrinsic periodicity and “smoothness,” low-pass filtering can be applied to remove noise. To this end, s is detrended by subtracting its least squares linear fit, then low-pass filtered by keeping low-frequency components in the Fourier domain and zero filling the rest (Figure 1E). To simplify the wording, we term DWI using spiral sampling combined with subsequent low-pass filtering “low-pass DWI (lop-DWI).” We applied lop-DWI on datasets acquired with uneven- and even-spiral sampling schemes, passing approximately 25% of low-frequency components in the Fourier domain (cut-off value of 11). This cut-off value was chosen by comparing results of all possible cut-off values based on the local reconstruction accuracy, of which the details will be given in Section “Evaluation of Digital Phantom Data.”

To compare the efficacy of our method, multiple DWI data sets were generated on a digital phantom and acquired on a human volunteer using three different sampling schemes:

(1) Even sampling: eighty-two uniformly distributed samples on a hemisphere obtained using the electrostatic repulsion technique (25);

(2) Even-spiral sampling: uniformly distributed samples based on an analytically exact spiral (27) (Figure 1A). To generate a uniform distribution of samples satisfying the criterion of antipodal symmetry, 164 samples were generated on a unit sphere and 82 of the points that were located on a hemisphere were selected, and arranged as shown in Figure 1A;

(3) Uneven-spiral sampling with periodic angular steps (24). An azimuth angular step of 20° was used to generate 82 directions on a hemisphere, using an in-house script, written in MATLAB (Mathworks, Natick, MA, USA).

Generation of Digital Phantom Data

Digital phantom data were generated using b-values of 1,000 and 3,000 s/mm². To test the performance of each scheme in denoising and local reconstruction, Rician noise was added to simulate SNR ranging from 10 to 100 in steps of 10. SNR of 5, 125, and 150 were also included as extreme cases. The phantom consists of a set of 27 fiber bundles contained within a spherical domain as depicted in Figure 2, from which DWI were simulated using Phantomas1. The signal in white matter simulates both intra- and extra-axonal diffusion, using a model similar to CHARMED (28).

FIGURE 2

Figure 2. Ground truth fiber bundle geometries of the phantom, from which Rician noise-corrupted DWI were simulated.

Evaluation of Digital Phantom Data

To estimate the amount of improvement in the signal, the MSE of each dataset of different SNRs was calculated taking the ground truth signal as reference. The mean and SD of MSE across all the voxels excluding background were computed for each case.

We compared the effects of our method on the accuracy of local reconstruction of fiber orientations with regularized Q-Ball Imaging (QBI) (20) and Constrained Spherical Deconvolution (CSD) (29). QBI was chosen because it also involves spherical low-pass filtering. A regularization constant of 0.006 was used as suggested in (20, 22, 26, 30). Dipy software (31) was used to estimate the Fiber Orientation Distribution (FOD) in each voxel using a spherical harmonic of order 8 and to identify the peaks of the FODs. It was assumed that each voxel contained a maximum of three fibers. We also applied the CSD method to the phantom data, as this technique was used to evaluate human brain data. We generated phantom data with the same SNR measured from the human brain data that we collected and applied CSD technique with suggested parameter settings from literature (22, 30), on both DWI and lop-DWI data. In what follows, we call QBI and CSD reconstruction of lop-DWI data lop-QBI and lop-CSD, respectively.

Two metrics were used to compare the results of local reconstruction between methods (32):

(1) Angular error of fiber bundle reconstruction. The angular error of estimated fiber bundle orientation was defined as the angular difference (in degrees) between the ground truth peaks and the closest estimated peaks in a voxel:

$\mathrm{Angular}\mathrm{error}=\frac{\mathrm{\pi }}{180}\arccos \left(\left|{\mathit{\text{p}}}_{\mathit{gt}}\cdot {\mathit{\text{p}}}_{\mathrm{i}}\right|\right),$ (2)

where p_gt and p_i are ground truth peak and an estimated peak of interest in each voxel, respectively. For a voxel containing multiple fiber bundles, the average angular error across fiber bundles was taken as the voxel-wise angular error.

(2) Error in the estimated number of fibers. The number of peaks in the FODs was counted as the number of fibers. Estimation error was then calculated by comparing the estimated number of fibers with ground truth. The percentage of voxels with over- or under-estimation of fiber number and the success rate, defined as the percentage of voxels in which all fibers were successfully reconstructed, are reported. Statistically significant differences (p < 0.05) in the two metrics between lop-QBI data and QBI data were identified using paired t-test.

Data Acquisition in a Human Subject

The medical research ethics committee of the University of Queensland approved the study, under the guidelines of the National Health and Medical Research Council of Australia. Informed consent was obtained from the volunteer. The image dataset was anonymized. A healthy volunteer was scanned on a Siemens 3 T Trio using a 32-channel head coil over a period of 3 h to acquire 12 sets of DWI data with the following sampling schemes: eight datasets with electrostatic even sampling (25), optimally ordered as described in Ref. (33); two datasets with even-spiral sampling (Figure 1A); and two datasets with uneven-spiral sampling (Figure 1B). For each dataset, 82 diffusion-weighted images with at least three unweighted (b₀) images were acquired. A twice-refocused bipolar diffusion spin-echo sequence (34) was used with 2.5 mm isotropic resolution, TE/TR = 112/9400 ms, partial Fourier filling of 6/8, matrix size of 100 × 100 × 55, b = 3,000 s/mm² (26), iPAT = 2 with GRAPPA reconstruction. The acquisition time to obtain each dataset was ~14 min.

The eight datasets with electrostatic even sampling were concatenated (not averaged) into a single dataset, resulting in a total of 656 DWIs and 48 b0s (The eight datasets have been made available at https://sites.google.com/site/lopdwi). Images were corrected for subject motion and residual eddy current-induced geometric distortions with the required B-matrix adjustments using ExploreDTI software (35). Lop-DWI data were obtained from datasets acquired using even (“even lop-DWI”) and uneven (“uneven lop-DWI”) spiral schemes using a cut-off frequency of 11.

Evaluation of Human Brain Data

The impact of lop-DWI was assessed in the human brain data by comparing SNR and the accuracy of estimates of fiber number. The corpus callosum was selected as the white matter structure of interest. It was segmented by choosing voxels with FA >0.7 and a principal eigenvector in the x-direction. SNR was calculated (31) by dividing the mean signal in the region of interest by the SD in background (36) voxels identified using median Otsu segmentation (31). For each dataset, we measured SNR across all gradient directions and the average value is reported here. The results from the eight evenly sampled datasets were averaged and are also included in the SNR comparison.

Local reconstruction accuracy was evaluated by comparing the estimation accuracy of fiber number across white matter regions. We used the CSD technique with suggested parameter settings in (22, 30). MRtrix (30) was used to apply CSD on all of the evenly sampled DWI data and lop-DWI data. Similar to phantom data, spherical harmonic order up to degree 8 was used to obtain FODs. For each dataset, voxels with FA >0.7 were identified and the spherical harmonic decompositions of all the resulting profiles were then averaged to estimate the response function. Fibers of each voxel were extracted by applying a threshold of 0.1 to FOD peak amplitude, assuming maximum fiber population of three for each voxel.

From the peaks, we calculated the percentage of voxels in white matter with single, two, and three (or more) fibers. To obtain region-specific statistics on fiber population, the Harvard–Oxford subcortical structural atlas (37) was used to segment out the white matter and the JHU–MNI–ss atlas (38) was used to segment white matter into 45 sub-regions. All masks were then eroded by one voxel to exclude boundary voxels with uncertainty. Masks were obtained using DSIStudio software (39).

Results

Phantom Data

Signal-to-noise ratio

Voxel-level MSE comparison across different datasets showed that mean MSE is significantly decreased when lop-DWI is applied (p < 1e−10), as summarized in Figure 3. For b-value = 1,000 s/mm², MSE of even lop-DWI data with SNR = 20 were lower than that from DWI with SNR = 30. For b-value = 3,000 s/mm², MSE of even lop-DWI data with SNR = 20 was only slightly higher than that from DWI with SNR = 30. For both b-values, MSE of lop-DWI for phantom data with SNR = 30 were as low as those from DWI data with SNR = 50. Note that even for SNR value as high as 150, significant improvement in signal was obtained (p < 1e−10, paired t-test). MSE of even lop-DWI and uneven lop-DWI data were almost equal (data are not shown here).

FIGURE 3

Figure 3. Mean square error over different SNRs for digital phantoms, simulated with (A) b-value = 1,000 s/mm², and (B) b-value = 3,000 s/mm². To aid visualization of MSE values, bar charts for SNRs above 30 were expanded and shown inside each plot.

Angular accuracy

Figure 4 shows that angular accuracy was improved when using lop-QBI, particularly when SNR was low. For b-value = 1,000 s/mm², decreases of angular error were significant for SNRs smaller than or equal to 50 (p < 1e−10). This held for SNRs smaller than or equal to 70, when b-value = 3,000 s/mm² was used (p < 1e−10). For high SNRs, improvements in angular estimation were slight or unnoticeable. However, for low SNRs, such as SNR = 20, the mean angular estimation was improved by almost 10°.

FIGURE 4

Figure 4. Mean and SD of angular error between fibers of digital phantom and the ground truth, measured from raw-DWI (black) and lop-DWI (red), across different SNRs. (A) is the results from simulation with b-value = 1,000 s/mm² and (B) with b-value = 3,000 s/mm². To aid visualization, mean and SD of angular error for SNRs of 20, 30, and 40 are illustrated separately inside each plot. In (A), all the mean differences ≤50 are statistically significant. In (B), all the mean differences ≤70 are statistically significant (p < 1e−10).

Fiber population

Figure 5 shows that the achieved success rate from lop-QBI was notably higher than that from QBI. For instance, for both b-values, success rate of lop-QBI when SNR = 50 were higher than or as high as the highest QBI success rate. For b-value = 3,000 s/mm² and for SNR = 20 and 30, success rate increased almost 35%, which can be explained by the substantial decrease in the over-estimation rate. For low SNR, spurious fibers were reconstructed almost in all the voxels, which were the consequence of fitting the model to the noise. This explains low under-estimation rate when SNR was low. Over- and under-estimation rate of lop-QBI were almost equal when SNR was high. However, QBI over-estimation rates were always higher than under-estimation. Under-estimation rate never exceeded 11% across all datasets. Under-estimation rate of lop-QBI were slightly higher than QBI, and reached similar value as SNR increased. Observed under-estimation rates for high SNR came from voxels with low fiber crossing angle. It is known that high angular resolution diffusion imaging (HARDI) techniques cannot resolve low angles. We note that under-estimation rates were very low when SNR was low (<30). This is because, when SNR is too low, high fluctuation in diffusion-weighted signal leads to reconstruction of spurious fibers, causing almost no under-estimation and high over-estimation.

FIGURE 5

Figure 5. Success, over-estimation, and under-estimation rates (see Method), measured from raw-DWI (dotted lines) and lop-DWI, across different SNRs. (A) is the results from simulation with b-value = 1,000 s/mm² and (B) with b-value = 3,000 s/mm².

When these local reconstruction analyses were performed on phantom data with SNR 40 (similar to our human data) using CSD, similar pattern was observed in the estimation of number of fiber bundles (Figure 6). However, for CSD, no significant improvement in angular estimation was apparent (not shown here).

FIGURE 6

Figure 6. Mean and SD of error in the estimation of number of fibers. Over-estimation errors of lop-CSD were significantly lower than that obtained from CSD (p < 0.01). Both CSD and lop-CSD had low under-estimation values, while lop-CSD had insignificantly higher under-estimation.

Even sampling vs. uneven sampling

Figure 7 shows that evenly sampled lop-QBI data had slightly higher SNR than unevenly sampled lop-QBI data. The differences in the mean angular error between the two techniques were significant only in four datasets (see asterisks in Figure 7). It should be noted that unevenly sampled lop-DWI outperformed DWI in all the above criteria, but the improvement was slightly weaker than evenly sampled lop-DWI. Moreover, results (not shown here) from raw data acquired with evenly sampled data and with evenly sampled spiral data were almost the same but slightly different from that with unevenly sampled spiral data, which confirms the importance of even distribution of sampling to avoid orientation variance in the distribution of noise (26).

FIGURE 7

Figure 7. Mean of angular error between fibers of digital phantom and the ground truth, measured from evenly sampled lop-DWI (Figure 1A) and unevenly sampled lop-DWI (Figure 1B), across different SNRs. (A) is the results from simulation with b-value = 1,000 s/mm² and (B) with b-value = 3,000 s/mm². Asterisks show statistically significant difference in mean values (p < 0.05).

Human Data

Signal-to-noise ratio

Table 1 indicates that SNR of diffusion-weighted images was improved significantly in the evenly sampled lop-DWI results. SNR was improved almost 28% using our approach. SNR of lop-DWI for both evenly and unevenly sampled spiral data was significantly higher than raw images (p < 1e−10). In addition, the results of paired t-test indicated that SNR of lop-DWI of evenly sampled spiral data was significantly higher than SNR of all of eight DWI repeats (p < 0.05). SNR of first eight scans was significantly higher than SNR of all raw images. Surprisingly, it was not significantly higher than any of the results of lop-DWI for evenly sampled spiral data (p ≈ 0.2). It should be noted that SNR of raw images from unevenly sampled spiral data were significantly lower than other raw images. This is because in this scheme, the samples are unevenly distributed over the sphere, which leads to an uneven propagation of noise. SNR of each image varies depending on the gradient direction and the region of interest. Spiral sampling has dense acquisition along one axis and coarse acquisition along the others. Therefore, the estimated mean SNR will be weighted by the SNRs of the region with the dense acquisition. In our experiment, more data were acquired in the direction that has low SNR (x-direction).

TABLE 1

Table 1. Mean and SD of SNR across all gradient directions of raw-DWIs and lop-DWIs of human brain data, obtained from different sampling schemes.

Qualitative comparison of lop-DWI based on a specific diffusion-gradient direction and raw images after correction for motion, eddy currents, and EPI distortions shows that quality of DWI was improved (Figure 8). Note that the improvement can be appreciated in the spatial domain, although lop-DWI was applied at voxel level in the gradient direction domain.

FIGURE 8

Figure 8. Images are instances of human brain DWI from a given gradient direction, from raw-DWI (first row) and lop-DWI (second row). Raw-DWI was corrected for motion, eddy current, and EPI distortion (see Method).

Fiber Population

When CSD was used, fiber bundle counts in lop-DWI data were notably similar to that obtained from the concatenation of eight DWI datasets (Figure 9A). Percentage of voxels with three fibers (or more) was much smaller in lop-DWI than DWI, but close to the ground truth, showing low over-estimation rate with the use of lop-DWI. Our results agree with those published by Jeurissen et al. (22). It should be noted that, same values were obtained from each one of eight repeats of the DWI data, or from each one of two repeats of lop-DWI. Therefore, only one pie chart was drawn for each method.

FIGURE 9

Figure 9. Fiber population count across voxels of the white matter, (A,B) percentages of voxels with single, two, or greater than three fiber populations across the whole white matter, estimated for gold standard (see Method), raw–DWI, and lop-DWI, using (A) CSD, and (B) QBI techniques.

Across all datasets processed using QBI (Figure 9B), large percentages of voxels were estimated to have only single fiber (more than 60%), which is in contrast with the study of Jeurissen et al. (22). This suggests that CSD is more powerful in extracting information on the nature of fiber bundles in voxels.

Region-specific fiber population estimates using CSD are reported in Table 2. Estimated values indicate that lop-DWI results were considerably closer to gold standard compared with those obtained from DWI. In different regions of corpus callosum, such as Tapetum, Genu, and Splenium, for the gold standard and lop-DWI cases, majority of voxels had single fiber. However, over-estimation was observed in the majority of voxels when raw images were used. The same trend held for other regions too. For example, most of the voxels of corticospinal tracts were estimated to have two fibers for the gold standard and lop-DWI cases, while most of the voxels of the raw-DWI (third column) were estimated to have three fibers. The mean error in estimation of number of fibers using lop-CSD was significantly lower than that using CSD (Figure 10). Mean difference between gold standard and lop-CSD in voxels with single fiber was insignificant (p = 0.06). However, for voxels with multi-fiber population, errors of estimation of lop-CSD were significantly lower than DWI.

TABLE 2

Table 2. Percentages of voxels with single, two, or greater than three fiber populations across regions of white matter, estimated for gold standard (see Method), DWI, and lop-DWI (extended from Figure 9).

FIGURE 10

Figure 10. Mean and SD of error in estimation of number of fibers, using CSD (white bars) and lop-CSD (black bars), compare with estimated values from gold standard data (see Method). * and ** indicate significant levels of p < 0.01 and p < 1e−10, respectively. Paired t-test was used.

Discussion

We outlined a new pre-processing technique to enhance diffusion-weighted data. The technique is based on spiral sampling of diffusion-weighted images in gradient direction domain that results in periodicity in the acquired signal in each voxel. In the acquired signal, fiber information concentrates in the low frequencies of the Fourier spectrum, and noise in high-frequencies. The analytically exact spiral scheme of Koay (27) was used to obtain evenly distributed, yet spiral samples over the sphere. The low-frequency part of the diffusion-weighted signal, when acquired in the gradient direction domain, contains information on the orientation of underlying fibers. Therefore, filtering in this domain removes noise and increases accuracy of fiber reconstruction (i.e., a lower rate of spurious fiber reconstruction).

Using a digital phantom and in vivo human brain imaging data, we showed that the spiral sampling-based filtering scheme significantly improves SNR, increasing the accuracy of fiber reconstruction from diffusion-weighted brain imaging data. SNR obtained from human brain data filtered using our method improves the DWI quality by as much as 25% (Table 1). SNR of filtered data (~50) was close to the SNR of the data obtained from average of eight scans (~55). It should be noted that SNR improvement with the lop-DWI method was achieved in gradient direction domain, and not with spatial filtering. This may lead to over smoothing of regions with high spatial variation, i.e., blurring of edges.

Local reconstruction of lop-DWI (evenly distributed spiral sampling scheme with low-pass filtering) was more accurate compared with some conventional techniques. In particular, over-estimation rate was significantly lower with the lop-DWI approach. Figures 9 and 10 show that estimated fiber population values in different regions of white matter, obtained from lop-DWI, were very close to that obtained from our gold standard data and, to that obtained in a similar experiment by Jeurissen et al. (22). Accurate estimation of the number of fiber bundles provides additional information for interpretation of structural changes of underlying tissue and increase specificity of DTI analysis. For example, it is possible to assess if an increase or decrease in FA value is due to the loss of one of the fibers. Moreover, robust estimation of number of fiber bundles is practical for HARDI techniques that aim to quantify structural characteristics of each fiber bundle of interest separately, one of which is the apparent fiber density technique (40).

To apply a voxel-wise signal processing technique such as lop-DWI to data acquired with conventional schemes, uniform samples can be rearranged in the gradient direction domain and smoothing techniques applied. Regardless of the sampling scheme used, the diffusion-weighted signals in a given voxel form a 1D signal in the gradient direction domain. By re-arranging the 1D signal in this domain, one can benefit from prior knowledge about the signal and noise. For example, when the diffusion-weighted signals corresponding to proximal gradient directions are arranged so that neighboring directions are adjacent to each other, the 1D signal is smooth, allowing smoothing techniques such as wavelet denoising or quadratic smoothing to be applied. lop-DWI uses a periodic spiral scheme leading to a smooth and periodic signal in each voxel. This prior knowledge allows Fourier-based smoothing, which we found to be more effective than other smoothing techniques (the comparison of different denoising techniques is not reported in this manuscript).

Here, we used a digital phantom with a comprehensive representation of white matter microstructure model and additive Rician noise. While the phantom takes into consideration a variety of tissue characteristics, it does not represent the entire complexity of brain microarchitecture. For example, the phantom does not take into account spin exchange among intra-cellular and extra-cellular spaces, which occurs due to membrane permeability (41). An exact parametric representation of the highly complex structure of human brain is impossible. In this study, we used an existing model that has been previously tested in the literature to generate the phantom (28) and the widely accepted Rician noise distribution (42, 43).

For the phantom study, where ground truth signal was available, we assessed the amount of signal improvement by estimating MSE between the filtered datasets and the ground truth. Results showed that significant improvement could be obtained with lop-DWI (Figure 3). For in vivo human brain imaging, due to the large number of contributing factors, exact estimation of SNR is challenging. In DWI, SNR depends not only on the region of interest but also on the direction of the applied diffusion gradient. Therefore, we calculated SNR across all diffusion gradient directions using the same region of interest (corpus callosum) for all of the acquired datasets and reported the mean and SD. This analysis approach also provided the opportunity to test for significant changes in the means.

Regularized QBI and CSD were chosen for fiber reconstruction. QBI utilizes spherical low-pass filtering, which enables us to compare our method on one of the state-of-the-art techniques that involves low-pass filtering. CSD was chosen to replicate similar experiment (22) for human brain data to justify our gold standard values.

In the phantom, ground truth peaks of FOD were available, which enabled voxel-level assessment of local reconstruction accuracy of different techniques. For human data with the absence of ground truth, we created a gold standard by concatenating eight repeats of DWIs. Local reconstruction of the concatenated data was used as gold standard to assess the local reconstruction accuracy of each technique. In in vivo experiments, we observed that voxel-level comparison might be biased due to misalignment between datasets and motion, even after motion correction. Therefore, region-specific statistics were calculated for evaluation. These considerations for the in vivo experiment make it challenging to assess particular local reconstruction criteria, such as angular error. Therefore, a similar approach to that proposed by Jeurissen et al. (22) was adopted to assess local reconstruction accuracy of the competing techniques. Although some inferences can be made on the basis of concordance with prior knowledge or previous studies on metrics such as number of fiber bundles for in vivo imaging, more detailed analysis of white matter structure requires a gold standard obtained with techniques such as optical imaging or electron microscopy.

In the digital phantom, local reconstruction of evenly sampled spiral data outperformed unevenly sampled spiral data. In human data, for instance, we found that SNR of images obtained from unevenly sampled spiral data were lower than for the other techniques. From simulation and in vivo experiments, we found that no additional advantageous exist in using unevenly sampled spiral data (Figure 1B) compared to evenly sampled spiral data.

The influence of sample density on denoising and on the accuracy of reconstruction of the FOD function using our sampling scheme is yet to be explored fully. Too few sampling points in the diffusion gradient direction domain ultimately limit the ability to capture signal periodicity. Fourier-based low-pass filtering in gradient direction domain enables the interpretation of the signal from spirally sampled data in terms of a series of Fourier coefficients. This provides the opportunity for future studies to explore interpolation in the diffusion gradient direction domain as a method to decrease the number of samples, and consequently, the acquisition time without a corresponding trade-off in accuracy of fiber reconstruction.

Features embedded in the acquired signal shape such as the value and location of signal extrema in each voxel may contain biologically relevant information, which can be explored to obtain data-driven metrics from the lop-DWI results. Low over-estimation rates with the use of the lop-DWI method can potentially improve tractography results, i.e., less false positive tracts (44). Furthermore, filtering in the spatial domain can be combined with lop-DWI to improve the signal further. In addition, continuity in the acquisition of spiral sampling scheme provides the possibility to apply techniques for smoothing continuous data such as dynamic non-local means based denoising (45). Our scheme cannot be applied to existing evenly sampled HARDI datasets. To apply the method to DWI data already acquired using even sampling schemes, an interesting study would be to find an optimum arrangement of acquired images from evenly sampled datasets to form a pseudo-spiral scheme.

Conclusion

A new pre-processing method for DWI data based on an evenly distributed spiral sampling scheme in the gradient direction domain has been evaluated. The evenly distributed spiral sampling scheme results in voxel signals that concentrate to the low-frequency range of the Fourier spectrum. This knowledge about signal formation was exploited for denoising diffusion-weighted images at voxel level. The method was quantitatively validated on phantom and human data, and improvements in SNR and local reconstruction were observed. The evenly distributed spiral sampling scheme provides a new way to investigate patterns in sets of diffusion-weighted images at the voxel level of biological relevance.

Conflict of Interest Statement

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.

References

1. Assaf Y, Alexander DC, Jones DK, Bizzi A, Behrens TE, Clark CA, et al. The CONNECT project: combining macro- and micro-structure. Neuroimage (2013) 80:273–82. doi: 10.1016/j.neuroimage.2013.05.055

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

2. Le Bihan D, Johansen-Berg HH. Diffusion MRI at 25: exploring brain tissue structure and function. Neuroimage (2012) 61:324–41. doi:10.1016/j.neuroimage.2011.11.006

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

3. Jones DK, Cercignani M. Twenty-five pitfalls in the analysis of diffusion MRI data. NMR Biomed (2010) 23:803–20. doi:10.1002/nbm.1543

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

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

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

5. Pierpaoli CC. Artifacts in diffusion MRI. In: Jones DK, editor. Diffusion MRI: Theory, Methods and Applications. New York, NY: Oxford University Press (2011). p. 303–18.

Google Scholar

6. Jones DK, Knösche TRT, Turner RR. 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

7. Pajevic S. Statistical issues in diffusion tensor MRI. In: Jones DK, editor. Diffusion MRI: Theory, Methods and Applications. New York, NY: Oxford University Press (2011). p. 331–53.

Google Scholar

8. Pajevic S, Basser PJ. Parametric and non-parametric statistical analysis of DT-MRI data. J Magn Reson (2003) 161:1–14. doi:10.1016/S1090-7807(02)00178-7

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

9. Donoho DL, Johnstone JM. Ideal spatial adaptation by wavelet shrinkage. Biometrika (1994) 81:425–55. doi:10.1093/biomet/81.3.425

CrossRef Full Text | Google Scholar

10. Nowak RD. Wavelet-based Rician noise removal for magnetic resonance imaging. IEEE Trans Image Process (1999) 8:1408–19. doi:10.1109/83.791966

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

11. Bao L, Robini M, Liu W, Zhu Y. Structure-adaptive sparse denoising for diffusion-tensor MRI. Med Image Anal (2013) 17:442–57. doi:10.1016/j.media.2013.01.006

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

12. Coupé P, Manjón JV, Robles M, Collins DL. Adaptive multiresolution non-local means filter for three-dimensional magnetic resonance image denoising. IET Image Process (2011) 6:558–78. doi:10.1049/iet-ipr.2011.0161

CrossRef Full Text | Google Scholar

13. Coupé P, Yger P, Prima S, Hellier P, Kervrann C, Barillot C. An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images. IEEE Trans Med Imaging (2008) 27:425–41. doi:10.1109/TMI.2007.906087

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

14. Descoteaux M, Wiest-Daesslé N, Prima S, Barillot C, Deriche R. Impact of Rician adapted non-local means filtering on HARDI. In: Metaxas D, Axel L, Fichtinger G, Székely G, editors. Medical Image Computing and Computer-Assisted Intervention, Lecture Notes in Computer Science. Berlin: Springer (2008). p. 122–30. doi:10.1007/978-3-540-85990-1_15

CrossRef Full Text | Google Scholar

15. Kim Y, Thompson PM, Toga AW, Vese L, Zhan L. HARDI denoising: variational regularization of the spherical apparent diffusion coefficient sADC. Inf Process Med Imaging (2009) 21:515–27.

Pubmed Abstract | Pubmed Full Text | Google Scholar

16. Krissian K, Aja-Fernandez SS. Noise-driven anisotropic diffusion filtering of MRI. IEEE Trans Image Process (2009) 18:2265–74. doi:10.1109/TIP.2009.2025553

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

17. Liu H, Yang C, Pan N, Song E, Green R. Denoising 3D MR images by the enhanced non-local means filter for Rician noise. Magn Reson Imaging (2010) 28:1485–96. doi:10.1016/j.mri.2010.06.023

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

18. Tong M, Kim Y, Zhan L, Sapiro G, Lenglet C, Mueller BA, et al. A variational model for denoising high angular resolution diffusion imaging data. 9th IEEE International Symposium on Biomedical Imaging (ISBI). Barcelona (2012). p. 530–3. doi:10.1109/ISBI.2012.6235602

CrossRef Full Text | Google Scholar

19. Wiest-Daessle N, Prima S, Coupé P, Morrissey SP, Barillot C. Non-local means variants for denoising of diffusion-weighted and diffusion tensor MRI. In: Ayache N, Ourselin S, Maeder A, editors. Medical Image Computing and Computer-Assisted Intervention, Lecture Notes in Computer Science Berlin: Springer (2007). p. 344–51. doi:10.1007/978-3-540-75759-7_42

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

20. Descoteaux M, Angelino E, Fitzgibbons S, Deriche R. Regularized, fast, and robust analytical Q-ball imaging. Magn Reson Med (2007) 58:497–510. doi:10.1002/mrm.21277

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

21. Deriche R, Calder J, Descoteaux M. Optimal real-time Q-ball imaging using regularized Kalman filtering with incremental orientation sets. Med Image Anal (2009) 13(4):564–79. doi:10.1016/j.media.2009.05.008

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

22. Jeurissen B, Leemans A, Tournier JD, Jones DK, Sijbers J. Investigating the prevalence of complex fiber configurations in white matter tissue with diffusion magnetic resonance imaging. Hum Brain Mapp (2013) 34:2747–66. doi:10.1002/hbm.22099

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

23. Stejskal EO, Tanner JE. Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. J Chem Phys (1965) 42:288–92. doi:10.1063/1.1695690

CrossRef Full Text | Google Scholar

24. Sepehrband F, Choupan J, Vegh V, Tieng Q, Reutens DC, Yang Z. Periodic sampling in the gradient direction domain facilitates noise reduction in diffusion weighted imaging. Proceeding of 21th International Symposium on Magnetic Resonance in Medicine. Salt Lake City, UT (2013).

Google Scholar

25. Jones DK, Horsfield MA, Simmons A. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magn Reson Med (1999) 42:515–25. doi:10.1002/(SICI)1522-2594(199909)42:3<515::AID-MRM14>3.3.CO;2-H

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

26. Tournier JD, Calamante F, Connelly A. Determination of the appropriate b value and number of gradient directions for high-angular-resolution diffusion-weighted imaging. NMR Biomed (2013) 26:1775–86. doi:10.1002/nbm.3017

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

27. Koay CG. Analytically exact spiral scheme for generating uniformly distributed points on the unit sphere. J Comput Sci (2011) 2:88–91. doi:10.1016/j.jocs.2010.12.003

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

28. Assaf Y, Basser PJP. Composite hindered and restricted model of diffusion (CHARMED) MR imaging of the human brain. Neuroimage (2005) 27:48–58. doi:10.1016/j.neuroimage.2005.03.042

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

29. 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:1459–72. doi:10.1016/j.neuroimage.2007.02.016

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

30. Tournier JD, Calamante F, Connelly A. MRtrix: diffusion tractography in crossing fiber regions. Int J Imaging Syst Technol (2012) 22:53–66. doi:10.1002/ima.22005

CrossRef Full Text | Google Scholar

31. 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

32. Daducci A, Canales-Rodríguez EJ, Descoteaux M, Garyfallidis E, Gur Y, Lin YC, et al. Quantitative comparison of reconstruction methods for intra-voxel fiber recovery from diffusion MRI. IEEE Trans Med Imaging (2014) 33:384–99. doi:10.1109/TMI.2013.2285500

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

33. Cook PA, Symms M, Boulby PA, Alexander DC. Optimal acquisition orders of diffusion-weighted MRI measurements. J Magn Reson Imaging (2007) 25:1051–8. doi:10.1002/jmri.20905

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

34. Reese TG, Heid O, Weisskoff RM, Wedeen VJ. Reduction of eddy-current-induced distortion in diffusion MRI using a twice-refocused spin echo. Magn Reson Med (2003) 49:177–82. doi:10.1002/mrm.10308

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

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

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

36. Descoteaux M, Deriche R, Le Bihan D, Mangin JF, Poupon C. Multiple q-shell diffusion propagator imaging. Med Image Anal (2011) 15:603–21. doi:10.1016/j.media.2010.07.001

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

37. Desikan RS, Ségonne F, Fischl B, Quinn BT, Dickerson BC, Blacker D, et al. An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. Neuroimage (2006) 31:968–80. doi:10.1016/j.neuroimage.2006.01.021

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

38. Oishi K, Faria A, Jiang H, Li X, Akhter K, Zhang J. 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:486–99. doi:10.1016/j.neuroimage.2009.01.002

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

39. Yeh FC, Tseng WYI. NTU-90: a high angular resolution brain atlas constructed by q-space diffeomorphic reconstruction. Neuroimage (2010) 58:91–9. doi:10.1016/j.neuroimage.2011.06.021

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

40. Raffelt D, Tournier JD, Rose S, Ridgway GR, Henderson R, Crozier S, et al. Apparent fibre density: a novel measure for the analysis of diffusion-weighted magnetic resonance images. Neuroimage (2012) 59:3976–94. doi:10.1016/j.neuroimage.2011.10.045

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

41. Tanner JE. Transient diffusion in a system partitioned by permeable barriers. Application to NMR measurements with a pulsed field gradient. J Chem Phys (1978) 69:1748–54. doi:10.1063/1.436751

CrossRef Full Text | Google Scholar

42. Basu S, Fletcher T, Whitaker R. Rician noise removal in diffusion tensor MRI. In: Larsen R, Nielsen M, Sporring J, editors. Medical Image Computing and Computer-Assisted Intervention, Lecture Notes in Computer Science Berlin: Springer (2006). p. 117–25. doi:10.1007/11866565_15

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

43. Coupé P, Manjón JV, Gedamu E, Arnold D, Robles M, Collins DL. Robust Rician noise estimation for MR images. Med Image Anal (2010) 14:483–93. doi:10.1016/j.media.2010.03.001

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

44. Sepehrband F, Choupan J, Gal Y, Reutens DC, Yang Z. Periodic spherical spiral sampling in gradient direction domain improves diffusion MRI accuracy. Proceeding of HARDI Reconstruction Challenge of IEEE 10th International Symposium on Biomedical Imaging. San Francisco, CA (2013).

Google Scholar

45. Gal Y, Mehnert AJH, Bradley AP, McMahon K, Kennedy D, Crozier S. Denoising of dynamic contrast-enhanced MR images using dynamic nonlocal means. IEEE Trans Med Imaging (2010) 29:302–10. doi:10.1109/TMI.2009.2026575

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