Training a Neural Network for Gibbs and Noise Removal in Diffusion MRI

TL;DR

This paper presents a neural network approach for removing Gibbs artifacts and noise in diffusion MRI, improving image quality and artifact correction flexibility in clinical settings.

Contribution

Introduces CNN models for artifact removal in diffusion MRI, including complex image processing, trained on synthetic data, enhancing artifact mitigation capabilities.

Findings

CNNs effectively reduce artifacts in diffusion-weighted images

Complex image CNN reduces artifacts in partial Fourier acquisitions

Method is applicable slice-by-slice for clinical use

Abstract

We develop and evaluate a neural network-based method for Gibbs artifact and noise removal. A convolutional neural network (CNN) was designed for artifact removal in diffusion-weighted imaging data. Two implementations were considered: one for magnitude images and one for complex images. Both models were based on the same encoder-decoder structure and were trained by simulating MRI acquisitions on synthetic non-MRI images. Both machine learning methods were able to mitigate artifacts in diffusion-weighted images and diffusion parameter maps. The CNN for complex images was also able to reduce artifacts in partial Fourier acquisitions. The proposed CNNs extend the ability of artifact correction in diffusion MRI. The machine learning method described here can be applied on each imaging slice independently, allowing it to be used flexibly in clinical applications.

Equations

$$\hat{\theta }=\underset{\theta }{\mathrm{argmin}}\sum _{d=1}^D{\Vert {\mathbf{m}}_d-f_{\theta }\left({\mathbf{x}}_d\right)\Vert }_2^2,$$

$$l(s)=c+{\sum }_{i=1}^3\frac{{\alpha }_i}{1+\text{exp}\left(\frac{s-{\beta }_i}{{\gamma }_i}\right)},$$

$$h_t(\mathbf{f})=\sqrt{\frac{\mathcal{S}\left(f_{\hat{\theta }}\left({\mathbf{x}}_t\right)\right)}{\mathcal{S}\left({\mathbf{m}}_t\right)}}.$$

$$p_{\text{sim}}(\mathbf{r})=\text{exp}\left(i\sum _{s=1}^S\sum _{b=1}^{B_s}{a}_{b,s}{e}^{-\frac{1}{z_s}{\Vert \mathbf{r}-{\mathbf{r}}_{0,b,s}\Vert }_2^2}\right),$$

$$S$$

$$B_s$$

$$z_s$$

$$a_{b,s}$$

$$r_{0,b,s,j}$$

Full text

Methods: A convolutional neural network (CNN) was designed for artifact removal in diffusion-weighted imaging data. Two implementations were considered: one for magnitude images and one for complex images. Both models were based on the same encoder-decoder structure and were trained by simulating MRI acquisitions on synthetic non-MRI images.

Results: Both machine learning methods were able to mitigate artifacts in diffusion-weighted images and diffusion parameter maps. The CNN for complex images was also able to reduce artifacts in partial Fourier acquisitions.

Conclusion: The proposed CNNs extend the ability of artifact correction in diffusion MRI. The machine learning method described here can be applied on each imaging slice independently, allowing it to be used flexibly in clinical applications.

1 Introduction

The Gibbs phenomenon [1, 2] and noise are artifacts that affect all magnetic resonance imaging scans. The traditional solution for both low signal-to-noise ratio (SNR) and Gibbs artifacts has been to apply a smoothing filter to the image. Although smoothing can reduce the variance of noise, it cannot directly correct for the Rician-biased signal present in low-SNR regimes, nor can it accurately reproduce images that are Gibbs-artifact-free. Methods using Gegenbauer polynomials have been shown to be able to detect edges and accurately estimate images up to the edge for the purpose of Gibbs removal [3]. Subvoxel shifting is another method that can substantially reduce Gibbs ringing in acquisitions without partial Fourier encoding [4]. Model-based reconstruction methods [5, 6, 7, 8] can improve SNR and/or reduce Gibbs ringing. In parametric protocols with repetitive acquistions such as diffusion MRI, correlations across the diffusion directions can be used for denoising. Recently, a random matrix theory-based MP-PCA method has been developed that takes advantage of these correlations for noise estimation [9] and removal [10]. Also for diffusion, the DESIGNER pipeline [11] based on MP-PCA and subvoxel shifting has been demonstrated and validated.

