Rician Noise Removal via a Learned Dictionary Rician Noise Removal via a Learned Dictionary

,


Introduction
Image denoising is one of the most fundamental issues in image processing [1,2]. In the past decades, many methods are proposed for removing Gaussian white noise [3,4], impulse noise [5][6][7][8][9], Poisson noise [10][11][12], and multiplicative noise [13][14][15][16][17][18][19][20][21]. With the Magnetic Resonance Imaging (MRI) being widely used, people are becoming gradually concerned with another vital noise, Rician noise. In this paper, we mainly study Rician noise and propose a new model that can be more efficient to remove it. Mathematically, the image degraded by Rician noise can be given by where is the original image and 1 , 2 ∼ (0, 2 ). Our goal is to find the unknown true image from the degraded image as well as possible.
In recent years, there have been many methods proposed for denoising on the images corrupted by Rician noise, such as filtering methods, wavelet methods, nonconvex methods, and convex methods. Based on the anisotropic diffusion process proposed by Perona and Malik [22], Gerig [23] presented the nonlinear anisotropic filtering method for Rician noise removal with the edge information preserved, but some other details may be lost. Studying the nonlocal means algorithm, in [24], Prima et al. proposed a nonlocal means variants method for denoising the diffusion-weighted and diffusion tensor MR image. And Manjn et al. [25] also improved a new filter (below we call it NLM model) which is based on the nonlocal means to reduce the noise in MR images. Inspired by it, Wiest-Daessle et al. [26] adapted the nonlocal means filter to data corrupted by Rician noise and presented nonlocal means filtering method for better respecting the image details and structures. In [27], Nowak came up with wavelet domain filtering method which adapts to variations in both the signal and the noise. From a single Rician-distributed image, Foi [28] presented a stable and fast iterative procedure for robustly estimating the noise level and also proposed a variance-stabilization method for efficiently removing Rician noise. In [29], Wood and Johnson proposed wavelet packet denoising method for Rician noise removal at low SNR. Wavelet methods have been efficient to remove Rician noise with the image details and edge preserved, but the problem that small dots influence image analysis process remains unresolved. In the mean time, the maximum 2 Mathematical Problems in Engineering a posteriori (MAP) estimation model was proposed, which is considered from the feature of noise-free image and includes the data fidelity term. The MAP model is as follows: where 0 is the modified Bessel function of the first kind with order zero [30]. But it is a nonconvex function and leads to a difficult problem to solve. In view of the MAP model, GTV model [31] was put forward by Getreuer et al., which is a convex approximation of the MAP model and can be easily solved, but the fidelity item of GTV model is a complicated piecewise function. Chen [32] proposed a new convex model that added a statistical property of Ricain noise into the MAP model, leading to a new strictly convex model under mild condition that can be easily solved by primal-dual algorithm, and below we call it CZ model.
In this paper, we study the Rician noise and propose a new reasonable and efficient model for Rician noise removal. As we know, natural images have a vital feature that is sparseness, and dictionary learning is being widely used for image denoising. Dictionary learning has been demonstrated to be efficient for various noise removal. Aharon and Elad [33,34] proposed the K-SVD algorithm for designing dictionary with sparse representation, and it is proven to be effective for additive white Gaussian noise removal. Inspired by the K-SVD algorithm, Huang [15] proposed a new model that combined the "AA" model [35] and K-SVD algorithm to remove multiplication noise and also presented a log − 0 minimization approach to solve it. Xiao [11] and Ma [10] also proposed new model via dictionary learning for Poisson noise removal and Poisson image deblurring. In addition, Liu et al. [36] applied two-level Bregman method to dictionary updating and proposed an efficient algorithm for reconstructing MR images. In [37], integrating total variation (TV) and dictionary learning, Liu et al. also proposed a novel gradient for image recovery. Similarly, integrating total generalized variation and adaptive dictionary learning, Lu. et al. [38] presented a novel dictionary learning model for MRI reconstruction. So we attempt to apply the sparse representation and dictionary learning into the MAP model for Rician noise removal. Owing to the nonconvexity of the MAP model, we add the sparse representation term to overcome the drawback, so we can use the classical primaldual algorithm to solve the model.
The following is the outline of our paper. In Section 2, we first briefly introduce the dictionary learning and sparse representation, and then we propose the new model that combines the MAP model with sparse representation term. Also, we give and elaborate the two-step algorithm for solving our model. In Section 3, we demonstrate that our model outperforms the other methods for Rician noise removal with numerical results. In the end, we draw our conclusion in Section 4.

