Introduction

In Magnetic Resonance Imaging (MRI), the denoising process is an important step because it impacts the diagnostic accuracy. In many cases the signal-to-noise-ratio (SNR) of the acquired images must be improved. This is carried out by using denoising methods and by preserving most of the image features. In MRI, the images are contaminated with Rician noise with a Probability Density Function (PDF) that depends on the level of the original signal. The problem of noise filtering in MRI has been addressed in different works (^{Pizurica, 2003}; ^{Casaca, 2010}; ^{Pu, 2016}; ^{Weaver, 1991}). In the signal processing literature, several popular denoising algorithms are based on thresholds for coefficients of some transformation (^{Sendur, 2002}). These approaches attempt to preserve significant features of the image while removing most of the noise.

If the wavelet transform is applied on the MR magnitude data directly, both the wavelet and the scaling noisy coefficients are biased estimates of the noise-free coefficients. In ^{Nowak (1999)}, it was suggested that the application of the wavelet transform on the squared MR magnitude image data (which is non-central chi-square distributed) would result in wavelet coefficients being no longer biased estimates of their noise-free counterparts. Although the bias still remains in the scaling coefficients, in many models it is not signal-dependent and the noise can, therefore, be easily removed (^{Pizurica, 2003}; ^{Nowak, 1999}). Similarly, denoising methods based on anisotropic diffusion have been proposed (^{Perona and Malik, 1990}; ^{Samsonov, 2004}) The difficulty with wavelet or anisotropic diffusion algorithms is the risk of over-smoothing the details, particularly in low SNR images (^{Tisdall, 2005}).

Also other strategies have been proposed. In ^{He (2009)}, the estimation of the noise-free coefficients of the MR's magnitude was carried out in a non-local mode using a log-likelihood function. The images are modeled as random fields where pixels with similar neighborhoods come from the same distribution. In, ^{Sijbers (1998)} the authors showed that the Maximum-Likelihood (ML) estimation, which is known to yield optimal results asymptotically, enhances the conventional estimation methods. The proposed method is unbiased for high signal-to-noise ratios (SNRs) and yields physical relevant results for low SNRs.

The concepts of non-local means and dictionary reconstruction have been used in ^{Nair (2014)}; ^{Aarya (2013)}; ^{Elad (2006)}. In ^{Nair (2014)}, a Non-Local Means Maximum Likelihood (NLMML) estimation method was used to estimate the noise-free coefficients from the MR's magnitude images by focusing on preserving edges and tissue boundaries. The method is an improvisation over non-local means maximum likelihood approach for Rician noise reduction in MR images. In ^{Aarya (2013)}, an adaptive filter for Rician noise based on the probability distribution function of the noise and the SNR information of the image was proposed. The filter uses the local statistics of the neighborhood within a mask to perform denoising. Thus, the filter adapts according to the local SNR of the neighborhood. In ^{Elad (2006)}, two training options are considered: using the corrupted image itself, or using a high-quality image database. The K Singular Values Decomposition (KSVD) was used to design an adaptive dictionary trained for natural real images, as well as on patches of the noisy image were also used.

Our method combines two important features, the Non-Local Maximum Likelihood (NLML) estimation and denoising via sparse and redundant representation in an unique process.

The KSVD method is introduced to construct an adaptive dictionary on the noisy image. The patches are grouped using the K-means algorithm and the Principal Components Analysis (PCA) to associate the characteristics (Annexed A). Afterward, the Rician distribution guides each group by using local statistics to perform the denoising process in a Maximum A Posteriory (MAP) formulation. As a restriction over the space of search we estimate the noise coefficients over a residual image resulting from the noise image and an initial filtered image.

The problem of recovering a high-quality image from one or several degraded (e.g, noisy, blurred, and/or down-sampled) versions can be generally formulated by:

Where

y = available distorted image

H = degradation matrix

x = original image vector and

n = noise vector