A drawback of many of these methods is that they rely on models that are not fully general. MP-PCA requires the distribution of the noise to be unaltered by processing, as well as many diffusion directions in order to take advantage of the corresponding correlations [9, 10], but many clinical protocols only use a small number of diffusion directions. Meanwhile, subvoxel shifting does not perform well on partial Fourier acquisitions [4]. Partial Fourier imaging is a complex case due to the interaction between blurring effects from the undersampling and ringing effects from the Fourier series approximation of the image. The combination of these effects is difficult to model explicitly due to the presence of ringing artifacts of varying frequency and phase [4].

Machine learning offers a framework for building models that can be applied to simultaneously remove a variety of artifacts in an acquisition-flexible manner. One potential class of models suitable for these tasks is convolutional neural networks (CNNs) [12, 13]. CNNs are models composed of simple building blocks — often convolution kernels only 3 pixels wide followed by an activation function. When many convolution kernels are used together, they can approximate highly nonlinear functions of many variables. The use of simple components makes training such models ideally suited for parallelization and efficient implementation on modern computer architectures. Several recent papers have demonstrated the applicability of CNNs for medical imaging in areas including image reconstruction [14, 15, 16, 17, 18], image quality transfer [19], super-resolution [20], segmentation [21], and artifact correction [22].

Diffusion MRI is an intriguing test application for Gibbs and noise removal models. Both effects can have drastic effects on diffusion parameter maps. The calculation of higher order diffusion parameters can require the application of b-values up to $b=2500\,$s/mm2 [23], which, coupled with strong diffusion gradients and long echo times, leads to low signal-to-noise ratio (SNR) [24]. The well-known $\leq 9$% intensity variation due to Gibbs ringing [1, 2] may translate into $\sim 100$% errors in the estimated signal moments or cumulants, in particular leading to the so-called “black voxels” (masked outliers) in kurtosis maps [6]. Thus, systemic defects — such as Gibbs or effects from machine learning-based processing — are substantially amplified in the corresponding parameter maps.

Here, our aim is to develop and evaluate a machine learning approach for Gibbs and noise removal. We show that the framework of CNNs allows the simultaneous removal of these artifact sources in an acquisition-flexible manner — specifically, we extend the ability for Gibbs artifact removal to the partial Fourier setting and noise removal to the setting with few diffusion directions. We accomplish this by training the model entirely on simulations of MRI acquisitions with images from the ImageNet data set [25]. Prior to the evaluation stage of our experiments, networks trained this way have never seen an MR image, mitigating overfitting risk. Our results suggest that, with CNN processing, high-quality parameter maps can be calculated from abbreviated acquisitions with a little more than half of the encoding time required for traditional full acquisitions, opening up new possibilities for MRI. The method can be applied independently to each imaging slice, allowing application to clinical protocols with few repetitive scans. Finally, the simulation framework can be adapted to other acquistion protocols, allowing generalization to other applications.

2 Methods

In MRI, images are reconstructed by sampling a Fourier series at a discrete set of locations where the Fourier series coefficients are the acquired k-space points. This suffers from the Fourier series approximation error that manifests as the Gibbs phenomenon [1, 2]. In practice, the MRI scanner does not produce the true, noise-free Fourier series coefficients, but a noise-corrupted version of them. Here, we denote the true image at a given resolution as $\mathbf{m}_{d}$ and its noisy, Gibbs-corrupted version as $\mathbf{x}_{d}$ where $d$ is an integer that indicates the $d$th image. Our goal is to learn a functional, $f_{\theta}(\cdot)$, that estimates $\mathbf{m}_{d}$ from $\mathbf{x}_{d}$ where $\theta$ are the parameters of the functional. We accomplish this by solving the following optimization problem:

[TABLE]

where $D$ is the size of the training data set.

