# Rician Noise Removal via a Learned Dictionary

Authors: Jian Lu; Jiapeng Tian; Lixin Shen; Qingtang Jiang; Xueying Zeng; Yuru Zou
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

This paper proposes a new effective model for denoising images with Rician noise. The sparse representations of images have been shown to be efficient approaches for image processing. Inspired by this, we learn a dictionary from the noisy image and then combine the MAP model with it for Rician noise removal. For solving the proposed model, the primal-dual algorithm is applied and its convergence is studied. The computational results show that the proposed method is promising in restoring images with Rician noise.


---

## Body

## 1. 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–9], Poisson noise [10–12], and multiplicative noise [13–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 y degraded by Rician noise can be given by(1)y=x+η12+η22,where x is the original image and η1,η2~N(0,σ2). Our goal is to find the unknown true image x from the degraded image y 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 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:(2)argminx⁡12σ2∫Ωx2dt-∫Ωlog⁡I0xyσ2dt+λ∫ΩDxdt,where I0 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−l0 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 primal-dual algorithm to solve the model.The following is the outline of our paper. In Section2, 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.


## 2. 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 imagey∈RN×N, regarding the image patch of size n×n, we order it as column vector Y∈Rn lexicographically. And we define a dictionary of size D∈Rn×k to simply construct the Sparseland model, where k>n 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:(3)α^=argminα⁡α0subjecttoDα≈Y.That is, α0≪n, where α0 denotes the number of the nonzero entries in α. And •0 is used to constrain the sparsity of representation.For simplicity, we substituteY-Dα2≤ϵ for Dα≈Y. Also, replacing the constraint with a penalty term, we can get the equivalent problem of (3):(4)α^=argminα⁡Y-Dα22+μα0,and a suitable choice of μ can make problem (3) equivalent to problem (4).Now we consider the entire imagey, that is, to consider all the patches of image y. Then we can construct the sparse representation model for the noisy image y. First, we index the image y with Ω=1,2,…,N2; then the image patches of size n×n located in Ω can be indexed by Γ={1,2,…,N-n+1}2. From patches to image, problem (4) becomes the following problem:(5)α^ij=argminαij⁡∑i,j∈ΓRijy-Dαij22+∑i,j∈Γμijαij0,where Rij is an n×N matrix that we can use to extract the (i, j) patch from the image.Looking back at the assumption that the dictionaryD 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 propose a new model for Rician noise removal which is as follows:(6)argminx∈SΩ⁡β2σ2∫Ωx2dt-β∫Ωlog⁡I0yxσ2dt+12∑Dαij-Rijx2+∑μijαij0+λ∫ΩDxdt,where S(Ω)≔{v∈BV(Ω):0≤v≤255} and ∫ΩDxdt is the total variation (TV) of x. The first and second terms of our model stem from the MAP model, which is data-fidelity 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 x that we need to solve, the dictionary D, and the sparse coefficients αij. 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 imagey, we give the initial dictionary and get the sparse representation coefficients αij; then use αij to learn dictionary D and update corresponding coefficients αij.(2) Use the primal-dual algorithm to get the recovered image that we wanted.


### 2.1. Dictionary Learning

In the first step of our algorithm, we will use the degraded imagey to train a dictionary, and all the image patches can be sparsely represented by the trained dictionary with the corresponding coefficients αij. 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:(7)D^,α^i,j=argminD,αi,j⁡12∑i,jRijx-Dαi,j2+∑i,jμi,jαi,j0.For solving the difficult problem, we have the following specific steps.(1) Initialization. Set x=y, D = overcomplete DCT dictionary.(2) Iteration. F o r   m = 0   t o   N   d o( a ) Given x and D, we get the sparse representation coefficients αij through solving the following problem:(8)α^i,j=argminαi,j⁡12∑i,jRijx-Dαi,j2+∑i,jμi,jαi,j0,and we can efficiently and effectively solve model (8) by using the orthogonal matching pursuit (OMP) method [33, 34].( b ) Given x and αij, we can update the dictionary D=[d1,d2,…,dk] column by column [33, 34]. For each column dl, l=1…k, we update it as follows.( i ) For those patches represented by dl, we write down and denote it by ζl={i,j∣αi,j(l)≠0}.( i i ) For each index (i,j)∈ζl, we compute the corresponding representation error through(9)eijl=Rijx-∑m≠ldmαi,jm,and then we use the columns {eijl}(i,j)∈ζl to define a matrix El.( i i i ) For the matrix El, we apply the singular value decomposition (SVD) into El and get El=UΔVT.Let the first column ofU be d~l to update dl, and multiply the first column of V by Δ(1,1) to update {αi,j(l)}(i,j)∈ζl.End forHere, the dictionary training and sparse representation are completed.