Our Proposed Model
Generally speaking, we consider that every signal instance from the family can be represented as a linear combination of few columns from a redundant dictionary. For the degraded image ∈ R √ ×√ , regarding the image patch of size √ × √ , we order it as column vector Y ∈ R lexicographically. And we define a dictionary of size D ∈ R × to simply construct the Sparseland model, where > implies that the dictionary D is redundant. Meanwhile, we should also make an assumption that the dictionary D is known and fixed. Then, column vector Y can be sparsely linearly represented by few atoms selected from the dictionary D.
That is to say that there is a sparse solution of the following problem:̂= That is, ‖ ‖ 0 ≪ , where ‖ ‖ 0 denotes the number of the nonzero entries in . And ‖•‖ 0 is used to constrain the sparsity of representation. For simplicity, we substitute ‖Y − D ‖ 2 ≤ for D ≈ Y. Also, replacing the constraint with a penalty term, we can get the equivalent problem of (3): and a suitable choice of can make problem (3) equivalent to problem (4). Now we consider the entire image , that is, to consider all the patches of image . Then we can construct the sparse representation model for the noisy image . First, we index the image with Ω = {1, 2, . . . , √ } 2 ; then the image patches of size √ × √ located in Ω can be indexed by Γ = {1, 2, . . . , √ −√ +1} 2 . From patches to image, problem (4) becomes the following problem: where R is an × matrix that we can use to extract the (i, j) patch from the image.
Looking back at the assumption that the dictionary D is known and fixed, we have the following question: how to choose and settle the dictionary? In [33], Aharon put forward the K-SVD algorithm for designing dictionary. Giving the initial dictionary, Aharon applied the singular value decomposition (SVD) into updating dictionary. In [34], Elad and Aharon applied the MAP estimator to problem (5) and had a comparison of different dictionaries (overcomplete DCT, global trained dictionary, and adaptive dictionary trained on patches from the noisy image), whose results show that the adaptive dictionary training is best.
After studying the sparse representation and dictionary learning of K-SVD algorithm [33], it inspires us to apply the sparse representation to the MAP model (2), and thus we Mathematical Problems in Engineering 3 propose a new model for Rician noise removal which is as follows: where (Ω) fl {V ∈ (Ω) : 0 ≤ V ≤ 255} and ∫ Ω | | is the total variation (TV) of . The first and second terms of our model stem from the MAP model, which is datafidelity caused by the statistical properties of Rician noise, and the third and forth terms are inspired by the sparse representation. The last TV regularization term can make the denoised image smooth.
In model (6), there are three unknown variables: the noise-free image that we need to solve, the dictionary D, and the sparse coefficients . Similar to [33,34], in order to solve problem (6) effectively and efficiently, here we have the following two-step algorithm.
(1) Based on the degraded image , we give the initial dictionary and get the sparse representation coefficients ; then use to learn dictionary D and update corresponding coefficients .
(2) Use the primal-dual algorithm to get the recovered image that we wanted.
. . Dictionary Learning. In the first step of our algorithm, we will use the degraded image to train a dictionary, and all the image patches can be sparsely represented by the trained dictionary with the corresponding coefficients . The whole process is just using the orthogonal matching pursuit (OMP) and K-SVD algorithms [33,34], and the procedure is to solve the following optimization problem: For solving the difficult problem, we have the following specific steps.

( ) Iteration.
= 0 ( ) Given and D, we get the sparse representation coefficients through solving the following problem: and we can efficiently and effectively solve model (8) by using the orthogonal matching pursuit (OMP) method [33,34].
( ) For those patches represented by , we write down and denote it by = {( , ) | , ( ) ̸ = 0}. ( ) For each index ( , ) ∈ , we compute the corresponding representation error through and then we use the columns { } ( , )∈ to define a matrix E .
( ) For the matrix E , we apply the singular value decomposition (SVD) into E and get E = UΔV T .
Let the first column of U bẽto update , and multiply the first column of V by End for Here, the dictionary training and sparse representation are completed.
. . Primal-Dual Algorithm. After the first step of our algorithm, we get the spare dictionary representation D from each patch R , and we now minimize (6) with respect to ; that is, Proposition 1. Let be a bounded function such that inf Ω > 0; then the objective function in ( ) is strictly convex with the constraint 4 Proof. Using the notations = ∑ ( , )∈Γ R R and = ∑ ( , )∈Γ R D , model (10) can be rewritten as ( ) is strictly convex. We know that → ‖ ‖ is convex, so the objective function (11) is strictly convex; that is, function (10) is strictly convex.
In particular, the solution of the dual problem 2(a) can be expressed as and | | = √ 2 + 2 + , = 1, . . . , . Moreover, the solution of the primal problem 2(b) can be obtained by Newton's method.
Inspired by [39], we provide a bias correction technique for (13) so that the mean of the restored image * equals that of the observed image . In numerical practice, the following step is implemented in order to preserve the mean of : after updating by Newton method, which also ensures that ≥ 0.
Referring to Theorem 1 in [42], we will get the convergence properties of our algorithm as follows.
Inspired by [42], the condition of convergence of our algorithm can be simplified into 2 < 1/8, based on the fact that ‖∇‖ 2 ≤ 8 with the unit spacing size between pixels. For simplicity, in our numerical experiments, we let = 8/ and = 0.015/ and just adjust , and the convergence of the algorithm is satisfied. Before we present numerical results, here we give the basic property of model (10) rewritten in the continuous settings; that is, where (Ω ) ∈ is a finite set of small patches covering Ω and R Ω is the restriction on Ω .