Solving (1) requires a data set of image pairs $(\mathbf{m}_{d},\mathbf{x}_{d}),d\in[1,...,D]$ or a means of generating them. We choose to simulate both $\mathbf{m}_{d}$ and $\mathbf{x}_{d}$ from a latent, high-resolution image that comes from the ImageNet data set [25].

Figure 1 shows our simulation pipeline. We describe the details of the steps in Section 2.2.1. Section 2.1 describes the model architecture for $f_{\theta}(\cdot)$. In Section 2.3.2 we describe test experiments that were conducted on data sets completely separate from that used for training, including real-world diffusion validation experiments in which we test the performance of two different networks and compare to state-of-the-art methods across partial Fourier factors. We show results from these experiments in Section 3.

2.1 Model Architecture

There are many potential model architectures for accomplishing deGibbsing. Architectures such as variational [14] and cascaded [15] networks offer the most flexibility for explicitly modeling system physics. However, these models often require access to the raw data and system parameters, e.g., for sensitivity map estimation. Use of such models may require accessing proprietary information or signing research agreements with vendors. Furthermore, the use of models that rely on raw data can incur substantial data storage requirements. Finally, they preclude the use of the model on previously-acquired data.

Due to these considerations, we address the Gibbs phenomenon with an image-to-image architecture. Our model comes from the class of U-Net [26] architectures. Figure 2 shows a diagram of the architecture. At the midpoint of the U-Net (shown as the bottom of the U in Figure 2), convolution operations affect almost the entire image, enabling the network to learn global features such as conjugate k-space symmetry.

We process complex images by inputting the real and imaginary components into the network as two separate channels [15]. After each convolution layer we apply batch normalization [27]. After batch normalization, we apply rectified linear unit (ReLU) activation functions, as is standard practice [28]. Unlike the standard U-Net, our model includes wavelet transforms for downsampling and upsampling. The theoretical properties of these modifications were analyzed previously [16] and demonstrated to better-preserve high-resolution features. Based on the intuition that it is easier to learn residuals rather than signals [29], we use a skipped connection over the entire architecture for the real and imaginary channels. For cases with complex inputs, we apply the magnitude operation as a final step after the skipped connection.

2.2 Training

Code for training our models is available on Github at <link>. We trained the network by simulating MRI acquisitions on photographs from the ImageNet data base [25], using resized, noise-free versions of the original photographs as target images $\mathbf{m}_{d}$. The data set contained 1,281,167 images in the training set and 50,000 images in the validation set. The “best” model was selected based on the minimal loss over the validation data set during the training epochs. Use of ImageNet for training requires specification of 1) an encoding simulation pipeline (Section 2.2.1) and 2) an optimization procedure (Section 2.2.2).

Although the images in ImageNet are not true MR images, we take this approach due to the fact that it is impossible to acquire large volumes of high-resolution diffusion images for the purpose of simulating Gibbs artifacts. An alternative would be to train from high-resolution anatomical data, but diffusion-weighted images can have substantially different features from their high-resolution anatomical counterparts. Thus, training a network on simulated Gibbs artifacts from high-resolution anatomical data runs the risk of having the network reconstruct signals that would not be present in the corresponding diffusion data. Furthermore, an anatomical data set would be highly homogeneous, limiting the generalizability of the resulting network. Conversely, ImageNet is a highly heterogeneous data set, enabling the network to learn more generic features and potentially enabling application of this approach in other MRI modalities that suffer from Gibbs artifacts. A key determinant of the feasibility of this approach is the simulation of encoding operations as they occur in vivo, which we discuss below.

2.2.1 Simulation

The simplest Gibbs simulation could be accomplished by cropping the k-space representation of a high-resolution image. However, in practice MR acquisitions are affected by many factors other than the Gibbs phenomenon. We simulate a subset of these effects. Figure 1 shows a flowchart of the simulation pipeline. The simulation begins with raw grayscale ImageNet image (converted via BT.601 color standards), resized to a standard size of $256\times 256$ using bilinear interpolation [30]. The simulation then randomly flips and/or transposes the image, each with a 50% probability. After flipping and transposing, the simulation proceeds to the random phase step. The random phase simulation generates a random set of Gaussian radial basis functions for each simulated image - this is described in detail in Appendix A. An example of simulated phase compared to a true in vivo phase case is shown in Figure 3.

