Denoising Methods for MRI Imagery

Email: kozaitis@fit.edu Abstract: There are several new denoising algorithms that may be useful for magnetic resonance imaging. We compared the performance of some of the latest denoising algorithms with those that have been developed specifically for MRI imagery. We found that the latest approaches show impressive performance in terms of mean squared error, but a measure based on perceptual quality may be needed to determine the best approach.


Introduction
There have been many recent advances in Magnetic Resonance Imaging (MRI) technology that enable faster and more detailed scans. With this rate of progress, clinicians will soon have the ability to generate even more detailed scans to assist them to make more decisive diagnoses. With these new advances, images of areas such as the lung and cardiac regions will be possible which is encouraging due to their important role in diagnosing disease. MRI imagery unfortunately contains noise that degrades image quality and can hinder diagnosis. Therefore, better methods are needed to reduce or remove noise from these images.
Noise degrades an image by introducing artifacts or reducing the effective resolution and there have been many methods suggested for removing noise in MRI imagery (Mohan et al., 2016). In one method (Yang et al., 2015), noise is transformed into additive noise, then filtered by a conventional Gaussian filter before application of the Non-Local Means (NLM) method (Buades et al., 2005). Next, an inverse transformation is used to reconstruct the final image. The NLM method in this approach compares neighborhoods of voxel values with many other neighborhoods and a filter is created whose coefficients come from a function of the similarity of the neighborhoods and the distance between them. Another approach based on the NLM method uses multiple images from different individuals to collaboratively denoise an MRI image (Geng et al., 2016). This approach reduces the effect of regions that are difficult to compare with the standard NLM filter by increasing the domain so more regions may be available.
In another approach based on the NLM filter (Manjon et al., 2014), the application of Principal Component Analysis (PCA) decomposition over a set of similar patches using a sliding window scheme is performed. The resulting filtered image is used as a guide image to accurately estimate the voxel similarities within a rotationally invariant NLM strategy. This guide image significantly improves the overall denoising performance. Thresholding is then performed by automatically estimating the local noise level from the eigenvalues of the PCA decomposition.
Another method exploits Markov random fields in order to implement a 3D maximum a posteriori estimator of the image (Baselice et al., 2017) that follows the approach in (Martin-Fernandez et al., 2004) developed for diffusion tensor MRI. This approach consists of defining a 3D local Gaussian Markov Random Field that adapts a model consisting of a hyper parameter map to the local behavior of the unknown image by analyzing the 3D neighborhood of each voxel. The model describes the spatial correlation between each pixel and its neighborhood that allows tuning to regularize smooth areas while preserving edges and small details in an unsupervised way. The smoothing effect is automatically tuned by the MRF model in order to find the optimal trade-off between noise reduction and details preservation.
Although the Block-Matching 3D (BM3D) method has been very successful in denoising, it and its variants have not been widely applied to MRI denoising. This method exploits redundancy of patterns in an image (Dabov et al., 2007). In general, an image is divided into small patches and then similar patches are grouped together in a 3D stack. By exploiting the correlation between the image patches within a group, a spare representation can be found and the data effectively filtered, originally with a Wiener filter. A recent MRI reconstruction algorithm that decouples iterations between denoising using the standard BM3D approach without a Wiener filter and reconstruction using an optimization formulation, has been developed (Eksioglu, 2016). Small advances in the standard BM3D method for direct denoising has also been reported (Elahi et al., 2014).
In addition to some MRI denoising methods, we evaluated two recent methods and applied them to MRI imagery. We used a noise model for MRI imagery and compared methods and a variation of the PSNR metric. We used MRI imagery of the brain to compare results. In the next section we describe the noise model and imagery used followed by the denoising algorithms and the results.

Materials and Methods
We considered some of the current leading methods for denoising both MRI and general imagery and compared them against one another.

NLM
The NLM filter uses a spatial domain approach that relies on similar neighborhoods to determine the value of a pixel. The denoised values are based on the mean of the values of all points whose neighborhood looks like the neighborhood of a particular pixel (Buades et al., 2005). The estimated value of a pixel is a weighted average of all pixels in an image that can be described as: where the weights w(x,y,x'y') of the filter depend of the similarity of the neighborhoods around the pixels x, y and x',y' and have previously been described (Buades et al., 2005). The summation is over only those neighborhoods being compared. The variables Z 0 and Z represent the original and denoised images respectively.

BM3D
The BM3D method is a clever way of using selfsimilarity in an image to reduce noise. The first step groups similar blocks of an image such that the entire image is assigned to different blocks. The blocks are arranged to essentially form a three-dimensional (3D) stack, so the image is reduced to many stack of similar blocks. Independently, a 3D transform of each stack is performed that exploits the high degree of correlation in each stack. The resulting sparse transform can then be denoised. The original method used a Wiener filter for the actual noise reduction. After the data has been denoised, patches are reconstructed and they are returned to their original locations. Generally, constant basis functions for all stack are used, independent of the data.

SSC-GSM
One method that has achieved outstanding results modeled sparse coefficients with a Gaussian Scale Mixture (GSM) that allowed the generalization of sparse coefficients to the specification of a sparse prior (Dong et al., 2015). The approach combined this notion with similar patches of an image which were characterized by the same prior. The basic idea behind this GSM-Simultaneous Sparse Coding (GSM-SSC) method is that groups of coefficients from similar patches are exploited for image restoration with the local and nonlocal dependencies.
In this approach basis vectors represent a signal as a linear combination of coefficients with a GSM and a scalar multiplier. The optimization problem reduces to the joint estimation of a Gaussian vectors and scalar multipliers that typically satisfy a l 1 norm constraint that uses the same multiplier for similar image patches. An estimate of the denoised image is used by some method, NLM in this case and variables to be optimized are done so recursively in an efficient manner.