### 2.2. Primal-Dual Algorithm

After the first step of our algorithm, we get the spare dictionary representationDαij from each patch Rijy, and we now minimize (6) with respect to x; that is,(10)argminx∈SΩ⁡β2σ2∫Ωx2dt-β∫Ωlog⁡I0yxσ2dt+12∑Dαij-Rijx2+λ∫ΩDxdt.Proposition 1.
Lety be a bounded function such that infΩy>0; then the objective function in (10) is strictly convex with the constraint σ4/β>supΩy2.Proof.
Using the notationsW=∑(i,j)∈ΓRijTRij and M=∑(i,j)∈ΓRijTDαij, model (10) can be rewritten as(11)argminx∈SΩ⁡β2σ2x,x-βlog⁡I0yxσ2,1Ω+12Wx,x-M,x+λxTV.Also, we define g(x)=β/2σ2x,x-βlog⁡I0yx/σ2,1Ω+1/2Wx,x-M,xx>0 and h(x)=log⁡I0(x)(x>0). According to [32], we have 0<h′′(x)=1-2/xI1(x)/I0(x)-I1x/I0(x)2<1. Then g′′(x)=β/σ2-βlog⁡I0yx/σ2′′+W>β/σ2-βy2/σ4+W. Due to 1≤Wij≤n, g′′(x)>0 holds when σ4/β>supΩy2, and g(x) is strictly convex. We know that x→xTV is convex, so the objective function (11) is strictly convex; that is, function (10) is strictly convex.Referring to Proposition1, we know that model (10) is convex; then we can apply the primal-dual algorithm for solving the minimization problem, which enjoys nice saddle-point structures and has good performance on nonsmooth convex optimization [39–41]. Moreover, we discuss a bias correction technique and the convergence of our algorithm.Briefly, we defineF(x)⁡≔β/2σ2x2-βlog⁡I0yx/σ2,1+1/2∑‍Dαij-Rijx2, and model (10) can be translated to the following format:(12)min0≤x≤255⁡Fx+λ∇x1.Because theTV-model has the duality property, we can change the primal-dual formulation of optimization problem of (12) into the following format:(13)maxp∈P⁡min0≤x≤255⁡Fx-λx,div⁡p,where P=p∈R2d:maxi∈{1,…,d}⁡pi2+pi+d21/2≤1.Here,p is the dual variable, and div=-∇T. The specific algorithm is as in Algorithm 1.Algorithm 1
1. Letκ and τ be given. Set x0=y, x-0=y, and p0=0,…,0T∈R2d. 2. Updatepn, xn, x-n iteratively as follows: (a)pn+1=argmaxp∈P⁡λx-n,div⁡p-12κp-pn22; (b)xn+1=argminx∈S(Ω)⁡F(x)-λx,div⁡pn+1+12τx-xn22; (c)x-n+1=2xn+1-xn. until some stop criterion is satisfied.In particular, the solution of the dual problem 2(a) can be expressed as(14)pin+1=π1λκ∇x¯ni+pini=1,…,2d,where π1(qi)=qi/max⁡1,qi, π1(qd+i)=qd+i/max⁡(1,|qi|), and qi=qi2+qi+d2, i=1,…,d. 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 x∗ equals that of the observed image y. In numerical practice, the following step is implemented in order to preserve the mean of y:(15)xin≔∑j=1dyj∑j=1dmax⁡xjn,0max⁡xin,0i=1,…,d,after updating xn by Newton method, which also ensures that xn≥0.Referring to Theorem 1 in [42], we will get the convergence properties of our algorithm as follows.Proposition 2.
Suppose thatκτλ2∇2<1. Then the iterates {(xn,pn)} of our algorithm converge to a saddle point of (13).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,(16)infx∈SΩ⁡β2σ2∫Ωx2dt-β∫Ωlog⁡I0yxσ2dt+12∑Dαs-RΩsx2+λxTV,where (Ωs)s∈I is a finite set of small patches covering Ω and RΩs is the restriction on Ωs.Theorem 3.
LetΩ be a bounded, open subset of R2 with Lipschitz boundary. Let y be a positive, bounded function. For x satisfying x∈BV(Ω), x>0, let(17)Ex=β2σ2∫Ωx2dt-β∫Ωlog⁡I0yxσ2dt+12∑Dαs-RΩsx2+λxTV,where λ>0. Suppose that (Ωs)s∈I is a finite local coverage of Ω and ∑s∈I∞Dαs2<+∞. Then E(x) has a unique minimizer when σ4/β>supΩ⁡y2.Proof.
LetE(x)=λxTV+Jx and J(x)=β/2σ2∫Ωx2dt-β∫Ωlog⁡I0yx/σ2dt+1/2∑Dαij-Rijx2. By Proposition 1, we know that J(x) is convex. Thus E(x) is convex and bounded from below. Consider a minimizer sequence (xn) for E(x). Then xnTV and J(xn) are bounded from above. Hence (xn) is bounded in BV(Ω). As BV is compact in L1(Ω), there is x0∈BV such that a subsequence (xnk) converges x0 in L1(Ω) and BV(Ω), and we may assume that xnk→x0 pointwise almost everywhere. As the BV norm is l.s.c., x0TV≤liminf⁡xnkTV. Since J(xnk) is bounded from below, by Fatou’s Lemma, we have J(x0)≤liminf⁡Jxnk. Thus Ex0≤liminf⁡Exnk=infBV(Ω)⁡E(x), and x0 minimizes E(x). Consequently, E(x) admits at least one minimizer x0. E(x) is strictly convex function; hence, E(x) admits a unique minimizer.