The random phase simulation requires a number of probabilistic parameters for specifying the distribution of the simulated phase - these are reported in the appendix.

Since most images have a no-signal background, following random phase simulation we apply random ellipsoid cropping to the images. We also include a 10% probability that no ellipsoid cropping occurs to handle cases that may not have a no-signal background. After image-domain cropping, the processing pipelines split. For the data, $\mathbf{x}_{d}$, we apply the FFT and crop the image from a $256\times 256$ grid to a $100\times 100$ grid in k-space to simulate the Fourier series truncation effect. The simulation adds noise at levels from 1 to 32 times the mean absolute k-space value, varying the noise level within this range on a uniform base-2 logarithmic scale (average noise ratio of 10).

After the inverse FFT, both processing pipelines are merged into an image normalization layer where both images are divided by the maximum absolute value in $\mathbf{x}_{d}$. The absolute value operation is then applied to $\mathbf{m}_{d}$. To simulate an artifact-free target, $\mathbf{m}_{d}$, we interpolate $\mathbf{m}_{d}$ from a $256\times 256$ grid to a $100\times 100$ grid using cubic splines. Finally, to provide the complex-valued CNN with more information, we concatenate a standard partial Fourier reconstruction method [31] to the complex-valued $\mathbf{x}_{d}$ output in the channel dimension prior to input to the neural network.

When magnitude images are used, the inputs have Rician noise rather than Gaussian. We maintain most of the above processing stream aside from three modifications. 1) The absolute value operation is applied after inverse FFT in the $\mathbf{x}_{d}$ stream, and 2) followed by dividing by the maximum value of the $\mathbf{x}_{d}$ image for both the $\mathbf{m}_{d}$ and $\mathbf{x}_{d}$ final images. This implicitly builds Rician bias correction into the network training.

2.2.2 Optimization

Our optimization routine followed standard practice. We used the ImageNet data set [25] with standard training and validation splits. Our convolution weights were initialized with a uniform distribution [32]. ReLUs were initialized with zero parameters. We used the Adam algorithm [33] with a learning rate of $1\times 10^{-3}$ over 10 epochs to minimize the mean-squared error cost function in (1). At each epoch, we calculated the loss over the entire validation data set and saved the model with the best validation loss for the in vivo test experiments. Each model was trained independently on each partial Fourier factor. We implemented our models in the PyTorch deep learning framework [30]. Training was performed on an IBM GPU computing cluster with Power9 8335-GTH compute nodes. Each node in the compute cluster had four Nvidia V100 GPUs with 16 GB of memory. Our specific trainings utilized one GPU with a batch size of 55. Training time took approximately 10 days.

2.3 Testing

In the test phase we evaluated the proposed CNN methods for their ability to process diffusion-weighted images (DWIs) on different data sets. In all our test experiments, we define “SNR” to be the mean absolute image value divided by the complex Gaussian standard deviation. Note that this is different from the k-space ratio used during training. We conducted canonical signal processing experiments and in vivo experiments. To provide context for the performance of the CNN methods, we also processed the DWIs with a state-of-the-art (SoA) method [11], as well as a method designed for partial Fourier imaging [31]. In all cases with diffusion parameter maps, the raw DWIs were first processed by the designated method prior to parameter map estimation. To clearly show that our test data sets are distinct from our training data sets, we index test signals with $t$ instead of $d$.

Raw: These parameter maps were calculated from the raw diffusion-weighted images with no extra processing before the parameter estimation stage.

State-of-the-Art (SoA) [11]: Prior to parameter map estimation, the images were denoised using random matrix theory-based MP-PCA method [10] using a 5 x 5 x 5 voxel kernel extent with Rician bias correction [10], followed by Gibbs removal via subvoxel shifting [4] where each voxel was discretized into 20 subvoxel elements. Unlike the CNN methods, the SoA method is able to use information from low-$b$ diffusion-weighted images to denoise images from high-$b$ diffusion gradients.

