Designer Background

Table of contents

  1. DWI Denoising with MPPCA
    1. Adaptive patching
    2. Symmetric PCA Thresholding
    3. Singular Value Shrinkage
    4. Denoising Complex Data
  2. DWI Gibbs Correction
  3. Rician bias correction
  4. EPI distortion correction and Eddy current and motion correction
  5. b0 normalization
  6. b1 inhomogeneity correction

For most single/multi-shell diffusion weighted imaging data, we recommend running designer using the following preprocessing methods. We summarize and provide references for the methods below. For their usage in Designer please see the usage page.

DWI Denoising with MPPCA

Designer includes the -denoise optional argument, which implements dMRI noise level estimation and denoising based on random matrix theory. The method exploits data redundancy in the patch-level PCA domain (Veraart2016b, CorderoGrande2019, Olesen2022). The method uses the prior knowledge that the eigenspectrum of random covariance matrices is described by the universal Marchenko-Pastur (MP) distribution.

The method used in Designer contains additional functionality compared to the original implementation. This includes adaptive patching, eigenvalue shrinkage, a modified threshold on the principal component domain, and eigenvalue shrinkage.

Adaptive patching

Adaptive patch example

In contrast to the original MPPCA approach where patches are squares or cubes around a given voxel (e.g., $5\times5\times5$ voxels), we enhanced the spatial redundancy, in addition to the redundancy in the measurements (here in the diffusion $q$-space), by pre-selecting voxels that have similar ground truth. Our goal is to minimize the number $p$ of components, and to maximize their contributions $s_i$, such that they describe most of the signal — this is the assumption of the underlying noise-free signal to be of low-rank.

Hence, for a voxel at $\br_0$, we would like to include the signals $S_m(\br_n)$ such that both $S_m$ are close to each other for different $\br_n$, and the voxels $\br_n$ are not too far from $\br_0$. The “distance” between signals can be described as: $w_{\alpha,\beta}[S(\br), S(\br’)] = |\br -\br’|^\alpha \cdot |S(\br)-S(\br’)|^\beta$.

Here, $|\br -\br’|$ is the Euclidean distance between voxels in 3 dimensions, and $|S(\br)-S(\br’)| = \sqrt{\sum_m |S_m(\br) - S_m(\br’)|^2}$ is the Euclidean distance between signals (the norm over the measurement index $m$).

The balance of preferring the distance between voxels and between signals is tuned by the choice of exponents $\alpha$ and $\beta$. Here we choose $\alpha=1$ and $\beta=2$ based on an empirical observation of improved denoising performance when $\beta > \alpha$ (preferring similarity of signals to the distance from $\br_0$). When $\alpha \gg \beta$ this method converges to the original MPPCA patching implementation (local signal-independent patch around $\br_0$). Based on the above distance function, we choose a patch around voxel at $\br_0$ as $N$ voxels (including $\br_0$), for which $w_{1,2}[S(\br_n), S(\br_0)]$ is the smallest. In Designer $N$ is 100 by default.

Symmetric PCA Thresholding

A data matrix $X = {X_{mn} } \equiv { S_m ( \br_n ) }$ of size $M\times N$ is formed by $M$ measurements from $N$ voxels in a patch. Let $M’ = \min(M,N)$ and $N’ = \max(M,N)$. Low-rank denoising corresponds to keeping $p$ largest components in the SVD of $X=\sum_{i=1}^{M’} s_i u_i\otimes v_i$, where $u_i$ and $v_i$ are the left and right singular vectors.

For the pure-noise case of $p=0$, all components of $XX^\top$ form the MP distribution, which has two independent properties: $\sigma^2 = \sum x_i /(MN)$, and $\sigma^2 = (x_i-x_{M’})/4\sqrt{MN}$. MPPCA uses these two properties to define two functions $\sigma_{1,2}(p)$ accounting for the noise variance from the bottom $M’-p$ components, with $p$ being the solution for $\sigma_1(p)=\sigma_2(p)$. Here, we use the following definitions based on considerations from Olesen et al..

\[\begin{align} \sigma_1^2(p) &= \frac{1}{(N'-p)(M'-p)}\sum_{i=p+1}^{M'} x_i \\\\ \sigma_2^2(p) &= \frac{x_{p+1}-x_{M'}}{4\sqrt{(N'-p)(M'-p)}} \end{align}\]

These are symmetric in $M’$ and $N’$, whereas the original MPPCA formulation had $N’$ instead of $N’-p$.

Singular Value Shrinkage

To overcome the eigenvalue repulsion due to noisy components, we reduce the sample singular values $s_i$ when recombining selected principal components into a low-rank matrix $\hat{X_\eta} = \sqrt{M’}\sigma \sum_{i=1}^p\, \eta(s_i/(\sqrt{M’}\sigma))\, u_i \otimes v_i$. According to Gavish et al., the optimal shrinker function based on minimizing MSE of the Frobenius norm can be expressed as

\[\eta(s) = \begin{cases} \frac{1}{s}\sqrt{(s^2-\gamma-1)^2-4\gamma} & s > 1 + \sqrt{\gamma} \,, \\ \\ 0 & s \leq 1 + \sqrt{\gamma} \,, \end{cases}\]

where $\gamma = M’/N’\leq 1$.

\[p(\lambda) = \begin{cases} \frac{\sqrt{(\lambda-\lambda_{-})(\lambda_{+}-\lambda)}}{2\pi\lambda\sigma^2(M/N)} & \lambda_{-}<\lambda<\lambda_{+} \,, \\ \\ 0 & \text{otherwise} \end{cases}\]

Denoising Complex Data

