Diffusion parameter EStImation with Gibbs and NoisE Removal (DESIGNER)

abstract: We propose a new post-processing pipeline (DESIGNER) for diffusion weighted images that includes Marchenko Pastur denoising, Gibbs artifact removal, Rician correction, EPI distortion estimation, and outlier detection and thereby improves the precision and accuracy of diffusion and kurtosis tensor parameter estimation. In particular, results show few notorious black voxels on kurtosis maps while the original resolution is maintained in contrast to state-of-the-art processing methods that apply smoothing. In addition this pipeline has improved test-retest reliability compared to competing techniques.


Diffusion Kurtosis Imaging (DKI) is the intuitive extension of Diffusion Tensor Imaging (DTI), a commonly implemented clinical MR imaging technique(Jensen 2010). DTI is most commonly used to estimate fractional anisotropy (FA) and mean diffusivity (MD) within cerebral white matter of the brain. These parameters have proven to be viable biomarkers for neurodegenerative disease and white matter pathology. While the DTI model can give a broad overview of the distribution of water molecules in the brain, the model is flawed by a dependence on a Gaussian approximation of diffusion through tissue microstructure. In other words, the DTI model assumes an unrestricted flow of molecules in all directions. Cellular histology informs us that membranes and organelles create a non-Gaussian distribution of of water diffusion through tissue(Falangola 2014). This is particularly evident in the brain where the myelination of axons strongly restricts molecular diffusion. By estimating the kurtosis of the water diffusion PDF we can measure the peakedness and tail size of a Gaussian distribution with the same variance, and therefore measure the tissue properties responsible for impeding perfectly Gaussian diffusion \cite{}.

Recent studies confirm that kurtosis measures may be more specific to disease pathology than DTI measures, e.g. grading gliomas (Cauter 2012) or assessing stroke-related changes (Hui 2012). In GM, high-resolution DKI has proved to be more specific to the extracellular volume than DTI. These findings support the importance of DKI for characterizing microstructural properties of WM and GM. GM-DKI might be particularly of interest for neuroscience and clinical research. For example, amyloid plaques in the cortex as found in Alzheimer’s disease (Zhang 2012) could lead to non-Gaussian diffusion and thus be particularly efficiently detected by DKI. Biophysical models have been introduced in order to better classify the different sources of restricted diffusion in diffusion MRI (dMRI). Such models include the combined hindered and restricted model of diffusion (CHARMED) (Assaf 2004), and the fully restricted multi-compartment white matter tract integrity (WMTI) model (Fieremans 2011). These models allow for the quantification of metrics such like intra-axonal water volume fraction and intra/extra cellular diffusivity, which may be correlated with disease progression, for example, WMTI has given insight into the degeneration of white matter relative to tract myelination in Alzheimer’s disease (Benitez 2014).

The main hindrance to the widespread implementation of clinical DKI is poor conditioning of the DKI model and the requirement of expert knowledge to estimate post-tensor metrics. Signal attenuation induced by long T2-relaxations which result from the long echo-time necessary to accommodate strong gradient pulses causes the signal-to -noise ratio (SNR) of dMRI data to plummet (Veraart 2013). Higher gradient amplitudes exacerbate motion and physiological noise, eddy current artifacts, and gradient field inhomogeneities in addition to the intrinsic low SNR of raw diffusion images, causing further problems when attempting to fit the kurtosis model. Fitting the DKI model to raw dMRI data often results in physically and mathematically implausible values, for example the negative diffusion and kurtosis coefficients often observed as black voxels in some kurtosis maps (See Figure 1\label{fig1}). These physically implausible kurtosis parameters bias statistical analysis, modeling and interpretation of dMRI data. The outliers particularly originate from the combination of low radial diffusivities with thermal noise and/or Gibbs ringing (citation not found: 3). Low SNR is compensated for by increasing the slice thickness and decreasing the spatial resolution, however this compensation decreases anatomic gray matter and white matter specificity. To help solve the ill-conditioned DKI problem, images are acquired using a greater number of b-shells and in more diffusion directions.