Magnitude-input CNN (MCNN): The Raw diffusion images after the magnitude operation, $\mathbf{x}_{t}$, are processed by a neural network trained with the simulation pipeline in Figure 1 with the target image $\mathbf{m}_{t}$ being a magnitude, deGibbsed image without PF. This particular network attempts to restore the missing partial Fourier information. This network has 31,646,338 parameters.

Standard Partial Fourier (Standard PF): The Raw diffusion images are processed by using a standard partial Fourier reconstruction method that utilizes the phase maps from the symmetric region of k-space [31] prior to parameter estimation.

Complex-input CNN (CCNN): The Raw complex-valued diffusion images, $\mathbf{x}_{t}$, are processed by a neural network trained with the simulation pipeline in Figure 1 with the target image $\mathbf{m}_{t}$ being a magnitude, deGibbsed image without PF. This particular network attempts to restore the missing partial Fourier information. This network has 31,646,338 parameters.

2.3.1 Canonical Experiments

To characterize the deGibbsing and denoising performance of the network, we performed a set of canonical signal processing experiments with two data sets. The first data set consisted of a simple edge phantom, while the second data set consisted of the test split of the 2013 ImageNet challenge [34].

Edge-spread function: The edge phantom, meant to study the response of the networks to an ideal unit-step, consisted of a 2D digital image with half of the image filled with 1s and the other half filled with 0s. This image was generated at a 1024 $\times$ 1024 matrix size (cropped to a 100 $\times$ 100 matrix size during simulation). The location of the edge was varied by rotating it. Figure 6 includes a composite image of the phantom generated from $f_{\hat{\theta}}(\mathbf{x}_{t})$ at many different SNR levels. We performed two sets of experiments with the edge phantom: contrast-to-noise (CNR) and deGibbsing experiments. In a previous study [4], it was noted that deGibbsing performance could vary depending on the angle of the edge. To assess the performance of the method on removing Gibbs artifacts at different angles, we generated the phantom at a set of angles beetween 0 and 45 degrees at a 1024 SNR level and computed the CCNN deGibbsed image.

We also performed contrast-to-noise ratio (CNR) experiments. The goal of the CNR experiments was to quantify the full-width-half-maximum (FWHM) of the network’s line-spread function (analogous to a point-spread function in 1D) across a variety of noise levels. We accomplished this by applying the simulation pipeline in Figure 1 without random flipping and transposing at a variety of CNRs between 0 and 10, resulting in images at 100 $\times$ 100 matrix sizes. For each CNR, we applied both the CCNN model (generating $f_{\hat{\theta}}(\mathbf{x}_{t})$ realizations for each one). This resulted in a series of edge-spread function images. To compensate for noise in calculating the FWHM values, for each row of the resulting image, we fit a summation of logistic functions [35]:

[TABLE]

where $c$, $\alpha_{i}$, $\beta_{i}$, and $\gamma_{i}$, $i\in[1,2,3]$ are fitting parameters. We computed the derivative of this function and used it to calculate a line-spread function for each row of the image. From the line-spread functions, we calculated the FWHMs and averaged across the rows. For each CNR level, we repeated and averaged this experiment 50 times.

ImageNet Test Data: The goal of the tests on the ImageNet data was to assess the network’s performance across a broad variety of signals. For this we used the test split from the 2013 ImageNet challenge [34]. This data was not used at all in training and is distinct from the validation data set. The test set included 40,152 images. We applied the simulation pipeline in Figure 1 and computed the power spectral ratio on the resulting images [6]. Let $\mathcal{S}(\cdot)$ be the operator that computes the power spectral density of the image. For each image-target pair in the validation data set, $(\mathbf{x}_{t},\mathbf{m}_{t})$, we computed

[TABLE]

Then, we averaged over all $t\in[1,...,n_{\textit{test}}]$ to compute a mean $\bar{h}(\bm{f})$ for each partial Fourier factor and SNR level and computed profiles across $\bm{f}$ in the partial Fourier direction to assess resolution performance.

2.3.2 In Vivo Test Experiments