Numerical Experiments
In this section, we will present our numerical experiments to evaluate the approximation accuracy and computational efficiency of our proposed algorithm. We compare our method for the denoising cases with the MAP model (2), the GTV model [31], the CZ model [32], and the NLM model [25]. Here, we have a brief overview of the models that we compared. GTV model: where = 0.8426; 1 is the modified Bessel function of the first kind with order one [30]. Because it is an approximate equality, we can not write down the final explicit restoration model. CZ model: where (Ω) fl {] ∈ (Ω) : 0 ≤ ] ≤ 255}. And because (Ω) [35,43] is the subspace of function ∈ 1 (Ω), the following quantity is finite: NLM model: Before introducing the NLM model, we have a description of the sign that will be used. For a user-defined radius , we define a square neighborhood window centered around pixel as . And a Gaussian weighted Euclidian distance of all the pixels of each neighborhood is defined as where is a normalized Gaussian weighting function with zero mean and standard deviation (usually set to 1). By giving more weight to pixels near the center, we can use to penalize pixels far from the center of the neighborhood window. And, based on the similarity between the neighborhoods and of pixels and , we calculate the similarity ( , ) as ( ) is the normalizing constant, and ℎ is a exponential decay control parameter.
Then, given an image , using the NLM method, we can calculate the filtered value at a point by the following formula: The parameter values of MAP model, GTV model, CZ model, and our proposed model are listed in Table 1. In order to preserve the mean of the observed image, the bias correction technique (15) is utilized for the CZ method and our method. All the experiments are performed under Windows 7 and MATLAB R2017a running on a PC equipped with 2.90GHz CPU and 4G RAM.
In our tests, we choose images "Lena" and "Barbara" with size of 512 × 512, "House" and "Monarch" with size of 256 × 256, "Brain" with size of 181 × 217, and "Mouse intestine" (briefly, we call it "Mouse" throughout this paper) with size of 248 × 254, which are shown in Figure 1. We evaluate the quality of recovered images obtained from various denoising algorithms by using the structural similarity index (SSIM) [44] and the peal-signal-to-noise ratio (PSNR) [45] defined by PSNR ( ,̂) = 20 log 10 ( where and̂denote the original image and the recovered image, respectively. First, we choose the test image of "Lena" and degrade it by Rician noise with = 10 and use different methods to recover the noisy image. The numerical results are enumerated in Table 2, where the first column gives the test image names, the second column gives the method names (here, Algorithm 1 denotes our method presented in Section 2), and columns 3-5 (resp., 6-8 and 9-11) are the PSNR (dB), SSIM, and CPU time (s) for = 10 (resp., = 15, = 20). We observe that the PSNR values of the restored images by Algorithm 1 are higher than those of MAP, GTV, CZ, and NLM method. Moreover, you can easily find that the PSNR values of the restored images by Algorithm 1 are always the highest among all those five methods for all images. As far as "Barbara" is concerned, the PSNR value of our method is more than 3 dB higher than CZ model at = 10. The characteristics of SSIM values are almost consistent with those of PSNR values. However, the SSIM values are composed of luminance comparison, contrast comparison, and structure comparison, which makes differences at some points.
The images are corrupted, respectively, by Rician distribution with = 10, = 15, and = 20, and the restored images of the above five algorithms for the images "Lena", "Barbara", "House", and "Monarch" are shown in Figures 2-5, respectively. We clearly find that much noise remains in the    results of MAP, GTV, and CZ method. For example, in Row 4 of Figure 3, we can see from the background that "Barbara" image restored by our method is more clearer than those restored by other methods, and textures of Barbara's trousers and scarf are kept better in the restored image by our method. The restored images by our method also preserve more significant details than MAP, GTV, CZ, and NLM methods on the hat tidbits of the "Lena" image. What is more, the flowers of the lower left of "Monarch" also indicate that our method is superior to other methods, because our method is a patchbased method, which is a different framework with TV-based methods (i.e., MAP method and CZ method). The images obtained from our method provide smoother regions and better shape preservation (e.g., the backgrounds in "Lena" and "Barbara"). Figures 6-11 are results of experiments on "Brain" and "Mouse" images. In particular, for the "Brain" and "Mouse" images, we not only present the recovered images but also use residuals to make a comparison with the original images. The first row is original image and recovered images of MAP method, GTV method, CZ method, NLM method, and our method. The second row is degraded image with various and residuals of those methods, respectively. In Figures 6-8, the recovered images of of MAP method, GTV method, CZ method, and NLM method are still unclear, so we can conclude that the recovered images by our method are best through our visions and residuals. As for the "Mouse" image, the recovered results are very similar; in this case, we use residuals to see which result is better. In general, the blurrier residual is, the better according recovered image is. We can find that the outline of residual by our method is blurriest.
In conclusion, the results of our numerical experiments demonstrate that our proposed method performs better than the MAP method, GTV method, CZ method, and NLM method.

Conclusion
In this paper, we proposed a new effective model via a learned dictionary for denoising images with Rician noise. More specifically, based on the dictionary training, we add a dictionary penalty term to the original nonconvex MAP model to establish a new denoising model and develop a two-step algorithm to solve our model. Also we carry out

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there are no conflicts of interest. Laboratory of Intelligent Information Processing, Shenzhen University, China (518060).