H can be the identity matrix, a blurring operator or a blurring and down-sampling operator. The studies have been conducted to restore x from y, (^{Yang, 2010}; ^{Sendur, 2002}). For an effective approach, the appropriate prior knowledge of the natural images must used (^{Rudin, 1992}), similarly to the ^{Tikhonov (1963)} and Total Variation (TV) regularization methods (^{Rudin, 1992}). However, these methods tend to over-smooth the restored images due to the rigid rules of a piecewise function. An alternative is the use of sparsity-based regularization (^{Sendur, 2002}). The model assumes that a signal x ∈ R^N can be represented as x ( φα, where φ ∈ R^NxM (N < M) is an over-complete dictionary and most entries of the coding vector α are zero or close to zero. The l₁ -norm based sparse coding problem is generally formulated in the following Lagrangian form:

where the constant λ denotes the regularization parameter. With an appropriate selection of λ, we can obtain a good approximation error and an appropriate sparsity of α. In the scenario of Image Restoration (IR), what we observed is the degraded image signal y. To recover x, first y is sparsely coded with respect to φ for solving the minimization

Therefore, y is reconstructed by y^=ϕαy. The local sparsity constraint ‖α‖1 in Eq. 3 may not lead to an accurate image reconstruction because image sparse coding coefficients α are not randomly distributed. It is due to the internal correlations or redundancies in natural images, (^{Casaca, 2010}; ^{Pu, 2014}). In ^{Pu (2014)}, a group sparse coding scheme to code similar patches simultaneously was proposed achieving impressive denoising results. The KSVD method reported in ^{Aharon (2006)} bring better results than the cosine transform (DCT) in the dictionary construction.

In our work, the Probability Density Function (PDF) of the observed image is linked with an initial noise estimation. The coefficients of the residual noise decomposition are used as a linker between the noisy and the reconstructed image. These two aspects are strong considerations to adapt the denoising method to the Rician taking advantage of an initial estimation, the processing of the residual noise in the domain of the sparse representation and the MAP formulation.

In this paper, the Rician noise distribution is analyzed. The noise reduction is based on the estimation of coefficients of the unknown real values of the signal. For fast processing and accuracy, the process is carried out with a sparse representation of the image. Similar patches of the image are grouped in clusters and treated independently. The main contribution of this work are:

The paper is organized as follows: Describes the patches decomposition of the signal and the sparse representation; exposes the treatment of the noise in the sparse representation and the PDF model of the noisy image; explains the proposed denoising method; discusses the results; and draws the conclusions.

Sparse representation and RICIAN model

In the proposed model, the input image y is first filtered by using a low pass filter with output y_f. Then, the residual image, considering the initial noise approximation, is obtained by yr=y-yf. The patch at i is extracted from an image of the form x_i = R_ix. The R_i matrix extracts the patch x_i of dimension n×n. Given a dictionary φ ∈ R^NxM (N < M). Each patch can be sparsely represented as xi=ϕiαx,i by solving an l₁ minimization problem αx,i=argminαi(‖xi-αi‖22+λ‖αi‖1), the image x can be represented by the set of sparse codes α_x,i with overlapped patches. A redundant patch-based representation of x is obtained by α_x,i as an over-determined system (^{Elad, 2006}). The Eq. 4 yields an overall image reconstructed by averaging each patch.

Noise in the sparse representation domain

The clean and the noisy image can be decomposed as a dictionary and the coefficients as the procedure developed in (^{Aharon, 2006}). In a real situation, only the coefficients of the residual and noisy image can be acquired. The coefficients of the residual image exhibit similar characteristics to the Rician distribution. We propose a Gaussian filter with variance σ ² to obtain the residual image y_r . The dictionaries are learned from the example image patches by using the KSVD (^{Aharon, 2006}). However, it has been shown that sparse coding with an over-complete dictionary is unstable (^{Elad, 2009}). Especially, in the scenario of image restoration. A solution to this is to select k clusters with similar patches and then learn a KSVD sub-dictionary for each one. The similar patches are grouped in consecutive rows of clusters in the matrix P_k. This is carried out on the noisy image, the filtered and the residual images.

In the next step, a β factor in the sparse plane is estimated from the residual image decomposition and then subtracted from the coefficients of the noisy image decomposition.

The residual image y_r is characterized by the previous knowledge of the Rician distribution. With a good estimation of β in y_r, the term α^y≈αx can be used instead of a_x in the Lagrangian model of reconstruction. Then, the coefficients of the clean image can be more approximated by using a minimization of Eq. 2, wherein a negative gradient method is proposed.

Rician noise distribution

If the real and imaginary parts of a signal are corrupted with zero-mean, equal variance and uncorrelated Gaussian noise, the envelope of the magnitude signal will follow a Rician distribution. This kind of noise is apparent in many practical situations. The signal magnitude M is expressed as:

where A is the clean signal level, n₁ and n₂ are uncorrelated Gaussian noise variables with zero mean and equal variance. The PDF of such image is a Rician distribution is expressed as:

with I₀ the 0^th -order modied Bessel function of the first kind and u the Heaviside step function. The moments of the distribution are difficult to calculate, the even-order (non-central) moments are simple polynomials (e.g. second-order moment). The signal with Rician noise is modeled with a probability density function (PDF) that changes the aperture with respect to the real level of signal A.

Proposed method

The connection between the MAP estimator and sparse representation was established in (^{Sendur, 2002}). In our model, the term θ = α_y - β is defined. Using Eq. 5, for a given β, the MAP estimation of θαy≈αx can be formulated as:

According to Eq. 3 and knowing that P(αx,β|αy)P(αy)=P(αy|αx,β)P(αx,β), considering that for little values of β, P(αx,β) =P(αy) and P(αx,β|αy)=P(αy|αx,β) then the prior knowledge with a Gaussian distribution can be modeled as P(θ|αy).

The likelihood term is characterized by the Gaussian distribution, θ and β are assumed to be independent. In the prior probability P(θ|αy), θ reflects the variation of α_y from the estimation of β. If we take β as a good estimation, then θy=αx+β are basically the coefficients α_y of the noisy image. Thus, we can assume that θ follows an independently and identically distributed (i.i.d). process, and the joint prior distribution of the image to be the Rician noise P (θ). The estimated image x^ can be modeled using Eq. 7 for small values of β.

The MAP approximation considers the coefficients of the unknown clean image x in Eq. 10 and the Rician distribution considers α_y for modeling the observed image using the approximation of Eq. 5. Therefore, the proposed MAP formulation is:

The approximation to the MAP model of Eq. 10 is a classical Gaussian function. Then, we can use Eqs. 7 and 8 to approximate to P(αy|αx,β) in the kind of Rician noise and the corrupted image can be modeled as:

If the predicted values of β are considered independent of y, the term P(α_x, β) can be approximated using Eq. 5 in a discrete form as:

Where α_y are the coefficients of the noisy image y, and ( the estimate of the noise. Then, a good estimation of α_yi and β_i, can be computed as the weighted average of the sparse codes associated with the non-local similar patches (including patch i) to patch i. For each patch x_i, we calculate a set of similar patches, denoted by β_i. Then, Eq. 12 is applied for each cluster β_i. The problem is how to obtain a well approximation β_i,j at position j of the group i. Applying logarithms to the MAP formulation in the corresponding localization at the cluster β_i, a good estimation of β maximizes the probability in the function 8 respect to α_x, a discrete formulation of the approximation using the Eq. 9 and 5 is

Determination of β

The parameter β_i is determined on each cluster of α_yr, β can be computed from the sparse codes of the residual output split in patches, the Eq. 14 is used to estimate α_x. The term ρ denotes the sparse codes related to the estimation of β as:

Where G(.) is a low pass filter operator and ρ is the sparse representation of y_r, ρ is characterized by the Rician distribution , then y_r = φρ.

where ω_i,k is a weight. Similar to the non-local means approach, we set the weights to be inversely proportional to the distance between the coefficients ρ_i of the patch y_ri and the coefficients ρ_i,k of the patches y_ri,k in the k^th cluster of patches.

αyri=φρ^i and αyri,k=φρ^i,k are the estimates of the patches ρ_i and ρ_i,k, h is a predetermined scalar and W is the normalization factor. With the non-local redundancies of natural images, we are able to achieve a good estimation of the unknown sparse vectors x_i with a good estimate of β, in the model of Eq. 14.

Stop criterion

The image estimation process is repeated until an approximation to the solution with error ξ is reached. The estimate of Eq. 14 is minimized with respect to β. In this case, we propose a numerical solution as the negative gradient. The negative differential of Eq. 14 is:

−Dαx=α^xl-α^xl-1.

Then a minimum differential is found numerically by using the conjugated gradient.

l is the number of iterations l = 1.ite. When l = 1, β is calculated using the residual y_r and a new assignation is made. The estimated image is assigned: α^y=α^xl.

The algorithm

Results

The proposed method was applied on different MRI studies with relaxation time T2. The images are normalized to 256 gray levels. The Tables 1 and 2 give quantitative results of the comparison of our method respect to other referenced procedures. The quality factors considered to evaluate the results were the ratio of signal peak over noise (PSNR) and the ratio of structural similarity (SSIM). In both cases, the proposed procedure is the most effective filtering the Rician noise with σ =30 in Table 1 and σ =40 in Table 2. Whith σ =30 our method enhance the PSNR in 0.3 dB and SSIM in 0.01 as average respect to the model NCSR. Whith σ =40 our method enhance the PSNR in 0.1 dB and SSIM in 0.1 as average respect to the model NCSR.

Comparisons

The first column of the Figure 1 contain the visual comparison of the result for the slice 3 in the Table 1. The Figure 2 shows the results of detail reproduction of our method with respect to other procedures that produces homogeneous areas. In the Figure 2c the procedure tends to homogenize the areas and eliminate fine detail. It is a model based in threshold over coefficients and dependences between decomposition bands. The Figure 2e shows a model that provides best results in the details reproduction. It employs also clustering techniques but the probabilistic inference is made with a fixed Gaussian function. The results of the Figure 2f lost fine details and the reconstruction tends to expand the structures. It is a model in which the group of coefficients are selected with an adaptive procedure for the estimation. In the case of the Figure 2d the result is a more detailed reconstruction of fine structures in areas from low to median-high levels of signal and in general a reconstruction more near to the original image. The estimation using a Rician noise PDF that change the width proportionally to the level of the signal makes more difficult the estimation with wide PDF. It is the reason for the presence of noise in high levels of signal. This effect could be avoided using equalization, normalizing or changing the ranges of the image’s levels, this technique could be treated in future work.

Figure 1 Visual comparison of results in Table 1, second column left to right 

Figure 2 Visual comparison of results in Table 1, a) Original Image, b) Noisy image, c) Wavelet LMS method, d) Proposed NLRD method, e) NCSR method, f) ASDS method σ = 40 