## 2.1. Dictionary Learning

In the first step of our algorithm, we will use the degraded imagey to train a dictionary, and all the image patches can be sparsely represented by the trained dictionary with the corresponding coefficients αij. 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:(7)D^,α^i,j=argminD,αi,j⁡12∑i,jRijx-Dαi,j2+∑i,jμi,jαi,j0.For solving the difficult problem, we have the following specific steps.(1) Initialization. Set x=y, D = overcomplete DCT dictionary.(2) Iteration. F o r   m = 0   t o   N   d o( a ) Given x and D, we get the sparse representation coefficients αij through solving the following problem:(8)α^i,j=argminαi,j⁡12∑i,jRijx-Dαi,j2+∑i,jμi,jαi,j0,and we can efficiently and effectively solve model (8) by using the orthogonal matching pursuit (OMP) method [33, 34].( b ) Given x and αij, we can update the dictionary D=[d1,d2,…,dk] column by column [33, 34]. For each column dl, l=1…k, we update it as follows.( i ) For those patches represented by dl, we write down and denote it by ζl={i,j∣αi,j(l)≠0}.( i i ) For each index (i,j)∈ζl, we compute the corresponding representation error through(9)eijl=Rijx-∑m≠ldmαi,jm,and then we use the columns {eijl}(i,j)∈ζl to define a matrix El.( i i i ) For the matrix El, we apply the singular value decomposition (SVD) into El and get El=UΔVT.Let the first column ofU be d~l to update dl, and multiply the first column of V by Δ(1,1) to update {αi,j(l)}(i,j)∈ζl.End forHere, the dictionary training and sparse representation are completed.


## 2.2. Primal-Dual Algorithm

