A partial differential equation-based general framework adapted to Rayleigh's, Rician's and Gaussian's distributed noise for restoration and enhancement of magnetic resonance image.

The proposed framework is obtained by casting the noise removal problem into a variational framework. This framework automatically identifies the various types of noise present in the magnetic resonance image and filters them by choosing an appropriate filter. This filter includes two terms: the first term is a data likelihood term and the second term is a prior function. The first term is obtained by minimizing the negative log likelihood of the corresponding probability density functions: Gaussian or Rayleigh or Rician. Further, due to the ill-posedness of the likelihood term, a prior function is needed. This paper examines three partial differential equation based priors which include total variation based prior, anisotropic diffusion based prior, and a complex diffusion (CD) based prior. A regularization parameter is used to balance the trade-off between data fidelity term and prior. The finite difference scheme is used for discretization of the proposed method. The performance analysis and comparative study of the proposed method with other standard methods is presented for brain web dataset at varying noise levels in terms of peak signal-to-noise ratio, mean square error, structure similarity index map, and correlation parameter. From the simulation results, it is observed that the proposed framework with CD based prior is performing better in comparison to other priors in consideration.

Introduction

The objective of this manuscript is to present the design and development of a general framework for image restoration and enhancement in magnetic resonance imaging (MRI). In MRI Rician noise is one of the prominent noises; however, Gaussian and Rayleigh noise are also present. These type of noises in the MRI can be identified by measuring signal-to-noise ratio (SNR) of image data. In literature, a variety of methods have been described for MRI de-noising, but for the first time, Henkelman[1] presented a method to estimate the noiseless magnitude of MR image from its noisy version data degraded with Rician noise. Estimation of the noise variance from MRI is often of key importance as an input parameter for image post-processing tasks. The estimated noise variance gives a measure of the quality of the MR data. Moreover, the noise variance is often a crucial parameter in image processing algorithms such as noise reduction, segmentation and parameter estimation or clustering.[2]

In literature, for estimation of the noise level in MRI, several methods have been proposed, viz. filtering approach, transform domain approach, and statistical approach. Filtering approach consists of linear filtering and non-linear filtering. While the spatial filter[3] and temporal filter[3] fall under the former category, non-linear filtering includes anisotropic diffusion filter (ADF),[4] adaptive ADF filter,[5] noise driven ADF filter,[6] noise adaptive ADF filter, fourth-order partial differential equation (PDE) filter,[7] adaptive fourth order PDE filter[8] and fourth order complex PDE filters.[9] Nonlocal means (NLM) filter,[10] fast NLM filter,[11] block-wiseoptimized NLM filter,[12] unbiased NLM filter,[13] dynamic NLM filter,[14] enhanced NLM filter,[15] and adaptive NLM filter,[16] the combination of domainand range filters,[17] bilateral domain and range filters,[18] trilateral domain, and range filters[19] are also non-linear filtering.

Examples of transform domain[20] approaches are curvelet,[21] contourlet[22] and wavelet,[23] adaptive multiscale product thresholding,[24] multiwavelet,[25] and undecimated wavelet.[26] Examples of statistical approach are maximum likelihood estimation,[2027] linear minimum mean square error (MSE) estimation,[28] phase error estimation,[29] nonparametric estimation,[30] and singularity function analysis[3132] have been described. Besides these, other noise removal methods proposed in the literature include machine learning-based approaches,[13334353637] discrete cosine transform-based filter,[36] principal component analysis-based technique,[38] and conventional approaches.[39]

The Rayleigh distributed method estimates the noise level in the background. Unfortunately, these methods proved to be useless for images because of unavailability of background information. Except in the brain imaging, background data may not be available in other imaging like in cardiac or in lung, for example in the case that the field of view (FOV)[40] is small, such that noise assumptions based on Rayleigh distribution fail. Most of the noisy background is also eliminated by new scanning techniques and software. These techniques may also affect the methods based on Rayleigh model, which require a certain amount of background pixels for proper estimation of the noise level.[41] Zero-mean Gaussian probability density function (PDF) illustrates the raw complex MR data acquired in the Fourier domain. After the inverse Fourier transform, the noise distribution in the real and imaginary components will still Gaussian due to linearity and the orthogonality of the Fourier transform. However, due to the subsequent transform to a magnitude image, the noise distribution will no longer be Gaussian, but Rician distributed. If I is the original signal amplitude, then the PDF of the reconstructed magnitude image M will be:

where I denotes amplitude of a noise-free image, σ2 denotes the Gaussian noise variance, J0 (.) shows the modified zero order Bessel function, ∈(.) is the unit step Heaviside function and M is the magnitude MRI. The Rician PDF is only valid for nonnegative values of M.[20] In the image background, where the SNR is low (SNR ≈ 0), the Rician PDF reduces to a Rayleigh distribution[42] with PDF:

When SNR is high (>3 dB), then the Rician distribution becomes Gaussian distribution[43] with mean and variance σ2 given as follows:

For the estimation of noise variance, a method based on the local computation of the skewness of the magnitude data distribution was proposed by Rajan et al. in 2010.[42] It is to be concluded that Rician distribution is always in between the moments of Rayleigh and Gaussian distributions. The relationship between σ2 and the variance of a Rician distribution at low and high SNR can be written as:

And σ2 = σ2R      (5)

respectively. In general, σ2 in terms of can be written as:

σ2 = σ2R × Ψ      (6)

where Ψ is a correction factor in the range [1;], i.e., when the Rician distribution approaches a Rayleigh distribution (at low SNR), the correction factor tends to be and when the Rician distribution approaches a Gaussian (at high SNR), the correction factor tends to be 1.

In view of the above discussion and limitations of the existing method, such as noncapability of removal of this type of noise and lower restoration accuracy, in this paper, we proposed a PDE based general framework for restoration and enhancement of MR data.

The proposed method is capable of removing all possible type of noise that may be present in MR data. The manuscript is organized into four sections: Section 1 presents introduction; Section 2 presents methods and model, Section 3 presents the experimental setup, results, and discussions, and Section 4 presents the conclusion of the work.

Methods and Model

The Rician, Rayleigh, and Gaussian noise removal and regularization of MRI data are obtained by minimizing the following nonlinear energy functional of the image I within a continuous domain Ω, using the variational framework:[44]

Where L{p(I/M)} shows the negative likelihood term of Rician or Rayleigh or Gaussian distributed noise in MRI, given by equation (1–2–3). During the filtering process, log-likelihood term measures the dissimilarities at a pixel between M and its estimated value I. L{p(I/M)} acts as the data attachment term or the likelihood term in equation(7).

Maximization of log-likelihood or minimization of the negative log-likelihood leads to de-noising of image data, but is an ill-posed problem and hence regularization is needed. That is why the second term ϕ(∥∇I∥) in equation (7) is needed and it acts as a regularization or penalty function or prior term. In equation (7), λ is a regularization parameter, which has a constant value and makes a balance between the data attachment term and regularization function. The value of λ has been determined experimentally and is set to a value for which peak signal to noise ratio is maximum during the iteration process of filtration. The nonlinear complex diffusion (CD) based, anisotropic diffusion (AD) based and total variation (TV) based prior is a suitable choice for the energy term ϕ(∥∇I∥) based on the concept of the energy function.

ϕ(∥∇I∥) = f(I)      (8)

f(I) is the diffusion PDE based prior obtained by minimization of E (I); from equation (8), substituting the value of ϕ(∥∇I∥) in equation (7) reads:

In case of only Rician noise, after solving modified zero order Bessel function,[45] when we take log and differentiating equation (1) with respect to I, we get the log-likelihood term of Rician's PDF as:

where k1 represents positive integer.

In case of only Gaussian noise, we put value of unit step Heaviside function is one in equation (3), after taking logarithmic of equation (3) becomes:

Differentiating equation (11) with respect to I, we get the log-likelihood term of Gaussian's PDF as:

or