MRI-LMMSE
A method was introduced for MRI data using a Linear Mean Square Error (LMMSE) (Aja-Fernandez et al., 2008) for Rician noise. This approach is realized by noting that the even-order moments of the Rician distribution are simple polynomials as opposed to integral expressions in general. Therefore, using the square of the signal allows a closed-form solution to be realized. One advantage of this approach is that a closedform solution is much faster than iterative optimization methods. In addition, the extension to an arbitrary number of dimensions is straightforward.

Shrinkage Fields
Shrinkage Fields were developed because other highquality methods for denoising based on minimization do not necessarily scale well for large images. Softthresholding functions applied to wavelets have often been referred to as shrinkage functions because they "pull" a function closer to zero. In a current denoising approach (Schmidt and Roth, 2014), shrinkage functions are modeled as a combination of radial basis function kernels. The significance here is that training data is used to determine the model parameters for the shrinkage functions rather than have them assigned manually. By modeling the shrinkage functions, a single quadratic optimization can be used in each iteration of the algorithm. A shrinkage field then can be described whose parameters are determined by the observed image, model parameters and a point spread function. The method is fast because the shrinkage field can be described in closed form rather than relying on more complicated formulations. A cascade of shrinkage fields is referred to in the literature when predictions of shrinkage fields are strung together.

Noise in MRI Imagery
In a single coil system, an MRI image data is generally complex and the magnitude image is used to maximize the SNR. The changing magnetic field from a circularly polarized excitation field results in two signals that differ in phase by π/2 which form the real and imaginary signals. These signals are generally corrupted by Additive Gaussian White Noise (AGWN). After creating the magnitude image, which is simply the square root of the sum of two real and imaginary components, the noise probability distribution becomes Rician (Gudbjartsson and Patz, 1995). The real and imaginary images are described as: where, k is the vector containing the image position variables, N(k) represents the noise with zero mean and standard deviation σ and the magnitude is: The BrainWeb (Kwan et al., 1999) database provides realistic simulated MRI brain data volumes produced by a simulator that can be used to evaluate the performance of various image processing method in a setting where the truth is known. An anatomical model used to simulate MRI data of the brain. It consists of a collection of voxel values for various tissue classes. The values indicate the percentage of tissue classes in each volume. The dimensions are based on 1 mm 3 and each dataset has 181×217×181 pixels. We used a normal anatomical model of the brain with T1-, T2-and proton-density-(PD) weighted imagery using 1 mm slices, 20% level of intensity of non-uniformity and 0, 3 and 9% noise. The percent noise indicates that the standard deviation is a particular percentage of the maximum value of the image. Three simulated MRI images of the same slice for T1-, T2-and PD weighting are shown in Fig. 1-3a. and with 9% noise in Fig. 1

Results
To compare images we initially used the common measure of PSNR which is: where: with the summations in the horizontal and vertical directions over the entire image and M and N are the number of pixels in the horizontal and vertical directions, respectively. The higher the value of PSNR the lower the MSE between the results and original image.
We showed the PSNR results for the images in Table  1 for 3 and 9% noise. For a 3% value of noise, the SSC-GSM performed the best, followed by the BM3D and then the SF method, with each of the methods separated by about 2 dB. The results of the three methods were significantly higher than the results of the NLM and LMMSE methods, with the NLM approach outperforming the LMMSE method.
When 9% noise was used, the SSC-GSM method performed the best followed by the SF method, but the performance of the BM3D dropped significantly. It mostly outperformed the NLM method, which clearly had better results than the LMSSE method.
The expression of MSE in Equation 4 is widely used. It essentially represents the degree of overlap between the input and denoised images. But, if the denoised image is shifted or scaled in some manner compared to the original, the value of MSE will change. To better represent he distribution of energy in the denoised image, we normalized its energy to that of the original before calculating the MSE. Therefore, any energy in the denoised image that is different than the original will contribute to the MSE. When the PSNR was calculated after normalization, we referred to it as the PSNRn.
The results using the PSNRn measure are shown in Table 2. For 3% noise, the overall values were more similar than in the PSNR case. The values for the SSC-GSM and SF methods were lower when the PSNR was used, but the NLM and LMSSE values were higher. The BM3D performed best, closely followed by the SSC-GSM method. The NLM and SF methods performed similarly, about 2dB less and had values larger than the LMSSE method.
For 9% noise, the overall values were also more alike than in the PSNR case. The values for the SS C-GSM and SF methods were lower when the PSNR was used, but the NLM and LMSSE values they were higher as in the case for 3% noise. A notable difference from the 3% noise case is that BM3D case had lower values and had similar performance to the NLM case. The SSC-GSM had the highest values followed by the SF, NLM and BM3D methods and then the LMSSE method.
Although the PSNR and PSNRn provide objective measurements, ultimately visual interpretation of images is what is of interest. In all cases the BM3D results looked somewhat grainy and the NLM results looked the smoothest. The LMSSE had characteristic somewhere between the two despite having relatively low PSNR and PSNRn values. The SSC-GSM and SF results appeared similar, but the SSC-GSM method retained a few more details and was similar to the LMSSE results.  T1  T2  PD  T1  T2  PD  BM3D 36   T1  T2  PD  T1  T2  PD  BM3D 36

Conclusion
The latest methods such as the SSC-GSM and SF approaches performed well both objectively and visually, especially at higher noise values when compared to other approaches. The LMSSE method performed visually very well despite having lower PSNR and PSNR values.
Overall, it appears that the SSC-GSM method performed best by a small amount. A better metric is needed to compare the methods for MRI denoising because the LMSSE approach visually appeared to work well but had low performance values.