After the first step of our algorithm, we get the spare dictionary representationDαij from each patch Rijy, and we now minimize (6) with respect to x; that is,(10)argminx∈SΩ⁡β2σ2∫Ωx2dt-β∫Ωlog⁡I0yxσ2dt+12∑Dαij-Rijx2+λ∫ΩDxdt.Proposition 1.
Lety be a bounded function such that infΩy>0; then the objective function in (10) is strictly convex with the constraint σ4/β>supΩy2.Proof.
Using the notationsW=∑(i,j)∈ΓRijTRij and M=∑(i,j)∈ΓRijTDαij, model (10) can be rewritten as(11)argminx∈SΩ⁡β2σ2x,x-βlog⁡I0yxσ2,1Ω+12Wx,x-M,x+λxTV.Also, we define g(x)=β/2σ2x,x-βlog⁡I0yx/σ2,1Ω+1/2Wx,x-M,xx>0 and h(x)=log⁡I0(x)(x>0). According to [32], we have 0<h′′(x)=1-2/xI1(x)/I0(x)-I1x/I0(x)2<1. Then g′′(x)=β/σ2-βlog⁡I0yx/σ2′′+W>β/σ2-βy2/σ4+W. Due to 1≤Wij≤n, g′′(x)>0 holds when σ4/β>supΩy2, and g(x) is strictly convex. We know that x→xTV is convex, so the objective function (11) is strictly convex; that is, function (10) is strictly convex.Referring to Proposition1, we know that model (10) is convex; then we can apply the primal-dual algorithm for solving the minimization problem, which enjoys nice saddle-point structures and has good performance on nonsmooth convex optimization [39–41]. Moreover, we discuss a bias correction technique and the convergence of our algorithm.Briefly, we defineF(x)⁡≔β/2σ2x2-βlog⁡I0yx/σ2,1+1/2∑‍Dαij-Rijx2, and model (10) can be translated to the following format:(12)min0≤x≤255⁡Fx+λ∇x1.Because theTV-model has the duality property, we can change the primal-dual formulation of optimization problem of (12) into the following format:(13)maxp∈P⁡min0≤x≤255⁡Fx-λx,div⁡p,where P=p∈R2d:maxi∈{1,…,d}⁡pi2+pi+d21/2≤1.Here,p is the dual variable, and div=-∇T. The specific algorithm is as in Algorithm 1.Algorithm 1
1. Letκ and τ be given. Set x0=y, x-0=y, and p0=0,…,0T∈R2d. 2. Updatepn, xn, x-n iteratively as follows: (a)pn+1=argmaxp∈P⁡λx-n,div⁡p-12κp-pn22; (b)xn+1=argminx∈S(Ω)⁡F(x)-λx,div⁡pn+1+12τx-xn22; (c)x-n+1=2xn+1-xn. until some stop criterion is satisfied.In particular, the solution of the dual problem 2(a) can be expressed as(14)pin+1=π1λκ∇x¯ni+pini=1,…,2d,where π1(qi)=qi/max⁡1,qi, π1(qd+i)=qd+i/max⁡(1,|qi|), and qi=qi2+qi+d2, i=1,…,d. 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 x∗ equals that of the observed image y. In numerical practice, the following step is implemented in order to preserve the mean of y:(15)xin≔∑j=1dyj∑j=1dmax⁡xjn,0max⁡xin,0i=1,…,d,after updating xn by Newton method, which also ensures that xn≥0.Referring to Theorem 1 in [42], we will get the convergence properties of our algorithm as follows.Proposition 2.
Suppose thatκτλ2∇2<1. Then the iterates {(xn,pn)} of our algorithm converge to a saddle point of (13).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,(16)infx∈SΩ⁡β2σ2∫Ωx2dt-β∫Ωlog⁡I0yxσ2dt+12∑Dαs-RΩsx2+λxTV,where (Ωs)s∈I is a finite set of small patches covering Ω and RΩs is the restriction on Ωs.Theorem 3.
LetΩ be a bounded, open subset of R2 with Lipschitz boundary. Let y be a positive, bounded function. For x satisfying x∈BV(Ω), x>0, let(17)Ex=β2σ2∫Ωx2dt-β∫Ωlog⁡I0yxσ2dt+12∑Dαs-RΩsx2+λxTV,where λ>0. Suppose that (Ωs)s∈I is a finite local coverage of Ω and ∑s∈I∞Dαs2<+∞. Then E(x) has a unique minimizer when σ4/β>supΩ⁡y2.Proof.
LetE(x)=λxTV+Jx and J(x)=β/2σ2∫Ωx2dt-β∫Ωlog⁡I0yx/σ2dt+1/2∑Dαij-Rijx2. By Proposition 1, we know that J(x) is convex. Thus E(x) is convex and bounded from below. Consider a minimizer sequence (xn) for E(x). Then xnTV and J(xn) are bounded from above. Hence (xn) is bounded in BV(Ω). As BV is compact in L1(Ω), there is x0∈BV such that a subsequence (xnk) converges x0 in L1(Ω) and BV(Ω), and we may assume that xnk→x0 pointwise almost everywhere. As the BV norm is l.s.c., x0TV≤liminf⁡xnkTV. Since J(xnk) is bounded from below, by Fatou’s Lemma, we have J(x0)≤liminf⁡Jxnk. Thus Ex0≤liminf⁡Exnk=infBV(Ω)⁡E(x), and x0 minimizes E(x). Consequently, E(x) admits at least one minimizer x0. E(x) is strictly convex function; hence, E(x) admits a unique minimizer.