Standardized post-processing pipelines have been built over the last few years with the purpose of further enabling the clinical feasibility of DKI. These pipelines typically rely on isotropic smoothing to reduce the effects of noise, but can cause partial volume effects and the degradation of fine anatomic detail. Even pipelines that use non-local means smoothing, smooth along low image gradients, or along specific eigenvector directions can cause the loss of GM and WM signal in dMRI datasets. Imposing positivity constraints on the DKI fit are the current practice to clean up parametric maps(Veraart 2013). Spatial smoothing inherently lowers the spatial resolution of the image and introduces additional partial volume effects that lead to complications in further quantitative analyses or to biases in microstructural modeling. Furthermore, constrained parameter fitting is time consuming and, more importantly, may bias the estimator. We advocate for a processing pipeline that targets specific sources of image artifacts and spatial noise in dMRI and corrects for each without the use of filters that degrade data quality.

The purpose of this new pipeline is to maximize the quality, precision, and accuracy of Diffusion Kurtosis parameters derived from diffusion MRI (dMRI) data without compromising efficiency. Our pipeline, dubbed DESIGNER, decomposes the problem in 2 independent steps: (a) Marchenko-Pastuer denoising, including Rician bias correction \cite{} and (b) Gibbs correction based on shifting the zero-crossings of sinc functions in k-space \cite{}. These two pre-processing steps alleviate the need for smoothing and allow for unconstrained weighted linear least squares fitting (WLLS) when performing the kurtosis tensor fit routine without lowering the resolution of the data. We demonstrate here how our new pipeline enables the robust estimation of parametric kurtosis metrics on clinically acquired DKI datasets without any additional processing time.

Materials and Methods

The simple flowchart shown in Figure 3 describes the order of operations of DESIGNER and compares it to a typical processing pipeline used in dMRI. The major difference between our DESIGNER and more commonly used pipelines is the replacement of the initial smoothing step with denoising and Gibbs artifact correction as well as unconstrained diffusion kurtosis tensor fitting.

Denoising using MP-PCA:

The first and most important step is data denoising using MP-PCAver(Veraart 2015). MP-PCA is a fast and accurate denoising technique that reduces signal fluctuations solely rooting in thermal noise, rather than from anatomical detail. This technique is enabled by exploiting data redundancy in the PCA domain using universal properties of the eigenspectrum of random (cf. thermal noise) uncorrelated covariance matrices. It is important that this step is first as it relies on noise being uncorrelated both spatially and in time. Performing this step after steps that use interpolation to reconstruct images will result in correlated noise and the failure of the denoising computation.

Figure 5 shows an example dataset of the Corpus Callosum of the human brain at \(b=1000\) s/mm\(^2\) where the image has been denoised. The leftmost image is raw data, the middle image has been denoised, and the rightmost image is the difference between the previous two. It is visible in this residual map that MP-PCA denoising removes random guassian distributed noise without compromising signal resulting from anatomic detail.

Gibbs Artifact Correction:

The MR image is reconstructed from k-space which is a finite sampling of the signal subjected to inverse Fourier transform in order to obtain the final image. At high-contrast boundaries the Fourier transform corresponds to an infinite number of frequencies and since sampling is finite a discrepancy appears in the image in the form of a series of periodic signal occilations (Veraart 2015a). These can appear in both phase-encode and frequency-encoding directions.

The origin of the MRI data matrix can be seen in the truncation of k-space during MRI data-acquisition. Consequently, correction techniques like Gegenbauer reconstruction or extrapolation methods aim at recovering these missing data. DESIGNER uses the fact that truncation in k-space is equivalent to convolution with a sinc-function in image space (Kellner 2015). Hence, the severity of the artifacts depends on how the sinc-function is sampled. Images are reinterpolated based on local, subvoxel shifts to sample the ringing pattern at the zero-crossings of the oscillating sinc-function. With this, the artifact can effectively and robustly be removed.

Rician Bias Correction

Unbiased noise variance estimation enabled by MP-PCA denoising allows for accurate MR signal intensity estimation by removing the noncentral - \(\chi\) distributed noise bias using the methods of (Gudbjartsson 1995) or the ́method of moments described in (Koay 2006) for biased estimates of noise variance. We encourage readers to derive unbiased estimates of noise variance \(\sigma^2\) by truncating diffusion datasets to b \(\leq 1000\) s/mm\(^2\) and then estimating the rician corrected signal \(A\) using \(A = \sqrt{M^2 - \sigma^2}\) where \(M\) is the raw image intensity.