In case of only Rayleigh noise, the log-likelihood term of Rayleigh's PDF proposed by Srivastava and Gupta in 2010[33] is as follows:

Hence, when we combine equations (10), (13) and (14), we get the combined log-likelihood term as follows:

where L'{p (I/M)} shows the negative likelihood term of combined Rayleigh's, Rician's and Gaussian's distributed noise in MRI. When we put the value of likelihood term from equation (15) into equation (9), we get the proposed general framework [Equation (16)] using Euler-Lagrange minimization technique combined with gradient descent approach.

Therefore, the proposed general framework-based model adapted to Rayleigh's, Rician's, and Gaussian's distributed noise reads:

with initial condition

It=0 = I0      (16b)

where λ1, λ2, and λ3 are the constants to be set according to noise pattern, λ is the regularization parameter, and I0 is the noisy image data.

Restoration of magnetic resonance image for different noise distribution

Case 1: Gaussian noise distribution

When λ1 = λ2 = 0 and λ3 = 1 then equation (16) becomes adapted to Gaussian distribution.

with initial condition

It=0 = I0      (17b)

Case 2: Rician noise distribution

When λ1 = λ2 = 1 and λ3 = 0 then equation (16) becomes adapted to Rician distribution.

with initial condition

It=0 = I0      (18b)

Case 3: Rayleigh noise distribution

When λ2 = λ3 = 0 and λ1 = 1 then equation (16) becomes adapted to Rayleigh distribution.

with initial condition

It=0 = I0      (19b)

Figure 1 illustrates the operation of the proposed general framework for restoration and enhancement of MRI data:

Figure 1

Restoration of magnetic resonance image for different noise distribution with different priors

Selection of prior terms

The following three types of diffusion-based prior terms shown in the Figure 2 are used and examined for their efficacy in the proposed methods.

Figure 2

Selection of prior terms

TV based method

AD based method

The CD based method.

Total variation based method

The TV regularization approach was first proposed by Rudin et al.[46] to de-noise an image corrupted with additive white Gaussian noise. In TV based framework for Rayleigh's, Rician's, and Gaussian's noise, the regularization function is defined as:[46]

In discrete case, TV is defined as:

where

and

For numerical implementations, the derivatives can be discretized using standard centered difference approximations and the quantity |∇I|| is replaced with for some small positive value of eps such as 0.00000000001. The value of eps can be assigned to lowest machine number to avoid divide by zero conditions during implementations.

Anisotropic diffusion-based method

In AD-based framework for Rayleigh's, Rician's and Gaussian's noise, the regularization function is defined as:[4]

f(I) = ∇.(c(∥∇I∥)∇I)      (26)

where the diffusion coefficient c(∥∇I∥) is defined as,[29]

where γ is the threshold parameter.

The complex diffusion-based method

In nonlinear CD-based framework for Rayleigh's, Rician's, and Gaussian's noise, the regularization function is defined as:[36]

f(I)=div(c(Im(I))▽I)      (28)

The diffusion coefficient c(Im(I)) is defined as follows:[36]

Here, k is known as threshold parameter and for digital images[36] the value of k ranges from 1 to 1.5. Equation (28) describes the nonlinear CD process, where linear forward diffusion controls the evolution of real part of the images, and both the real and imaginary equations control the evolution of imaginary part of the image.

A qualitative property of edge detection, that is, second smoothed derivative is described by the imaginary part of the image for small value of θ, whereas real values depict the properties of ordinary Gaussian scale-space. For large values of θ, the imaginary part feeds back into the real part creating the wave like ringing effect which is an undesirable property. Here, for experimentation purposes, the value of θ is chosen to be . The adaptive value of edge threshold parameter is used in Equation (29). It is defined as negative exponential distribution:

where α and k0 are constants, usually 1.

Discretization of the proposed model

For digital implementations, the Equations (16a) and (16b) can be discretized using finite differences scheme.[47] For example, the discretized form of TV based proposed model reads:

Similarly, AD-based model can be discretized using finite difference scheme:

Similarly, nonlinear CD model can be discretized using finite difference scheme:

It=0 = I0      (34)

The von Neumann analysis[47] shows that condition requires , for the numerical scheme, given by equation (31–34) to become stable. If the size of the grid is set to be Δx = 1, after that Δt <¼, that is, Δt <0.25. Hence, for the stability of equation (31–34), the value of Δt is set to be 0.24.

Results and Discussion

Brain web database[48] is used for simulated (synthetic) and real (clinical) data sets of normal brain MRIs, to compare the effectiveness of the proposed technique. There are three modalities (pulse sequences) dataset present in the brain web databases[48] which are T1-, T2- and proton density (PD)-weighted. The proposed method and other standard methods used for comparison purposes were implemented using MATLAB R2014.

The performance of restoration results is analyzed for images artificially degraded by mainly Rician's noise and partially Gaussian's noise and Rayleigh's noise if image background is present. Linear minimum mean square error estimator (LMMSE),[49] RLMMSE,[50] and recursive version of signal-to-noise ratio-based nonlocal linear minimum mean square error estimator (RSNLMMSE),[51] are familiar existing techniques used for comparing the proposed method in the case of Rician noise. For Rician noise the best setups as proposed by the authors and the free parameters of these methods are used during experimentation. To obtain the best results the relevant values of the parameters are given below:

LMMSE:[49] window of size 5 × 5, linear minimum MSE estimator

RLMMSE:[50] Recursive version of linear minimum MSE estimator, window of size 5 × 5

RSNLMMSE:[51] window using a 5 × 5, recursive version of SNR-based nonlocal LMMSE.

The parameters are adjusted empirically for de-noising MRIs, and the setup of all the parameters using the proposed scheme is shown in Table 1. The ground truth MR data are artificially contaminated with a noise variance having the range 5%-30% to evaluate the quantitative metrics. Based on structure similarity index map (SSIM) and MSE average restoration results for Rician noise and based on peak SNR (PSNR), MSE, SSIM, and correlation parameter (CP) average restoration results for Gaussian and Rayleigh noise over 4–50 iterations or till the convergence of all these de-noising methods are computed.

Table 1

Parameters setup of the proposed method for de-noising magnetic resonance images

The performance analysis and comparative study of the proposed method with other standard methods are represented on the basis of quantitative results SSIM (MSE) for different levels of Rician noise in Table 2. The value of PSNR, MSE, SSIM and CP represented for different levels of Gaussian and Rayleigh noise in Tables 3 and 4, respectively. In the case of Rician noise SSIM and MSE values show that at low as well as at high rates of Rician noise, the proposed method has much better restoration results than existing methods.

Table 2

Quantitative comparison on simulated magnetic resonance data (brain web) for Rician noise using structure similarity index map (mean square error)

Table 3

Quantitative comparison of proposed method on simulated magnetic resonance data (brain web) for Gaussian noise using peak signal-to-noise ratio, mean square error, structure similarity index map and correlation parameter

Table 4

Quantitative comparison of proposed method on simulated magnetic resonance data (brain web) for Rayleigh noise using peak signal-to-noise ratio, mean square error, structure similarity index map and correlation parameter

Figures 3 and 4 illustrate detailed results, obtained with the close-up view of the restored images for better inspection, approaches, incorporates real image, noisy image, and the restored image. The visual results for simulated MR slice are corrupted with 10% level of Rician noise is presented in Figure 3, Gaussian and Rayleigh noise in Figure 4. On the basis of quantitative and visual results, it is apparent that the proposed approach has produced more accurate results such as more noise removing the ability, and preservation of edges and structural information, at all levels of Rician noise to compare the visual performance, existing and proposed.

Figure 3

Simulated T1-weighted magnetic resonance image with Rician noise. (a) Original image. (b) 10% noisy image. (c) Recursive version of linear minimum mean square error estimator. (d) Recursive version of signal-to-noise ratio-based nonlocal linear minimum mean square error estimator. (e) Linear minimum mean square error estimator. (f) Proposed with total variation. (g) Proposed with anisotropic diffusion. (h) Proposed with complex diffusion