phase denoising example Designer provides users with the option to include a phase dataset as an optional input. To perform denoising on complex data, Designer first denoises the complex data using MPPCA (symmetric thresholding) using a $15 \times 15$ 2d box-shaped kernel within each slice. The denoised and spatially smoothly-varying phase $\phi_{\rm dn}(\br)$ is then unwound according to: $S_{\rm real} = $Re$\,(S_{\rm complex} e^{-i\phi_{\rm dn}})$. Finally, the noisy phase-unwound signal $S_{\rm real}$ is denoised using an adaptive 3d moving patch, symmetric thresholding, and eigenvalue shrinkage. Phase unwinding helps remove spurious components arising due to shot-to-shot phase variations in dMRI. In addition to improved denoising performance (due to the two-pass approach), MP-Complex also reduces the Rician noise floor of the data, reducing bias and increasing the precision of downstream parameter estimation.

DWI Gibbs Correction

This application attempts to remove Gibbs ringing artifacts from MRI images using the method of local subvoxel-shifts proposed by Kellner et al. and extended to Partial Fourier acquired images by Lee et al.. This method is designed to work on images acquired with either full k-space coverage or partial k-space coverage, by setting the input argument -pf.

The augmented pipline to eliminate 6/8 Partial Fouerier induced Gibbs ringing is shown here (reproduced from Lee et al.): RPG pipelines example The PF acquisition and zero filling induce extra ringings of interval $2\Delta y$ in phase-encoding (y) direction. The ringing removal pipeline has two steps: First, the weighting filter $Gy$ is applied to the Fourier domain of the original image, smoothing out ringings of interval $\Delta x$ in x-direction. The remaining ringings of interval $\Delta y$ in y direction are subsequently removed by applying unring y to the image, and ringings of interval $2\Delta y$ are removed by applying unring y to sub-images consisting of odd columns and even columns respectively. Combining the two sub-images, we obtain the first corrected image. Secondly, the weighting filter $Gx$ is applied to the Fourier domain of the original image, smoothing out ringings of interval $\Delta y$ in y direction. The remaining ringings of interval $\Delta x$ in x direction are further removed by applying unring x to the image, and ringings of interval $2 \Delta y$ in y direction are kept untreated. Then we obtain the second corrected image. Finally, the average of two corrected images yields the output image. In this figure, each image is rotated by 90$\circ$, and the column in each image is in a horizontal direction.

Phantom example of PF induced gibbs: RPG phantom example

Rician bias correction

In addition to high noise levels inherent to MRI, noise in DICOM level MRI data (data that has gone through a standard reconstruction from k-space to image space) is no longer Gaussian in profile. Since MRI data is acquired in the frequency (k-space) domain it comes with two noise channels, a real noise channel and an imaginary noise channel. DICOM level MRI data consists of the absolute value of this complex information, which means that the Gaussian noise inherent to the raw data becomes skewed during image reconstruction. The image below shows the process of how signal and noise skew in the positive direction during reconstruction. This positive bias will also bias quantitative scalar maps. Rician bias example

Designer offers two methods for the reduction of Rician bias. The first is by including the -phase option along with the -denoise option in order to perform complex denoising.

The second method is through the -rician option which uses the method of Koay et al. to approximate the level of bias using $\sigma$ estimated using MPPCA. The form of this approximate bias correction is $S’ = \sqrt{S^2 - \sigma^2}$. Where $S$ is the magnitude DWI signal.

The -rician option is designed for data with SNR > 2. If your data does not meet this criteria we strongly recommend using the -phase option for performing noise floor corrections.

EPI distortion correction and Eddy current and motion correction

Designer makes use of topup, eddy, and the mrtrix3 script dwifslpreproc, for the majority of EPI and motion correction steps. Designer encapsulates some of the surrounding image data and metadata processing steps either by parsing the BIDS json sidecar or through the use of input parameters pe_dir and -pf. It is intended to simplify these processing steps for most commonly-used DWI acquisition strategies, whilst also providing support for some more exotic acquisitions. The usage and examples sections demonstrates the ways in which the script can be used based on the (compulsory) -rpe_* command-line options, akin in the implementation of dwifslpreproc.

The script will attempt to run the CUDA version of eddy; if this does not succeed for any reason, or is not present on the system, the CPU version will be attempted instead. By default, the CUDA eddy binary found that indicates compilation against the most recent version of CUDA will be attempted; this can be over-ridden by providing a soft-link “eddy_cuda” within your path that links to the binary you wish to be executed.

Additional motion correction can be performed by using the -pre_align option, which is designed to simply align the first b0 image from each input series, such that for input series with large motion between them will be better conditioned for input to eddy.

Users can also use the command -ants_motion_correction to perform volume-to-volume rigid registrations. This option is included for cases with such severe motion or unusual gradient schemes that eddy is not enough to completely eliminate motion artifacts from the data. It should be noted that each of the options - eddy, pre_align, and -ants_motion_correction will interpolate the dwi data and these options should not all be used by default. We recommend using eddy as the default motion correction strategy and using the other two options only in cases of severe motion.

b0 normalization

Designer performs normalization in a specific manner for cases where multi-shell data is sent as separate input series. It performs voxelwise rescaling by taking the ratio of smoothed (Gaussian kernel (mm) with 3 standard deviations) b=0 images from each series to rescale all subsequent diffusion-weighted images. This is of particular importance in the case where this are changes in amplifier gain between input series, or when certain scan parameters such as “prescan normalize” are on for some series and off for others.

b1 inhomogeneity correction

Perform B1 field inhomogeneity correction for a DWI volume series using ants N4 bias field correction.