In the proposed method, the approach of the exhaustive search and clustering is realized over α_yr. This give an increment of the SSIM and PSNR values over all the slices analyzed. The method is robust for high levels of noise, in our comparison it reach σ = 40. The results obtained in the Tables 1 and 2 shows an enhancement of the denoising techniques almost in chronological sense in relation with the publication of the concepts and procedures. The most recent methods use techniques for grouping and reduction of the information to be processed as in ASDS, NCSR and our proposed model Non Local Rician Denoise (NLRD).

The models under comparison use different techniques for noise treatment. The wavelet denoising method uses the background estimation for the level of noise and clusters similar coefficients. Then in the groups, the level of noise is subtracted. In this case, more fine bands are more affected by noise and the upper bands are assumed the desired signal. The PDF in the Rician case combines signal and noise in the model shown in the Eq. 7. Then, the MAP formulation 9 gives the possibility of α^x estimation when a relationship is established between the unknown original coefficients α_x and the noisy coefficients α_y through of (. The MAP estimation is carried out over the sparse domain of clusters of patches of the image. The separation of the image in different groups of patches allows to carry out a precise approximation of ( and the α^x coefficients.

Conclusion

The proposed method combines the PDF of the Rician noise and the redundancies of the structures of the image for denoising. The spatial redundancy of the patches is exploited using the PCA classifier and the KSVD dictionary learning technique. The noise estimation was carried out with a desired precision after some iterations. The Rician noise distribution is analyzed on the image decomposed in sub-regions and clusters. The mean of the group of patches is an important consideration because the Rician distribution changes the width according to the level of the real signal. Due to the width changes of the PDF the model must enhance the response for high levels of signal, despite of this the general result brig a good detail reproduction for the areas of the image. The utility of the sparse representation of the image to characterize the noise distribution was shown. For the estimation of local values, the redundancies around all the residual image are exploited. The method shows that the PDF of the Rician noise can be used in a MAP formulation using an estimation of the residual image as a link between the noisy and the estimated image. The levels of the noisy and estimated images are presents in the PDF. The methods of sparse representation, clustering and noise estimation selected for the procedure are robust for adaptation in the case of Rician noise. It also shows to outperform the results of methods in the state-of-the-art. In future works, the procedure will be tested with other models of noise distributions. Proposing a more precise initial reduction of the noise and the processing with multiresolution decomposition of the image. Other interest of this work is the automatic determination of the hyper-parameters of the PDF by using the input image and the proposed model. Also more work is required to avoid the residual noise in high levels of signal using equalization, normalizing or changing the ranges of the image’s levels.