We tested the method in vivo by scanning a volunteer with a diffusion imaging protocol under IRB approval, saving the raw data. We scanned the volunteer after all model and training parameters were complete, so that these data qualify as a true prospective test set. The volunteer was scanned on a 3T Prisma scanner with a 16-channel head coil (Siemens Healthineers, Erlangen, Germany). Diffusion-weighted images were acquired with a spin echo EPI sequence at 100 $\times$ 100 matrix size over 38 slices (2 mm isotropic resolution) and along 66 isotropically-distributed directions at $b$-values of 0, 1,000, and 1,500 s/mm2. We applied the diffusion-weighted acquisition four times: once for each PF factor of 5/8ths, 6/8ths, 7/8ths, and fully-sampled, with the intention that the fully-sampled DWIs with current state-of-the-art processing serve as a gold standard. The acquisition parameters were set based on the shortest values we could achieve for the fully-sampled acquisition, leading to a TR of 9 s and a TE of 105 ms. The partial Fourier acquisitions used a standard protocol.

The images were reconstructed with adaptive coil combination [36]. Parameter map estimation was performed on the processed image using a constrained weighted linear least squares approach [37]. Since the complex CNN can take into account smooth variations in the phase of the image to recover high-resolution information, we also included as a comparison a reconstruction using a standard partial Fourier method that estimates smooth phase from the symmetric region of k-space [31] applied with our own custom software to the raw data.

3 Results

3.1 Canonical Experimental Results

Figure 4 shows the performance of the network in a canonical deGibbsing experiment with a hard edge with the edge oriented at 0 and 45 degrees.

Both the CCNN and the MCNN methods performed similarly at both angles, with the main difference being the sharpness of the edges with the MCNN method both at the simulated edge and the image boundaries. These image boundary residual artifacts do not appear in the $45\degree$ plot since this profile does not reach the image boundaries.

We also examined the CCNN activations (i.e., images from each channel after a CNN layer) in the case of the edge phantom, which are shown in Figure 5.

Figure 5 shows one channel (out of 64) from the activations after three convolutional steps and prior to the first wavelet transform. In early layers, we found that the networks tended to enhance Gibbs artifacts for removal in the final layers. Other features of the image are detected in the other network channels, as shown in Figure 1 of the supplementary material. Some channels detect changes in the phase of the image, while other channels detect flat regions. The overall properties of the network depend on integrations of all of these channels. We also examined activations in other layers after wavelet downsampling operations, but the activations in these layers were less interpretable.

Figure 6 shows the response of the CCNN in the edge experiments across a range of contrast-to-noise ratios.

Above a CNR of 1, the CNN FWHM rapidly converges towards 1. Below a CNR ratio of 1, the network’s FWHM increases as the CNR approaches 0. Figure 6 also shows a composite image of CNN image estimates, $f_{\hat{\theta}}\left(\mathbf{x}_{t}\right)$, across CNR factors. Figure 6 shows visually how as the CNR decreases, the edge begins to blur.

Figure 7 shows profiles of the mean spectral response, $\bar{h}(\bm{f})$, in the partial Fourier direction of the CCNN model across different partial Fourier factors and SNR levels over the entire test set.

The spectral responses were calculated from the 2013 test split of ImageNet as described in Section 2.3.1. As the SNR decreases, the network applies more smoothing, resulting in a loss of high-resolution features of about 30% at an SNR value of 1. The MCNN method exhibited similar profiles, but with more smoothing at low SNR - the MCNN smoothing without partial Fourier at an SNR of 1 led to about 50% attenuation of the high frequencies. In all cases, we observed features in the spectral responses around spectral cutoff points - this included the ripples at the edge of the sampled k-space in the No PF setting, with further cusps developing near the PF cutoff points in the PF settings.

3.2 In Vivo Test Experiment Results

Figure 8 shows examples of non-diffusion-weighted images before and after processing. The raw images (Raw) served as the input for the magnitude deep learning (MCNN) and complex deep learning (CCNN) methods.

Both methods remove artifacts, but the MCNN method allows residual rippling artifacts to pass through in the presence of partial Fourier. These rippling artifacts are not present in the CCNN method.