Figure 4

Simulated T1-weighted magnetic resonance image with Gaussian and Rayleigh noise. (a) Gaussian noise corrupted magnetic resonance image. (b) Restored image with proposed method from Gaussian noise corrupted magnetic resonance image. (c) Rayleigh noise corrupted magnetic resonance image. (d) Restored image with proposed method from Rayleigh noise corrupted magnetic resonance image

Retaining the important structural information, such as texture and edges, is considered as an important task in image restoration during noise smoothing process. The detailed information, present in the image do not quantify MSE. A well-known quantitative measure SSIM is used for measure the detail preservation performance of the proposed filter shown in Table 2. The proposed technique is superior in terms of retaining structural information at all noise levels clearly shown in Figure 5a–c.

Figure 5

(a) Structure similarity index map based comparison of T1-weighted modality for Rician noise (b) Structure similarity index map based comparison of T2-weighted modality for Rician noise (c) Structure similarity index map based comparison of proton density-weighted modality for Rician noise

The solution can be computed in one single step (or a few steps for the RLMMSE filter), making it computationally efficient for large data sets; this is the main advantage of the LMMSE[49] filter (and to some extent for the RLMMSE filter). The proposed technique is a deterioration in terms of MSE comparison with the best performing RSNLMMSE[51] and RLMMSE,[50] at a low noise rate (5%), 5.08 for T1, 4.58 for T2, and 5.18 for PD. The result shows that the efficiency of proposed filter also increases as the noise rate increases.

At high noise rates, the proposed technique accurately differentiates the low and high noise regions; hence, the better result obtained. Similarly, in the case of SSIM, Table 2 shows that the proposed scheme outperforms the existing techniques. Comparison of the proposed filter using MSE values are shown in Figure 6a–c, respectively, using simulated data sets. The above figure clearly indicates that the proposed technique is superior at all noise levels.

Figure 6

(a) Mean square error based comparison of T1-weighted modality for Rician noise (b) Mean square error based comparison of T2-weighted modality for Rician noise (c) Mean square error based comparison of proton density-weighted modality for Rician noise

Performance analysis

In this paper, metrics for comparing the performance of various noise reduction schemes in considerations from MRI are defined as follows:[5253]

Mean square error

Where I is the original image without noise, I' is the filtered image, m × n is the size of the image and i = 1....... m, j = 1..........n.

Peak signal-to-noise ratio

Here, RMSE is the root MSE. For optimal performance, measured values of MSE should be small, and that of PSNR should be large.

Structure similarity index map

SSIM is used to compare luminance, contrast, and structure of two different images. It can be treated as a similarity measure of two different images. This similarity measure is a function of luminance, contrast, and structure. The SSIM of two images X and Y be calculated as:

Where µi (i = X or Y) is the mean intensity, σi (i = X or Y) is the standard deviation, and Ci (i = 1 or 2) is the constant to avoid instability when µx2 + µy2 is very close to zero and is defined as Ci = (KiL)2 in which Ki << l and L is the dynamic range of pixel values, for example, L = 255 for 8-bit grayscale image.

Conclusion

A PDE-based general framework filter adapted to Rician noise, Gaussian noise and Rayleigh noise was proposed for restoration and enhancement of MRIs. The proposed filter includes two terms namely data fidelity and prior. The data fidelity term, that is, likelihood term is derived from Rician PDF, Gaussian PDF, and Rayleigh PDF and TV based prior, AD based prior and a nonlinear CD based prior are used. Further, mathematical simplifications have been introduced for likelihood term for efficient implementation of the algorithm. The proposed method was tested on brain web data set for varying noise levels, and performance was evaluated in terms of MSE and SSIM for Rician noise. Similarly, the proposed method also removes Gaussian noise as well as Rayleigh noise. From obtained results and comparative analysis with other standard methods, it is observed that the proposed method is performing better. Further, visual results clearly indicate that the proposed technique has the capability of better noise removal.

Financial support and sponsorship

Nil.

Conflicts of interest

There are no conflicts of interest.