## 3. 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:(18)argminx⁡λ∫ΩGσx,y+∫ΩDxdt,(19)Gσx,y=Hσxifx≥cσ,Hσcσ+Hσ′cσx-cσifx≤cσ,Hσ′x=xσ2-yσ2Bxyσ2,Bs≡I1sI0s≈s3+0.950037s2+2.38944ss3+1.48937s2+2.57541s+4.65314,where c=0.8426; I1 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:(20)argminx∈SΩ⁡12σ2∫Ωx2dt-∫Ωlog⁡I0xyσ2dt+1σ∫Ωx-y2dt+λ∫ΩDxdt,where S(Ω)≔{ν∈BV(Ω):0≤ν≤255}. And because BV(Ω) [35, 43] is the subspace of function x∈L1(Ω), the following quantity is finite: (21)Jx=sup⁡∫Ωxtdiv⁡ξtdt∣ξ∈C0∞Ω,R2,ξL∞Ω,R2≤1.NLM model:Before introducing the NLM model, we have a description of the sign that will be used. For a user-defined radiusRsim, we define a square neighborhood window centered around pixel i as Ni. And a Gaussian weighted Euclidian distance of all the pixels of each neighborhood is defined as (22)Lm,n=GρyNm-yNnRsim2,where Gρ 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 Gρ to penalize pixels far from the center of the neighborhood window. And, based on the similarity between the neighborhoods Nm and Nn of pixels m and n, we calculate the similarity f(m,n) as (23)fm,n=1Cme-Lm,n/h2,Cm=∑∀ne-Lm,n/h2,C(m) is the normalizing constant, and h is a exponential decay control parameter.Then, given an imagey, using the NLM method, we can calculate the filtered value at a point m by the following formula: (24)NLMym=∑∀n∈yfm,nyn,0≤fm,n≤1,∑∀n∈yfm,n=1.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.Table 1
Parameterλ values for all testing algorithms.
Images σ = 10 σ = 15 σ = 20 MAP GTV CZ Algorithm1 (β=0.1) MAP GTV CZ Algorithm1 (β=0.1) MAP GTV CZ Algorithm1 (β=1) Lena 15 15 0.07 0.01 20 20 0.055 0.02 25 25 0.045 0.02  Barbara 22 23 0.045 0.01 30 30 4.91 0.02 37 37 0.03 0.02  House 14 14 0.075 0.01 19 19 1.00 0.5 24 24 0.045 10  Monarch 15 15 0.065 0.01 20 20 0.05 0.02 30 30 0.04 0.02  Brain 17 15 0.065 0.01 25 22 0.05 0.01 35 22 0.045 0.01  Mouse 17 17 0.06 0.01 24 24 0.05 0.01 31 32 0.035 0.01In our tests, we choose images “Lena” and “Barbara” with size of512×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(25)PSNRx,x^=20log10⁡2552x-x^2,where x and x^ denote the original image and the recovered image, respectively.Figure 1
The original images.
(a)
Lena (b)
Barbara (c)
House (d)
Monarch (e)
Brain (f)
MouseFirst, 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.Table 2
The comparison of different denoising methods.
Images Methods σ = 10 σ = 15 σ = 20 PSNR SIM Time PSNR SSIM Time PSNR SIM Time Lena Noisy 28.15 0.874  24.65 0.773  22.18 0.680 MAP 34.13 0.957 160.00 32.27 0.935 210.65 30.98 0.915 249.88 GTV 34.10 0.957 13.38 32.21 0.935 13.38 30.90 0.915 13.35 CZ 34.14 0.956 4.11 32.24 0.933 5.04 30.94 0.911 5.24 NLM 34.82 0.961 4.55 32.74 0.939 41.47 31.00 0.914 41.16 Algorithm1 35.39 0.964 480.01 33.36 0.945 433.97 31.59 0.924 443.05  Barbara Noisy 28.16 0.913  24.66 0.839  22.21 0.765 MAP 31.10 0.950 207.92 28.55 0.913 245.68 26.90 0.876 285.53 GTV 31.09 0.950 13.14 28.52 0.913 13.16 26.85 0.875 13.16 CZ 31.07 0.948 4.41 28.49 0.910 4.91 26.83 0.873 5.73 NLM 33.32 0.968 4.79 31.10 0.947 40.90 29.16 0.921 40.00 Algorithm1 34.43 0.972 852.05 32.06 0.956 609.78 29.94 0.932 547.64  House Noisy 28.15 0.605  24.63 0.445  22.14 0.346 MAP 34.42 0.883 34.62 32.65 0.859 42.28 31.36 0.838 57.97 GTV 34.33 0.882 3.43 32.46 0.858 3.35 31.08 0.839 3.35 CZ 34.35 0.881 0.85 32.55 0.858 1.00 31.29 0.836 1.15 NLM 35.08 0.886 10.22 33.23 0.846 10.22 31.39 0.802 10.14 Algorithm1 35.99 0.902 253.91 34.25 0.876 183.07 32.74 0.858 168.34  Monarch Noisy 28.13 0.735  24.67 0.611  22.19 0.517 MAP 33.04 0.934 52.41 30.54 0.907 68.69 29.02 0.865 83.11 GTV 32.99 0.934 3.44 30.44 0.907 3.43 28.87 0.865 3.42 CZ 33.06 0.936 0.99 30.57 0.905 1.20 29.04 0.885 1.39 NLM 32.56 0.940 10.08 30.75 0.908 10.15 29.24 0.871 10.16 Algorithm1 33.87 0.951 803.41 31.35 0.930 487.84 29.60 0.909 349.56  Brain Noisy 27.16 0.657  23.57 0.549  21.13 0.467 MAP 32.95 0.835 21.23 30.08 0.769 21.39 27.97 0.719 23.12 GTV 33.63 0.952 2.57 31.09 0.917 2.14 29.41 0.889 2.17 CZ 33.67 0.950 0.56 31.10 0.919 0.64 29.21 0.891 0.69 NLM 35.15 0.962 6.16 32.07 0.925 6.12 29.72 0.885 5.99 Algorithm1 35.40 0.967 194.60 32.47 0.942 129.26 30.07 0.911 108.05  Mouse Noisy 28.18 0.828  24.67 0.706  22.13 0.592 MAP 31.99 0.930 52.87 29.71 0.890 65.66 28.01 0.847 63.97 GTV 31.94 0.929 3.92 29.60 0.886 3.91 27.89 0.843 4.29 CZ 31.95 0.929 1.17 29.58 0.887 1.31 27.86 0.842 1.52 NLM 32.49 0.935 9.69 29.84 0.888 9.72 27.55 0.833 9.72 Algorithm1 33.71 0.952 700.02 30.63 0.914 397.23 28.12 0.865 295.06The 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 patch-based 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”).Figure 2
Results of “Lena” by different methods. (a) Degraded images (row 1: withσ=10, row 2: with σ=15, row 3: with σ=20), (b) MAP method (row 1: λ=15, row 2: λ=20, row 3: λ=25), (c) GTV method (row 1: λ=15, row 2: λ=20, row 3: λ=25), (d) CZ method (row 1: λ=0.07, row 2: λ=0.055, row 3: λ=0.045), (e) NLM method (row 1-row 3 are denoised image of σ=10,15,20, row 4: patch with σ=20), (f) ours (row 1: λ=0.01, row 2: λ=0.02, row 3: λ=0.02).
(a) (b) (c) (d) (e) (f)Figure 3
Results of “Barbara” by different methods. (a) Degraded images (row 1: withσ=10, row 2: with σ=15, row 3: with σ=20); row 4: patch with σ=20, (b) MAP method (row 1: λ=22, row 2: λ=30, row 3: λ=37); row 4: patch with σ=20, (c) GTV method (row 1: λ=23, row 2: λ=30, row 3: λ=37); row 4: patch with σ=20, (d) CZ method (row 1: λ=0.045, row 2: λ=0.035, row 3: λ=0.03); row 4: patch with σ=20, (e) NLM method (row 1-row 3 are denoised image of σ=10,15,20, row 4: patch with σ=20), (f) ours (row 1: λ=0.01, row 2: λ=0.02, row 3: λ=0.02); row 4: patch with σ=20.
(a) (b) (c) (d) (e) (f)Figure 4
Results of “House” by different methods. (a) Degraded images (row 1: withσ=10, row 2: with σ=15, row 3: with σ=20), (b) MAP method (row 1: λ=14, row 2: λ=19, row 3: λ=24), (c) GTV method (row 1: λ=14, row 2: λ=19, row 3: λ=24), (d) CZ method (row 1: λ=0.075, row 2: λ=0.055, row 3: λ=0.045), (e) NLM method (row 1-row 3 are denoised image of σ=10,15,20, row 4: patch with σ=20), (f) ours (row 1: λ=0.01, row 2: λ=0.5, row 3: λ=10).
(a) (b) (c) (d) (e) (f)Figure 5
Results of “Monarch” by different methods. (a) Degraded images (row 1: withσ=10, row 2: with σ=15, row 3: with σ=20), (b) MAP method (row 1:λ=15, row 2: λ=20, row 3: λ=30), (c) GTV method (row 1:λ=15, row 2: λ=20, row 3: λ=30), (d) CZ method (row 1: λ=0.065, row 2: λ=0.05, row 3: λ=0.04), (e) NLM method (row 1-row 3 are denoised image of σ=10,15,20, row 4: patch with σ=20), (f) ours (row 1: λ=0.01, row 2: λ=0.02, row 3: λ=0.02).
(a) (b) (c) (d) (e) (f)Figures6–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.Figure 6
Results of “Brain” withσ=10 by different methods. (a) is the original “Brain” image and (g) is the degraded image with σ=10. (b)-(f) are the denoised images of MAP method (λ=17), GTV method (λ=15), CZ method (λ=0.065), NLM method, and ours (λ=0.01), and (h)-(l) are residuals of those methods, respectively.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)Figure 7
Results of Brain withσ=15 by different methods. (a) is the original “Brain” image and (g) is the degraded image with σ=15. (b)-(f) are the denoised images of MAP method (λ=25), GTV method (λ=22), CZ method (λ=0.05), NLM method, and ours (λ=0.01), and (h)-(l) are residuals of those methods, respectively.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)Figure 8
Results of Brain withσ=20 by different methods. (a) is the original “Brain” image and (g) is the degraded image with σ=20. (b)-(f) are the denoised images of MAP method (λ=35), GTV method (λ=22), CZ method (λ=0.045), NLM method, and ours (λ=0.01), and (h)-(l) are residuals of those methods, respectively.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)Figure 9
Results of “Mouse” withσ=10 by different methods. (a) is the original “Mouse” image and (g) is the degraded image with σ=10. (b)-(f) are the denoised images of MAP method (λ=17), GTV method (λ=17), CZ method (λ=0.06), NLM method, and ours (λ=0.01), and (h)-(l) are residuals of those methods, respectively.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)Figure 10
Results of “Mouse” withσ=15 by different methods. (a) is the original “Mouse” image and (g) is the degraded image with σ = 15. (b)-(f) are the denoised images of MAP method (λ = 24), GTV method (λ = 24), CZ method (λ = 0.05), NLM method, and ours (λ = 0.01), and (h)-(l) are residuals of those methods, respectively.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)Figure 11
Results of “Mouse” withσ = 20 by different methods. (a) is the original “Mouse” image and (g) is the degraded image with σ = 20. (b)-(f) are the denoised images of MAP method (λ = 31), GTV method (λ = 32), CZ method (λ = 0.035), NLM method, and ours (λ = 0.01), and (h)-(l) are residuals of those methods, respectively.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)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.


## 4. 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 experiments on various images to demonstrate the effectiveness of our model. The numerical experiments show that our proposed model is promising in denoising images with Rician noise.



---

*Source: 8535206-2019-02-18.xml*