Figure 9 shows parameter maps for the different methods without partial Fourier (No PF) and with 5/8ths partial Fourier (5/8 PF). Mean diffusivities calculated from the raw DWI data (Raw) show notable noise and Gibbs ringing artifacts, while this is substantially removed with the state-of-the-art (SoA) method. However, the state-of-the-art method begins to lose its ability to compensate for the Gibbs ringing and resolution loss when partial Fourier is utilized in the acquisition. The effects of partial Fourier increases are primarily evident in the enlargement of the lateral ventricles and the presence of black lines in the vicinity of the lateral ventricles. The MCNN model is able to compensate somewhat for the ringing effects, but begins to introduce substantial artifacts at the 5/8ths partial Fourier factor, whereas the CCNN model continues to give high-quality mean parameter maps across all PF factors.

Similar trends are observed in the other diffusion parameter maps.

Figure 10 compares mean diffusivity maps across various partial Fourier factors between the CCNN and state-of-the-art methods.

The methods perform similarly without partial Fourier acceleration, but as partial Fourier acceleration increases, the image is continually degraded in the SoA method, with "black voxels" appearing around the lateral ventricles. The CCNN method mitigates the appearance of these artifacts in the parameter maps.

4 Discussion

Our results indicate that with carefully-crafted simulations, it is possible to train CNN models based solely on MR encoding operators that can be applied to in vivo MR imaging data. We characterized the performance of this method with canonical experiments, showing Gibbs removal capability with hard edges, as well as providing preliminary evidence that the network’s denoising capability depends on signal-adaptive smoothing, and that this smoothing largely takes place at contrast-to-noise ratios below 1. In principle there are no restrictions on applying the method in new applications other than retraining based on the application image matrix size. Our chosen in vivo setting here was diffusion parameter mapping due to the relevance of the artifacts and the sensitivity of diffusion parameter mapping to Gibbs effects. Nonetheless, we would recommend conducting independent performance studies for each application prior to using the method.

We found that it is possible to recover the high-resolution information missing in partial Fourier acquisitions while simultaneously reducing Gibbs artifacts and suppressing noise in a way that is robust enough for diffusion parameter estimation. Figures 6 and 7 suggest that the primary determinant of resolution with the proposed models is SNR rather than PF factor. The CCNN method performed better than the MCNN method in our experiments, both qualitatively and quantitatively. Although DWI images from both methods were superior to the raw, unprocessed images, in the presence of partial Fourier encoding the MCNN method introduced artifacts that were subtle in DWIs and readily apparent in the parameter maps in Figure 9. This failure could have arisen due to systematic artifacts introduced by the MCNN - an effect that is known to be possible with the use of CNNs [38]. The CCNN model, on the other hand, may be able to use the smoothness of the phase to estimate missing k-space areas - something which is not possible after the magnitude operation has been applied. It may be possible to adapt the MCNN method in the future to handle the ripple effects seen in Figures 8 and 9 - we leave this investigation to future work.

Practically, the success of the approach can be attributed to 1) having a sufficiently deep neural network able to adequately capture the difference between Gibbs-corrupted and ground-truth images and 2) to the heterogeneity of the synthetic training set and the richness of the encoding simulations providing ample training data. During training, the models saw over 10 million candidate phase maps simulated on over 1 million baseline images. Data sets of this size are difficult to acquire in medical imaging, and impossible to acquire for the purpose of simulating Gibbs artifacts in diffusion-weighted imaging. The use of millions of parameters in the design of CNNs, although standard practice in deep learning, includes overfitting risk. Training such a model requires careful consideration of the simulation process and its relationship with the test dataset. As we see in Figures 8 and 9, subtle or benign artifacts in the DWIs are substantially magnified in the corresponding parameter maps. To consider the effects of overfitting, we designed our experiments to have not only separate training, validation, and test data sets as different groups of images, but to have the test phase of the experiments consider images of a completely different class. Prior to our diffusion imaging experiments, we did not use MR images. Implicitly, this suggests that we would observe the same nature of artifacts in healthy subjects as those with pathology. Furthermore, training in this way mitigates the method’s susceptibility to selection bias in the construction of the dataset that would occur with clinical imaging.

Our results were achieved by processing each diffusion-weighted image separately. An advantage of this approach is that our CNN methods can be applied to each image independently for each MR imaging application. Our methods are also readily extendable to other non-diffusion applications that may benefit from Gibbs removal. However, a disadvantage is that the proposed methods do not fully leverage the repetitive correlations present in applications such as diffusion MRI. It is difficult to denoise images at high $b$-values without access to information from low $b$-value images. Such relations are considered in the random matrix denoising method [10]. This suggests that one avenue to explore in the future for diffusion would be to combine the random matrix theory-based denoising [10] as a pre-processing step, able to incorporate the correlations between different DWIs, with the subsequently improved performance of CCNN-based removal of Gibbs and other artifacts.

5 Conclusion

We have shown that Gibbs artifacts and noise can be substantially reduced by training a convolutional neural network with simulations of the image acquisitions. We demonstrated our method on canonical experiments and in vivo data, and showed its potential for artifact removal. The method can be applied independently on each slice of the imaging dataset, enabling its use in many clinical settings such as clinical diffusion MRI where few diffusion directions may be available. In the future, we will examine possible extensions, including consideration of 3D anatomical structures, correlations across repetitive acquisitions, explicit incorporation of physics into the model, and validation on other clinical applications.

Appendix A Phase Simulation

Inspired by results from functional approximation theory [39, 40] and machine learning [41], we approximate the phase as being a summation of Gaussian radial basis functions:

[TABLE]

where $p_{\textit{sim}}(\bm{r})$ is the phase map, $\bm{r}$ is a 2-length vector for spatial position, $a_{b,s}$ is the amplitude of the $(b,s)$th basis function, $\bm{r}_{0,b,s}$ is the center of the $(b,s)$th basis function, and $z_{s}$ governs the width of the basis functions in the $s$th subset. In functional approximation theory, formulations like (4) are used to build a candidate for $p_{\textit{sim}}(\bm{r})$ given a set of points, $(a_{b,s},\bm{r}_{0,b,s}),b\in[1,...,B_{s}]$ ($z_{s}$ is usually chosen to be fixed). However, rather than estimate $p_{\textit{sim}}(\bm{r})$ from a measured cloud of points, our goal is to randomly generate an example from a synthetic cloud of points. For each training example we randomly generate $p_{\textit{sim}}(\bm{r})$ according to the following model:

[TABLE]

where $\mathcal{N}$ indicates a Gaussian distribution, $|\mathcal{N}|$ is a folded Gaussian distribution, and $\mathcal{U}$ indicates a uniform distribution. The parameters $\lambda_{S}$, $\lambda_{B}$, $\mu_{z}$, $\sigma_{z}$, $\mu_{a}$, and $\sigma_{a}$ are set for training, and $r_{\textit{min},j}$ and $r_{\textit{max},j}$ are the minima and maxima for $\bm{r}$ in the $j$th imaging dimension. Contrary to practice in the approximation literature, we allow $z_{s}$ to vary to avoid the degenerate training case where networks create filters to detect Gaussian functions of a fixed width.

When $B_{s}$ and $S$ are not too large, $p_{\textit{sim}}(\bm{r})$ will typically be smooth as in real application settings; however, sharp changes can also occur during training when basis function centers are near each other. As a further form of data augmentation, we include a 1% probability that no phase at all is simulated. Figure 3 shows an example of simulated phase compared to real phase from in vivo data. For our experiments, we used $\lambda_{S}=12$, $\lambda_{B}=15$, $\mu_{z}=64\text{ pixels}^{2}$ $\sigma_{z}=100\text{ pixels}^{2}$, $\mu_{a}=0\degree$, and $\sigma_{a}=5\degree$. We set $r_{\textit{min},j}$ and $r_{\textit{max},j}$ to restrict the center of the Gaussian basis function to the image support in the $j$th dimension.

acknowledgements

We would like to thank E. Owens, M. Drozdzal, M. Tygert, and the rest of the fastMRI team for many useful discussions.

\printendnotes