Multiplicative Noise Removal via a Learned Dictionary
4534
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 21, NO. 11, NOVEMBER 2012
Multiplicative Noise Removal via a Learned Dictionary Yu-Mei Huang, Lionel Moisan, Michael K. Ng, and Tieyong Zeng
Abstract— Multiplicative noise removal is a challenging image processing problem, and most existing methods are based on the maximum a posteriori formulation and the logarithmic transformation of multiplicative denoising problems into additive denoising problems. Sparse representations of images have shown to be efficient approaches for image recovery. Following this idea, in this paper, we propose to learn a dictionary from the logarithmic transformed image, and then to use it in a variational model built for noise removal. Extensive experimental results suggest that in terms of visual quality, peak signal-to-noise ratio, and mean absolute deviation error, the proposed algorithm outperforms state-of-the-art methods. Index Terms— Denoising, dictionary, multiplicative noise, sparse representation, variational model.
I. I NTRODUCTION
I
MAGE denoising is one of the most fundamental issues in image processing. In the past decades, a great deal of methods emerged for removing additive white Gaussian noise, for example, total variation-based methods [1]–[3], dictionary-based methods [4], [5], wavelet-based methods [6], [7], nonlocal-based methods [8], etc. In recent years, multiplicative noise removal has attracted much attention, see, for instance [9]–[17]. Mathematically, the image f degraded by multiplicative noise is the multiplication of the clean image u with the noise η, f = uη
(1)
where the images u and f , noise η are functions from to R+ with ⊂ R2 . Multiplicative noise arise in many imaging applications, such as synthetic aperture radar (SAR), ultrasound imaging [18]–[20], electronic microscopy [21] and Manuscript received November 11, 2011; revised April 29, 2012; accepted May 18, 2012. Date of publication June 18, 2012; date of current version October 12, 2012. This work was supported in part by the NSFC Grant 11101195 and Grant 11171371, the Specialized Research Fund for the Doctoral Program of Higher Education of China under Grant 20090211120011, and the China Post-Doctoral Science Foundation funded Project 2011M501488. The work of M. K. Ng was supported in part by the RGC grants and HKBU FRGs. The work of T. Zeng was supported by the RGC Grant 211710, Grant 211911, and Grant HKBU FRGs. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Jose M. Bioucas-Dias. Y.-M. Huang is with the School of Mathematics and Statistics, Lanzhou University, Lanzhou 730000, China (e-mail: [email protected]). L. Moisan is with the Laboratoire MAP5 (CNRS UMR 8145), Université Paris Descartes, Paris 75006, France (e-mail: [email protected]). M. K. Ng and T. Zeng are with the Department of Mathematics, Hong Kong Baptist University, Kowloon, Hong Kong (e-mail: [email protected]; [email protected]). Digital Object Identifier 10.1109/TIP.2012.2205007
laser images. In this paper, we focus on the multiplicative noise η that follows the Gamma distribution (see [9], [11]) with the probability density function K K η K −1 −K η e (2) (K ) where K is an integer related to the noise level and is the classical Gamma function. In the literature, several methods have been proposed for multiplicative noise removal. Most of these methods can be classified into two categories: the methods derived from the classical MAP Bayesian rule (see [9], [11]), and the efficient methods that convert the multiplicative noise removal problem into an additive noise removal problem (see [22]). The following is a brief review of these two classes of methods. pdf(η) =
A. MAP-Based Multiplicative Noise Removal The image denoising problem consists in estimating the original image u from the observed image f . A classical statistical approach based on the maximization of the posterior distribution Prob(u| f ) leads to the MAP estimate uˆ = arg max Prob(u| f ). u
Equivalently, uˆ is a minimizer of − log Prob(u| f ), and using Bayes rule eventually leads to the variational problem: f γ inf logu + + φ(u) u>0 u Z where we assume that the image u follows the Gibbs prior γ (u) = Z1 e−γ φ(u) with the normalization factor Z . Inspired by the above MAP analysis, several variational methods have been proposed for Gamma multiplicative noise removal problems. A rather popular one is the Aujol-Aubert (“AA” ) method [9], which estimates the original image u through the following minimization problem: f + λuT V . inf log u + (3) u>0 u Here, λ is the positive regularization parameter and uT V is the Total Variation (TV) regularization term. In a continuous setting, and for a smooth image u, it reads |∇u(ξ )|dξ uT V :=
where |∇u| = (∂u/∂ξ1 )2 + (∂u/∂ξ2 )2 . The main drawback of the objective function in (3) is that it is not a convex function of u. The optimal solution
1057–7149/$31.00 © 2012 IEEE
HUANG et al.: MULTIPLICATIVE NOISE REMOVAL VIA A LEARNED DICTIONARY
4535
of (3) may not be unique. Moreover, classical optimization algorithms provide only a local minimum of the objective function, which may not necessarily be the global optimal solution and depends on initial guesses of iterative optimization methods. Inspired by the “AA” model, Huang et al. [11] reformulated the problem in the logarithm domain (log u) and obtained a convex optimization problem. Then they designed an alternating minimization method to solve this optimization problem. Numerical results showed that this convex model can provide denoised images with higher SNRs and better visual quality, independently of initial guesses. Recently, in [10], a split-Bregman algorithm was presented to solve (3) directly and it is numerically quite efficient. In [23], a simple and fast method using graph-cut minimization was studied.
The main goal of this paper is to propose an effective method for multiplicative noise removal which combines the MAP formulation, the sparse representation principle, and the idea of converting the multiplicative noise removal into the additive noise removal. The outline of the paper is as follows. In Section II, we first review the K-SVD method which is currently one of the best performing ones for additive noise removal. In Section III, a new model is proposed for multiplicative noise removal, and some basic properties of it are discussed in Section IV. In Section V, we show experimental results to demonstrate the quality of the denoised images and the efficiency of the proposed method. Finally, concluding remarks are given in Section VI.
B. Methods Based on Additive Noise Removal
In this section, let us briefly review the K-SVD√method √ in a discrete setting, that is, when is a grid of size N × N . This method was originally proposed by Aharon and Elad [4] for additive Gaussian √ √ noise removal. Suppose that a√noisy N × N indexed by = {1, 2, . . . , N }2 , image f of R written as a column vector f ∈ R N , results from a zeromean Gaussian noise b ∈ R N with a known standard deviation σ superimposed on an original image u ∈ R N . The basic assumption of √ the K-SVD method is that each image patch √ (of fixed size n × n, and rewritten as a column vector in Rn ) can be represented sparsely as a linear combination of atoms taken from a fixed dictionary D ∈ Rn×k . Here dictionary means a set of column vectors in Rn . Each column vector is called an atom and is usually of unit norm. The total number of these atoms is k. A typical dictionary may be an overcomplete DCT dictionary (sampling cosine waves in different frequencies to obtain a fixed number of atoms, see [4], [26] for details), or may be learned from a noisy image (see below). Let us assume that √ √ (6) P = {1, 2, . . . , N − n + 1}2 √ √ is the set of n × n image patches located in . Given a vector x = (x 1 , x 2 , . . . , x m ) ∈ Rm , we shall consider its classical l 2 norm x, and also its l 0 “norm”
II. A DDITIVE K-SVD D ENOISING A LGORITHM Since there have been a number of good methods to remove additive noise in the literature, it is natural to consider converting the multiplicative noise denoising problem into the additive noise denoising problem. The main idea is to take the logarithm of both sides of equation (1), leading to log f = logu + logη.
(4)
This is an additive noise removal problem, for which the probability distribution of the noise logη is given (see, e.g., [9]) by K K e K (logη−e pdf(logη) = (K )
logη )
.
(5)
In [14], Shi and Osher extended the inverse scale-space method to multiplicative noise removal models by converting multiplicative noise into additive noise (4), and their proposed model is strictly convex. Furthermore, the convergence of the flow for the multiplicative noise model was demonstrated, as well as its regularization effect and its relation to the Bregman distance. Durand et al. [22] proposed an efficient hybrid approach based on the log-image data and a reasonable underoptimal hard-thresholding on its curvelet transform. They minimized a criterion composed of an l 1 data-fidelity term (on the thresholded frame coefficients) and a TV regularization term in the log-image domain. The main interest of their method is that it combines the advantages of shrinkage and variational methods. Theoretical results on the hybrid criterion were presented, and a Douglas-Rachford splitting algorithm was used to minimize the energy function. The numerical results demonstrated the superiority to the other compared methods, especially for images containing tricky geometrical structures. Moreover, Steidl and Teuber [15] considered using the so called I-divergence as a data fitting term, and TV or nonlocal means as a regularization term in the cost function for removing multiplicative Gamma noise. A nonlocal filter for removing multiplicative noise was proposed by Teuber et al. [16]. Furthermore, Teuber et al. presented a new similarity measure for nonlocal filters in [24] and the denoising results revealed some good properties of the filters. Some other nonlocal-based filters are further elaborated in [25].
x0 := #{i |1 ≤ i ≤ m, x i = 0}, that counts the number of nonzero entries in x 1 , . . . x m . Under the sparsity assumption, the Gaussian noise removal task can be described as the following energy minimization task: ˆ u} {αˆ i j , D, ˆ = arg min λ f − u2 {D,αi j ,u}
+
1 Dαi j − Ri j u||2 + μi j αi j 0 . 2 (i, j )∈P
(i, j )∈P
(7) Rn×N
In this model, Ri j ∈√ √is a binary n × N matrix which extracts from u a n × n patch located in (i, j ) ∈ P (thus, Ri j u ∈ Rn ). The first term of (7) demands a proximity between the recorded image f and its denoised unknown version u. The second term demands that each patch Ri j u of the reconstructed image u can be represented by the dictionary
4536
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 21, NO. 11, NOVEMBER 2012
D ∈ Rn×k by coefficients vector αi j ∈ Rk with a bounded error. The third part is a (weighted) sparsity constraint, which demands that the number of coefficients required to represent any patch is small. Note that the values μi j are patch-specific weights that are determined hiddenly by the optimization procedure. The minimization (7) with respect to the unknowns D, αi j , u yields the denoising algorithm. The choice of D is of high importance for the performance of the algorithm. In [4], it is shown that dictionary training can be done by the minimization problem (7). The algorithm proposed there is an iterative block-coordinate relaxation method which fixes all the unknowns except for the one to be updated. Thus, (7) is solved with a three-step iterative method: learning typical patterns (dictionary) from the noisy image via K-SVD; representing patches sparsely via the Orthogonal Matching Pursuit (OMP) algorithm [27]; averaging the denoised patches to finally obtain a global denoised image. III. P ROPOSED M ODEL Inspired by K-SVD [4] for dictionary learning and the “AA” model [9], here we propose a log-l 0 minimization approach, that is, an approach based on the sparse representation (l 0 ) of the log-image over a certain unknown dictionary of image patches. Denoting by√ 1 the √ constant 1 over the discrete image domain (a N × N grid) and by ·, · the usual scalar product in Euclidean spaces, the proposed model reads f ˆ αˆ i j , u} { D, ˆ = arg min λ log u + , 1 + γ loguT V u {D,αi j ,u>0} 1 + Dαi j − Ri j log u||2 2 (i, j )∈P + μi j αi j 0 (8) (i, j )∈P
where λ, γ are positive regularization parameters, P is given in (6) and u ∈ R N is the estimated image. The · T V term is defined, in the discrete setting, by summing over the image domain the norm of ∇u, the classical 2-neighbors discrete gradient estimate (see, e.g., [1]). As in (7), Ri, j ∈ Rn×N is the matrix corresponding to the extraction of the patch located in (i, j ), and αi, j ∈ Rk is the sparse vector of coefficients to represent the patch Ri, j log u with the dictionary D ∈ Rn×k . The hidden parameters (μi j )(i, j )∈P are determined by the optimization procedure described in [4]. The interpretation of the new model (8) is rather clear. The first to a constant, is the discrete version of the term term, up f (log u+ ) u in model (3); the second term is inspired by [11], [14], [15] where the non-convex issue of (3) is avoided; the last two terms assume that the small patches of the transformed image log u can be represented sparsely over a certain adaptive dictionary D. As TV is efficient to recover cartoon images while dictionaries are well adapted to textures, we can expect that Equation (8), which defines a kind of “cartoon + texture” model, will lead to good results. Similarly to [4], we propose a two-phases algorithm to approximate this log-l 0 minimization problem:
1) Learn the dictionary D from the log image and represent the log image patches via the learned dictionary to get Dαi, j . 2) Use a fast iterative method to solve a reduced convex model to get the recovered image. A. Learning a Dictionary from the Degraded Image and Obtaining its Sparse Representation In the first step of our method, a dictionary D which contains typical patterns of the observed image log f is obtained through the K-SVD method, and each patch (i, j ) of log f is represented sparsely with the coefficient vector αi, j in this dictionary [4], [11]. This is just the OMP and K-SVD steps in the additive K-SVD denoising algorithm, and the procedure is as follows. For the set of image patches Ri, j log f , we seek a dictionary D ∈ Rn×k and αi, j that minimize the cost function 1 {D, αi, j } = arg min Ri, j log f − Dαi, j 2 {D,αi, j } 2 (i, j )∈P + μi, j αi, j 0 , (9) (i, j )∈P
where μi, j is the hidden parameter which indicates how those two forces should be weighted. An iterative algorithm is designed to solve (9) effectively [4]. Here we note that in the proposed model, the sparse representation is performed in the logarithm domain of the image. Since the log function is strictly monotone, f and log f have the same edges, so we can expect that the efficiency of the K-SVD method remains on log f [11]. B. Dictionary and Total Variation Based Multiplicative Noise Removal Using the sparse dictionary representation Dαi, j obtained for each patch Ri, j log f via K-SVD in the first phase, we now minimize (8) with respect to u, that is, f arg min λ log u + , 1
u u>0 1 + Ri, j log u − Dαi, j 2 + γ loguT V . 2 (i, j )∈P
(10) As we shall see later, the experimental results show that the proposed model (10) is effective. However, the criterion is not convex with respect to u, and thus the global minimizer is difficult to obtain. Inspired by [11], [14], let us introduce the transformation u = exp(y). Then we can obtain a new equivalent model, which is now convex with respect to the variable y: arg min λ y + f e−y , 1
y
+
1 Ri, j y − Dαi, j 2 + γ yT V . (11) 2 (i, j )∈P
Using the notations RiTj Ri j , W = (i, j )∈P
M=
(i, j )∈P
RiTj Dαi, j
HUANG et al.: MULTIPLICATIVE NOISE REMOVAL VIA A LEARNED DICTIONARY
(note that W and M have the same dimensions as u and f ), Equation (11) can be rewritten as arg min λ y + f e y
−y
, 1
1 + W y, y − M, y + γ yT V . (12) 2 Let us denote by i, j the coefficient (i, j ) of any matrix . We can notice that since for any (i, j ) ∈ , the weight Wi, j counts how many times pixel (i, j ) is used to construct √ the √ image patches of size n × n, we have 1 ≤ Wi, j ≤ n. Proposition 1: Assume that inf f > 0, then the objective function in (12) is strictly convex. Proof: The function t ∈ R → t + f i, j e−t is convex, hence so is y → y + f e−y , 1 ; y → − M, y is also convex; y → W y, y is strictly convex since Wi, j ≥ 1; and classically y → yT V is convex. As a sum of convex and strictly convex functions, the objective function (12) is strictly convex. Thanks to the Proposition above, the minimization problem (12) can be handled by many typical optimization algorithms. Below we will design an efficient dual approach to solve it. Denote X = R N . The gradient ∇u is a vector in the vector space Y = X × X. And for p = ( p1 , p2 ), q = (q 1 , q 2 ) ∈ Y , we define a scalar product by pi,1 j qi,1 j + pi,2 j qi,2 j . p, q Y = √ 1≤i, j ≤ N
We also consider the discrete divergence operator div : Y → X, which is the opposite of the adjoint of the discrete gradient operator [1]. Indeed, we have − div = ∇ ∗ , that is, ∀u ∈ X, ∀ p ∈ Y, ∇u, p Y = − u, div p X . Notice that classically, we have yT V = ∇ y1 = max ∇ y, p
p∞ ≤1
where similarly to [28], we define p∞ , p1 as | pi, j | p∞ = max | pi, j |, p1 = i, j
4537
Algorithm 1 Chambolle-Pock Algorithm Solving (12) 1: Initialization: y¯ (0) , y (0), p(0) . 2: For k = 0, 1, 2, . . ., repeat until a stopping criterion is reached p(k+1) := argmin p∞ ≤1 {− p, ∇ y¯ (k) + β p − p(k) 2 }, y (k+1) := argmin y∈X {G(y) − div p(k+1) , y
+βy − y (k) 2 }, y¯ (k+1) := 2y (k+1) − y (k) , 3: Output: yˆapprox := y (k+1) .
formula, while the second step can be handled with Newton’s method efficiently. Theoretically, if the positive parameter β is chosen big enough, then the above algorithm converges to the solution of (13). √ Proposition 2: If β > 2, then Algorithm 1 converges to the unique minimizer of (12). Proof: The problem (12) can be rewritten as min G(y) + ∇ y1 . y∈X
This is exactly the nonlinear primal problem studied in [28] [see (3) there] and one can carefully verify that the Program 1 here is the Algorithm 1 proposed there. By Theorem 1 of [28], the algorithm converges as soon as 1 2 1 · L 2 since L ≤ 2 2 (see [1], [28]). The functional proposed in (10) is a dictionary + TV-based multiplicative noise removal model. The dictionary approach permits to restore textures efficiently, while TV is known to preserve edges well. As we shall see, the proposed model combines both merits together and manages to efficiently restore images corrupted by a multiplicative Gamma noise.
i, j
IV. BASIC P ROPERTIES
with | pi, j | = ( pi,1 j )2 + ( pi,2 j )2 . Now let us consider 1 1 G(y) = λ y + f e−y , 1 + W y, y − M, y . γ 2 We can easily prove that the dual problem of (12) is max min G(y) − div p, y .
p∞ ≤1 y∈X
(13)
Inspired by [28], [29], an iterative splitting method can be applied to solve (13). Starting from an initialization y¯ (0) , y (0), p(0) the sequence of iterates y (1), p(1) , y (2) , p(2) , · · · can be obtained from Program 1. Note that in this algorithm, the first step can be updated by means of a closed-form
Before we present numerical results, let us study the basic properties of the model (10) rewritten in the continuous settings, that is, f log u + λ inf log u∈B V (),u>0 u 1 + Rs log u − Dαs 2 2 s∈I
+ γ loguT V
(14)
where (s )s∈I is a finite set of small patches covering , and Rs is the restriction on s . We need a simple lemma. Lemma 3: Let a > 0 and M0 = min{0, 2 + 2 log( a2 )}. Then ∀x ∈ R, x + ae−x ≥ |x| + M0 .
4538
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 21, NO. 11, NOVEMBER 2012
Proof: Set g(x) = x + ae−x . Then g(x) − x ≥ 0 for any x ∈ R, and g(x) + x = 2x + ae−x attains its minimal value 2 + 2 log a2 in x = log a2 . We thus have ∀x ∈ R, g(x) − |x| = min(g(x) − x, g(x) + x) ≥ M0 as announced. Proposition 4: Let be a bounded open subset of R2 with Lipschitz boundary, and f a bounded function such that inf f > 0. For u satisfying log u ∈ BV (), u > 0, let f E(u) := log uT V + λ log u + . u Then E(u) has a unique minimizer. Proof: Let v = log u, then E(u) = F(v) with F(v) := vT V + λJ (v) and J (v) := (v + f e−v ). By Lemma 3, we have (15) J (v) ≥ (v + (inf f )e−v ) ≥ v1 + M0 ||
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
for some constant M0 (depending on f ). Thus, F is convex and bounded from below. Consider a minimization sequence (v n ) for F. Then v n T V and J (v n ) are bounded from above. By (15), so is v n 1 . Hence, (v n ) is bounded in BV (). As BV is compact in L 1 (), there is a v ∈ BV such that a subsequence (v nk ) converges to v in L 1 , and we may assume (up to another extraction) that v nk → v pointwise almost everywhere. As the BV norm is l.s.c., vT V ≤ lim inf v nk T V . Since v nk + f e−v nk is bounded from below, by Fatou’s Lemma, we have J (v) ≤ lim inf J (v nk ). Thus F(v) ≤ lim inf F(v nk ) = inf B V () F, and v minimizes F. Consequently, E admits at least one minimizer, u = ev . Readily, F is a strictly convex function. Hence F (and consequently E) admits a unique minimizer. Theorem 5: Let be a bounded, open subset of R2 with Lipschitz boundary. Let f be a positive, bounded function. For u > 0 satisfying log u ∈ BV (), let E(u) : = γ log uT V f 1 +λ Rs u log u + + u 2 s∈I
−Dαs
2
(16)
where γ > 0. Suppose that (s )s∈I is a finite local coverage 2 of and s∈I Dαs < +∞. Then E(u) has a unique minimizer. Proof: Consider a minimization sequence (u n ) of E. Let v n = log u n . Similarly to the previous proposition, one can take a subsequence (v nk ) which converges both in L 1 and BV . V. N UMERICAL R ESULTS In this section, the denoising results are presented to demonstrate the efficiency of the proposed method. Comparisons are also made with the most efficient methods that have been proposed recently for multiplicative noise removal. These five methods are: the “BS” algorithm [30]; the “AA” method [9]; the “DFN” algorithm [22]; the “DRSM” method [15] and the updating nonlocal method “NL” [16]. The parameters of all
Fig. 1. Comparative results for the Cameraman image (noise level: K = 10). (a) Original. (b) Noisy (K = 10). (c) AA (24.48). (d) BS (24.43). (e) DFN (26.08). (f) DRSM (25.93). (g) NL (26.55). (h) Proposed (27.32).
(a)
(b)
(c)
(d)
(e)
(f)
Fig. 2. Sample of the recovered camera stand part of Cameraman corresponding to Fig. 1 when K = 10. (a) Original. (b) BS (24.43). (c) DFN (26.08). (d) DRSM (25.93). (e) NL (26.55). (f) Proposed (27.32).
the methods are selected with respect to the best recovered PSNR and MAE values (see definitions below). Except for the experiments appearing in [22], the parameters used in the “AA,” “BS” and “DFN” methods are the same as in this paper. As in [22], two quantities are to assess the quality √ used √ of the recovered images. For a N -by- N pixels noise-free image u and its denoised version uˆ by a given algorithm, we define the peak signal-to-noise ratio (PSNR, in decibels) by PSNR = 10 log10
N · |max(u) − min(u)|2 (in dB) uˆ − u2
where |max(u)−min(u)| measures the gray-scale range of the original image. The Mean Absolute-deviation Error (MAE)
HUANG et al.: MULTIPLICATIVE NOISE REMOVAL VIA A LEARNED DICTIONARY
(a)
(b)
(d)
4539
(a)
(b)
(c)
(d)
(e)
(f)
(c)
(e)
(f)
Fig. 5. Sample of the recovered yard part of Nîmes corresponding to Fig. 4 when K = 10. (a) Original. (b) BS (27.0). (c) DFN (27.80). (d) DRSM (27.95). (e) NL (28.42). (f) Proposed (28.85).
Fig. 3. Sample of the recovered camera stand part of Cameraman when K = 4. (a) Original. (b) BS (22.36). (c) DFN (22.98). (d) DRSM (23.67). (e) NL (24.20). (f) Proposed (24.91).
(a)
(a)
(b)
(c)
(d)
(e)
(f)
(b) Fig. 6. Sample of the recovered yard part of Nîmes when K = 4. (a) Original. (b) BS (24.92). (c) DFN (25.84). (d) DRSM (25.75). (e) NL (26.07). (f) Proposed (26.34).
(c)
(d)
(e)
(a)
(f)
(g)
(b)
(h)
Fig. 4. Comparative results for the Nîmes image (noise level: K = 10). (a) Original. (b) Noisy (K = 10). (c) AA (25.26). (d) BS (27.0). (e) DFN (27.80). (f) DRSM (27.95). (g) NL (28.42). (h) Proposed (28.85).
(c)
(d)
(e)
(f)
(g)
(h)
is defined by 1 uˆ − u1 . N The images Cameraman with size 256-by-256, Nîmes (a French city) with size 512-by-512, Peppers with size 256by-256 and Barbara (head part) with size 256-by-256 are used in our experiments. The images are contaminated by multiplicative Gamma noise with two different levels K = 10 and K = 4 respectively. The original images are shown in (a) of Fig. 1, 4, 7, 10, respectively. The corresponding corrupted images are shown in (b) of Fig. 1, 4, 7, 10, respectively, and the corresponding recovered images by the “AA,” “BS,” “DFN,” “DRSM,” “NL” and the proposed method are shown MAE =
Fig. 7. Comparative results for the Peppers image (noise level: K = 10). (a) Original. (b) Noisy (K = 10). (c) AA (24.40). (d) BS (25.09). (e) DFN (26.40). (f) DRSM (26.52). (g) NL (27.07). (h) Proposed (27.30).
in subparts (c) to (h) of these figures respectively. Certain corresponding zoomed samples of the recovered images by the “BS,” “DFN,” “DRSM,” “NL” and the proposed method
4540
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 21, NO. 11, NOVEMBER 2012
(a)
(b)
(d)
(c)
(e)
(a)
(b)
(c)
(d)
(e)
(f)
(f)
Fig. 8. Sample of the recovered base part of Peppers corresponding to Fig. 7 when K = 10. (a) Original. (b) BS (25.09). (c) DFN (26.40). (d) DRSM (26.52). (e) NL (27.07). (f) Proposed (27.30).
(a)
(b)
(c)
(d)
(e)
(f)
Fig. 9. Sample of the recovered base part of Peppers when K = 4. (a) Original. (b) BS (22.82). (c) DFN (23.81). (d) DRSM (24.19). (e) NL (24.32). (f) Proposed (24.78).
(a)
Fig. 11. Sample of the recovered shoulder part of Barbara corresponding to Fig. 10 when K = 10. (a) Original. (b) BS (22.55). (c) DFN (23.80). (d) DRSM (23.06). (e) NL (24.60). (f) Proposed (25.47).
(a)
(b)
(c)
(d)
(e)
(f)
(b)
Fig. 12. Sample of the recovered shoulder part of Barbara when K = 4. (a) Original. (b) BS (20.97). (c) DFN (21.60). (d) DRSM (21.63). (e) NL (22.70). (f) Proposed (23.03).
(c)
(d)
(e)
(f)
(g)
(h)
Fig. 10. Comparative results for the Barbara image (noise level: K = 10). (a) Original. (b) Noisy (K = 10). (c) AA (22.56). (d) BS (22.55). (e) DFN (23.80). (f) DRSM (23.06). (g) NL (24.60). (h) Proposed (25.47).
are shown in subparts (b) to (f) of Fig. 2, 3, 5, 6, 8, 9, 11, 12, respectively. The recovered PSNR values are provided in the subpart captions of each figure.
Since the readability of the whole recovered images is rather poor, we only show them in the case K = 10 and essentially base our observations on close-ups. Moreover, due to the inferior performance of the “AA” method, we do not show the corresponding images for this method in the close-up figures. Tables I to IV report two performance assessments (the PSNR and MAE measurements, with the best results in boldface) for each of the six considered methods, along with the computation time (in seconds) reported in the “TIME” column. For the proposed method, the computing time is expressed as the sum of the dictionary training time and the remaining computation time. Both in terms of visual effect and PSNR/MAE measurements, it clearly appears that the proposed method outperforms the five others.
HUANG et al.: MULTIPLICATIVE NOISE REMOVAL VIA A LEARNED DICTIONARY
TABLE I P ERFORMANCE E VALUATION FOR THE Cameraman I MAGE , T IME FOR
P ROPOSED M ETHOD = D ICTIONARY T RAINING + R EMAINING
4541
TABLE III P ERFORMANCE E VALUATION FOR THE Peppers I MAGE , T IME FOR P ROPOSED M ETHOD = D ICTIONARY T RAINING + R EMAINING
Recovered K
10
4
Method
Recovered
PSNR
MAE
TIME(s)
AA BS DFN
24.48 24.43 26.08
9.90 9.44 7.40
19.35 1.00 58.10
DRSM NL Proposed
25.93 26.55 27.32
7.24 6.82 6.63
AA BS DFN
21.97 22.36 22.98
DRSM NL Proposed
23.67 24.20 24.91
K
Method
PSNR
MAE
TIME (s)
AA BS DFN
24.40 25.09 26.40
10.91 9.90 8.49
21.25 0.99 54.17
5.62 4.67 16.93 + 2.90
DRSM NL Proposed
26.52 27.07 27.30
8.13 7.62 7.57
5.96 5.74 13.60 + 2.93
14.26 12.14 10.61
19.58 0.99 53.74
AA BS DFN
22.68 22.82 23.81
13.67 13.08 11.52
20.41 1.00 52.67
9.23 8.97 8.97
5.65 4.57 11.55 + 2.94
DRSM NL Proposed
24.19 24.32 24.78
10.47 10.15 10.14
6.11 5.72 10.66 + 1.61
10
4
TABLE II P ERFORMANCE E VALUATION FOR THE Nîmes I MAGE , T IME FOR
TABLE IV P ERFORMANCE E VALUATION FOR THE Barbara I MAGE , T IME FOR
P ROPOSED M ETHOD = D ICTIONARY T RAINING + R EMAINING
P ROPOSED M ETHOD = D ICTIONARY T RAINING + R EMAINING
Recovered K
10
4
Method
Recovered
PSNR
MAE
TIME (s)
AA BS DFN
25.26 27.0 27.80
8.83 7.69 7.21
113.0 3.9 243.2
DRSM NL Proposed
27.95 28.42 28.85
6.96 6.59 6.33
AA BS DFN
24.55 24.92 25.84
DRSM NL Proposed
25.75 26.07 26.34
K
Method
PSNR
MAE
TIME (s)
AA BS DFN
22.56 22.55 23.80
12.82 12.38 10.86
15.45 1.06 56.30
29.46 19.25 30.12 + 9.60
DRSM NL Proposed
23.06 24.60 25.47
11.65 9.82 8.99
6.79 3.95 16.39 + 1.24
10.06 9.87 9.09
120.4 3.7 240.3
AA BS DFN
18.69 20.97 21.60
20.05 15.16 14.21
15.31 0.92 55.04
8.93 8.75 8.60
29.49 18.07 18.31 + 9.43
DRSM NL Proposed
21.63 22.70 23.03
13.87 12.34 12.04
6.90 2.37 13.17 + 1.71
For the Cameraman image with K = 10 (see Fig. 1, Fig. 2 and Table I), the estimated PSNR for “DRSM” and “NL” methods are 1.39 and 0.77 smaller than with the proposed method, while the MAE are 0.61 and 0.19 larger. In the case K = 4 (again, for Cameraman image), the best PSNR and MAE scores are obtained with the proposed method too. Comparing the visual appearance of the recovered images in Fig. 2 and Fig. 3, we find that the “NL” method yields a good result too, but the camera stand is better recovered and the artifacts are smaller with the proposed method. We also observe that the “AA” method over-smooths the image and sometimes leaves outliers in the recovered images; “BS” and “DFN” give rise to more artifacts in the recovered images; the recovered effect obtained by “DRSM” is not as good as the results of the “NL” and the proposed method. For the larger size picture Nîmes, Fig. 5 and Fig. 6 clearly show that the proposed method can preserve some details in the yard part of the images, while the same part is more obscure in the images recovered with the other methods.
10
4
For the picture Peppers, the recovered results shown in Fig. 7, Fig. 8 and Fig. 9 display the same recovered qualities as for Cameraman and Nîmes images: less artifacts appear in the recovered images and the textures are best preserved with the proposed method. Concerning Barbara image (zoom on head part, as in [15]), Fig. 10, Fig. 11 and Fig. 12 further show the superiority of the proposed method to recover the texture of the original image. Fig. 13 show the recovered effect when TV is not used in the proposed method [γ is zero in (8)] for the restoration of Barbara image in the case K = 4. Clearly, the result obtained with TV has less artifacts than without TV, and this is especially clear for the flat region. Therefore, the combination of TV regularization and dictionary learning improves the efficiency of the proposed method to a certain degree. All the recovered results show that the proposed method can effectively remove the multiplicative noise in the degraded images, and it manages to preserve well both edges and textures at the same time.
4542
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 21, NO. 11, NOVEMBER 2012
R EFERENCES
(a)
(b)
(c)
Fig. 13. Sample of the recovered hand part of Barbara when K = 4, with and without the TV term in (8). (a) Original. (b) Without TV (22.65). (c) Proposed (23.03).
Fig. 14.
Dictionary learned for the (log) image Cameraman when K = 10.
Concerning the results obtained with the method proposed in [11], as expected they are not sensitive to the initialization of the iterative algortihm used, but we observed that they were just a little better than those obtained by the “AA” method. Also, the method proposed in [14] always takes a long time to complete the restorations. Both methods are not as good as the methods “DFN,” “DRSM,” “NL” and the proposed method. Note that √ in all experiments, the parameter λ was set to be λ = 1/ (1, K ), where is the logarithmic derivative of the Gamma function, β was set to 10, and γ = 1. The proposed method is much less sensitive to the parameter values. To illustrate the mechanism of the proposed method, we show the dictionary learned from the logarithm of the Cameraman image (K = 10) in Fig. 14. VI. C ONCLUSION In this paper, we considered the multiplicative noise removal problem. Writing a MAP formulation for the log-transformed image with a combination of a TV-regularization and a sparse representation in an adaptive dictionary of image patches, we obtained a new variational formulation that combines the edge-preserving efficiency of TV-regularization with the ability of dictionary-based methods to recover textures well. Using the Orthogonal Matching Pursuit and an adaptation of Chambolle-Pock algorithm, we built an algorithm that permits to efficiently remove multiplicative noise and outperforms the most recent ones, both visually and in quantitative terms. ACKNOWLEDGMENT The authors would like to thank S. Durand, J. Fadili, M. Nikolova, and Dr. Teuber, for kindly providing the datum and original codes of their papers. The authors would also like to thank the anonymous reviewers for their extremely helpful comments.
[1] A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imag. Vision, vol. 20, nos. 1–2, pp. 89–97, 2004. [2] T. Chan and J. Shen, Image Processing and Analysis: Variational, PDE, Wavelet, and Stochastic Methods. Philadelphia, PA: SIAM, 2005. [3] L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Phys. D, vol. 60, nos. 1–4, pp. 259–268, 1992. [4] M. Elad and M. Aharon, “Image denoising via sparse and redundant representations over learned dictionaries,” IEEE Trans. Image Process., vol. 15, no. 12, pp. 3736–3745, Dec. 2006. [5] S. Lintner and F. Malgouyres, “Solving a variational image restoration model which involves L ∞ constraints,” Inverse Problems, vol. 20, no. 3, pp. 815–831, 2004. [6] D. Donoho and I. Johnstone, “Ideal spatial adaptation via wavelet shrinkage,” Biometrika, vol. 81, no. 3, pp. 425–455, 1994. [7] J. Starck, E. Candes, and D. Donoho, “The curvelet transform for image denoising,” IEEE Trans. Image Process., vol. 11, no. 6, pp. 670–684, Jun. 2002. [8] A. Buades, B. Coll, and J. Morel, “A review of image denoising algorithms, with a new one,” J. Multi. Model. Simul., vol. 4, no. 2, pp. 490–530, 2005. [9] G. Aubert and J. Aujol, “A variational approach to remove multiplicative noise,” SIAM J. Appl. Math., vol. 68, no. 4, pp. 925–946, 2008. [10] J. Bioucas-Dias and M. Figueiredo, “Multiplicative noise removal using variable splitting and constrained optimization,” IEEE Trans. Image Process., vol. 19, no. 7, pp. 1720–1730, Jul. 2010. [11] Y. Huang, M. Ng, and Y. Wen, “A new total variation method for multiplicative noise removal,” SIAM J. Imag. Sci., vol. 2, no. 1, pp. 20–40, 2009. [12] K. Krissian, C. Westin, R. Kikinis, and K. Vosburgh, “Oriented speckle reducing anisotropic diffusion,” IEEE Trans. Image Process., vol. 16, no. 5, pp. 1412–1424, May 2007. [13] L. Rudin, P. Lions, and S. Osher, “Multiplicative denoising and deblurring: Theory and algorithms,” in Geometric Level Sets in Imaging, Vision, and Graphics. New York: Springer-Verlag, 2003, pp. 103–119. [14] J. Shi and S. Osher, “A nonlinear inverse scale method for a convex multiplicative noise model,” SIAM J. Imag. Sci., vol. 1, no. 3, pp. 294– 321, 2008. [15] G. Steidl and T. Teuber, “Removing multiplicative noise by DouglasRachford splitting methods,” J. Math. Imag. Vision, vol. 36, no. 2, pp. 168–184, 2010. [16] T. Teuber and A. Lang, “Nonlocal filters for removing multiplicative noise,” in Scale Space and Variational Methods. New York: SpringerVerlag, 2012, pp. 50–61. [17] Y. Yu and S. T. Acton, “Speckle reducing anisotropic diffusion,” IEEE Trans. Image Process., vol. 11, no. 11, pp. 1260–1270, Nov. 2002. [18] V. Dutt and J. Greenleaf, “Adaptive speckle reduction filter for logcompressed B-scan images,” IEEE Trans. Med. Imag., vol. 15, no. 3, pp. 802–813, Dec. 1996. [19] M. Isana, R. Wagner, B. Garra, D. Brown, and T. Shawker, “Analysis of ultrasound image texture via generalized rician statistics,” Opt. Eng., vol. 25, no. 6, pp. 743–748, 1986. [20] R. Wagner, S. Smith, and J. Sandrik, “Statistics of speckle in ultrasound B-scans,” IEEE Trans. Sonics Ultrason., vol. 30, no. 3, pp. 156–163, May 1983. [21] A. Kryvanos, J. Hesser, and G. Steidl, “Nonlinear image restoration methods for marker extraction in 3d fluorescent microscopy,” Proc. SPIE, vol. 5674, pp. 432–443, Mar. 2005. [22] S. Durand, J. Fadili, and M. Nikolova, “Multiplicative noise removal using L 1 fidelity on frame coefficients,” J. Math. Imaging Vision, vol. 36, no. 3, pp. 201–226, 2010. [23] L. Denis, F. Tupin, J. Darbon, and M. Sigelle, “SAR image regularization with fast approximate discrete minimization,” IEEE Trans. Image Process., vol. 18, no. 7, pp. 1588–1600, Jul. 2009. [24] T. Teuber and A. Lang, “A new similarity measure for nonlocal filtering in the presence of multiplicative noise,” to be published. [25] C. A. Deledalle, L. Denis, and F. Tupin, “Iterative weighted maximum likelihood denoising with probabilistic patch-based weights,” IEEE Trans. Image Process., vol. 18, no. 12, pp. 2661–2672, Dec. 2009. [26] F. Malgouyres and T. Zeng, “A predual proximal point algorithm solving a non negative basis pursuit denoising model,” Int. J. Comput. Vision, vol. 83, no. 3, pp. 294–311, 2009. [27] Y. Pati, R. Rezaiifar, and P. Krishnaprasad, “Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition,” in Proc. 27th Asilomar Conf. Signals Syst. Comput., vol. 1. 1993, pp. 629–639.
HUANG et al.: MULTIPLICATIVE NOISE REMOVAL VIA A LEARNED DICTIONARY
[28] A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging,” J. Math. Imag. Vision, vol. 40, no. 1, pp. 120–145, 2011. [29] M. Zhu and T. Chan, “An efficient primal-dual hybrid gradient algorithm for total variation image restoration,” Dept. Math., UCLA, Los Angeles, CAM Rep. 08-34, 2008. [30] C. Chesneau, M. Fadili, and J. Starck, “Stein block thresholding for image denoising,” Appl. Comput. Harmon. Anal., vol. 28, no. 1, pp. 67–88, 2010.
Yu-Mei Huang received the M.Sc. degree in computer science from Lanzhou University, Lanzhou, China, in 2000, and the Ph.D. degree in computational mathematics from the Hong Kong Baptist University, Kowloon, Hong Kong. She is currently an Associate Professor with Lanzhou University. Her current research interests include signal and image processing, scientific computing, and numerical linear algebra.
Lionel Moisan received the Ph.D. degree from Paris-Dauphine University, Paris, France, and the HDR degree from Paris XI University, in 1997 and 2003, respectively. He studied mathematics and computer science from the École Normale Supérieure, Paris. He was a CNRS Researcher with École Normale Supérieure de Cachan, Bruz, France, from 1998 to 2003, and is currently a Professor of applied mathematics with Paris-Descartes University. His current research interests include variational models for image processing (restoration, smoothing) and a-contrario statistical models for image analysis (gestalt and structure detection).
4543
Michael K. Ng received the B.Sc. and M.Phil. degrees from the University of Hong Kong, Hong Kong, in 1990 and 1992, respectively, and the Ph.D. degree from the Chinese University of Hong Kong, Hong Kong, in 1995. He was a Research Fellow with the Computer Sciences Laboratory, Australian National University from 1995 to 1997, and an Assistant and Associate Professor with the University of Hong Kong from 1997 to 2005. He has been a Professor with the Department of Mathematics, Hong Kong Baptist University since 2005. His current research interests include bioinformatics, image processing, scientific computing, and data mining. Dr. Ng serves on the editorial boards of international journals.
Tieyong Zeng received the B.S. degree from Peking University, Beijing, China, the M.S. degree from École Polytechnique, Palaiseau, France, and the Ph.D. degree from the University of Paris XIII, Paris, France, in 2000, 2004, and 2007, respectively. He was a Post-Doctoral Researcher with the Centre de Mathématiques et de Leurs Applications, ENS de Cachan, France. He is currently an Assistant Professor with Hong Kong Baptist University, Hong Kong. His current research interests include image processing, statistical learning, and scientific computing.
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 21, NO. 11, NOVEMBER 2012
Multiplicative Noise Removal via a Learned Dictionary Yu-Mei Huang, Lionel Moisan, Michael K. Ng, and Tieyong Zeng
Abstract— Multiplicative noise removal is a challenging image processing problem, and most existing methods are based on the maximum a posteriori formulation and the logarithmic transformation of multiplicative denoising problems into additive denoising problems. Sparse representations of images have shown to be efficient approaches for image recovery. Following this idea, in this paper, we propose to learn a dictionary from the logarithmic transformed image, and then to use it in a variational model built for noise removal. Extensive experimental results suggest that in terms of visual quality, peak signal-to-noise ratio, and mean absolute deviation error, the proposed algorithm outperforms state-of-the-art methods. Index Terms— Denoising, dictionary, multiplicative noise, sparse representation, variational model.
I. I NTRODUCTION
I
MAGE denoising is one of the most fundamental issues in image processing. In the past decades, a great deal of methods emerged for removing additive white Gaussian noise, for example, total variation-based methods [1]–[3], dictionary-based methods [4], [5], wavelet-based methods [6], [7], nonlocal-based methods [8], etc. In recent years, multiplicative noise removal has attracted much attention, see, for instance [9]–[17]. Mathematically, the image f degraded by multiplicative noise is the multiplication of the clean image u with the noise η, f = uη
(1)
where the images u and f , noise η are functions from to R+ with ⊂ R2 . Multiplicative noise arise in many imaging applications, such as synthetic aperture radar (SAR), ultrasound imaging [18]–[20], electronic microscopy [21] and Manuscript received November 11, 2011; revised April 29, 2012; accepted May 18, 2012. Date of publication June 18, 2012; date of current version October 12, 2012. This work was supported in part by the NSFC Grant 11101195 and Grant 11171371, the Specialized Research Fund for the Doctoral Program of Higher Education of China under Grant 20090211120011, and the China Post-Doctoral Science Foundation funded Project 2011M501488. The work of M. K. Ng was supported in part by the RGC grants and HKBU FRGs. The work of T. Zeng was supported by the RGC Grant 211710, Grant 211911, and Grant HKBU FRGs. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Jose M. Bioucas-Dias. Y.-M. Huang is with the School of Mathematics and Statistics, Lanzhou University, Lanzhou 730000, China (e-mail: [email protected]). L. Moisan is with the Laboratoire MAP5 (CNRS UMR 8145), Université Paris Descartes, Paris 75006, France (e-mail: [email protected]). M. K. Ng and T. Zeng are with the Department of Mathematics, Hong Kong Baptist University, Kowloon, Hong Kong (e-mail: [email protected]; [email protected]). Digital Object Identifier 10.1109/TIP.2012.2205007
laser images. In this paper, we focus on the multiplicative noise η that follows the Gamma distribution (see [9], [11]) with the probability density function K K η K −1 −K η e (2) (K ) where K is an integer related to the noise level and is the classical Gamma function. In the literature, several methods have been proposed for multiplicative noise removal. Most of these methods can be classified into two categories: the methods derived from the classical MAP Bayesian rule (see [9], [11]), and the efficient methods that convert the multiplicative noise removal problem into an additive noise removal problem (see [22]). The following is a brief review of these two classes of methods. pdf(η) =
A. MAP-Based Multiplicative Noise Removal The image denoising problem consists in estimating the original image u from the observed image f . A classical statistical approach based on the maximization of the posterior distribution Prob(u| f ) leads to the MAP estimate uˆ = arg max Prob(u| f ). u
Equivalently, uˆ is a minimizer of − log Prob(u| f ), and using Bayes rule eventually leads to the variational problem: f γ inf logu + + φ(u) u>0 u Z where we assume that the image u follows the Gibbs prior γ (u) = Z1 e−γ φ(u) with the normalization factor Z . Inspired by the above MAP analysis, several variational methods have been proposed for Gamma multiplicative noise removal problems. A rather popular one is the Aujol-Aubert (“AA” ) method [9], which estimates the original image u through the following minimization problem: f + λuT V . inf log u + (3) u>0 u Here, λ is the positive regularization parameter and uT V is the Total Variation (TV) regularization term. In a continuous setting, and for a smooth image u, it reads |∇u(ξ )|dξ uT V :=
where |∇u| = (∂u/∂ξ1 )2 + (∂u/∂ξ2 )2 . The main drawback of the objective function in (3) is that it is not a convex function of u. The optimal solution
1057–7149/$31.00 © 2012 IEEE
HUANG et al.: MULTIPLICATIVE NOISE REMOVAL VIA A LEARNED DICTIONARY
4535
of (3) may not be unique. Moreover, classical optimization algorithms provide only a local minimum of the objective function, which may not necessarily be the global optimal solution and depends on initial guesses of iterative optimization methods. Inspired by the “AA” model, Huang et al. [11] reformulated the problem in the logarithm domain (log u) and obtained a convex optimization problem. Then they designed an alternating minimization method to solve this optimization problem. Numerical results showed that this convex model can provide denoised images with higher SNRs and better visual quality, independently of initial guesses. Recently, in [10], a split-Bregman algorithm was presented to solve (3) directly and it is numerically quite efficient. In [23], a simple and fast method using graph-cut minimization was studied.
The main goal of this paper is to propose an effective method for multiplicative noise removal which combines the MAP formulation, the sparse representation principle, and the idea of converting the multiplicative noise removal into the additive noise removal. The outline of the paper is as follows. In Section II, we first review the K-SVD method which is currently one of the best performing ones for additive noise removal. In Section III, a new model is proposed for multiplicative noise removal, and some basic properties of it are discussed in Section IV. In Section V, we show experimental results to demonstrate the quality of the denoised images and the efficiency of the proposed method. Finally, concluding remarks are given in Section VI.
B. Methods Based on Additive Noise Removal
In this section, let us briefly review the K-SVD√method √ in a discrete setting, that is, when is a grid of size N × N . This method was originally proposed by Aharon and Elad [4] for additive Gaussian √ √ noise removal. Suppose that a√noisy N × N indexed by = {1, 2, . . . , N }2 , image f of R written as a column vector f ∈ R N , results from a zeromean Gaussian noise b ∈ R N with a known standard deviation σ superimposed on an original image u ∈ R N . The basic assumption of √ the K-SVD method is that each image patch √ (of fixed size n × n, and rewritten as a column vector in Rn ) can be represented sparsely as a linear combination of atoms taken from a fixed dictionary D ∈ Rn×k . Here dictionary means a set of column vectors in Rn . Each column vector is called an atom and is usually of unit norm. The total number of these atoms is k. A typical dictionary may be an overcomplete DCT dictionary (sampling cosine waves in different frequencies to obtain a fixed number of atoms, see [4], [26] for details), or may be learned from a noisy image (see below). Let us assume that √ √ (6) P = {1, 2, . . . , N − n + 1}2 √ √ is the set of n × n image patches located in . Given a vector x = (x 1 , x 2 , . . . , x m ) ∈ Rm , we shall consider its classical l 2 norm x, and also its l 0 “norm”
II. A DDITIVE K-SVD D ENOISING A LGORITHM Since there have been a number of good methods to remove additive noise in the literature, it is natural to consider converting the multiplicative noise denoising problem into the additive noise denoising problem. The main idea is to take the logarithm of both sides of equation (1), leading to log f = logu + logη.
(4)
This is an additive noise removal problem, for which the probability distribution of the noise logη is given (see, e.g., [9]) by K K e K (logη−e pdf(logη) = (K )
logη )
.
(5)
In [14], Shi and Osher extended the inverse scale-space method to multiplicative noise removal models by converting multiplicative noise into additive noise (4), and their proposed model is strictly convex. Furthermore, the convergence of the flow for the multiplicative noise model was demonstrated, as well as its regularization effect and its relation to the Bregman distance. Durand et al. [22] proposed an efficient hybrid approach based on the log-image data and a reasonable underoptimal hard-thresholding on its curvelet transform. They minimized a criterion composed of an l 1 data-fidelity term (on the thresholded frame coefficients) and a TV regularization term in the log-image domain. The main interest of their method is that it combines the advantages of shrinkage and variational methods. Theoretical results on the hybrid criterion were presented, and a Douglas-Rachford splitting algorithm was used to minimize the energy function. The numerical results demonstrated the superiority to the other compared methods, especially for images containing tricky geometrical structures. Moreover, Steidl and Teuber [15] considered using the so called I-divergence as a data fitting term, and TV or nonlocal means as a regularization term in the cost function for removing multiplicative Gamma noise. A nonlocal filter for removing multiplicative noise was proposed by Teuber et al. [16]. Furthermore, Teuber et al. presented a new similarity measure for nonlocal filters in [24] and the denoising results revealed some good properties of the filters. Some other nonlocal-based filters are further elaborated in [25].
x0 := #{i |1 ≤ i ≤ m, x i = 0}, that counts the number of nonzero entries in x 1 , . . . x m . Under the sparsity assumption, the Gaussian noise removal task can be described as the following energy minimization task: ˆ u} {αˆ i j , D, ˆ = arg min λ f − u2 {D,αi j ,u}
+
1 Dαi j − Ri j u||2 + μi j αi j 0 . 2 (i, j )∈P
(i, j )∈P
(7) Rn×N
In this model, Ri j ∈√ √is a binary n × N matrix which extracts from u a n × n patch located in (i, j ) ∈ P (thus, Ri j u ∈ Rn ). The first term of (7) demands a proximity between the recorded image f and its denoised unknown version u. The second term demands that each patch Ri j u of the reconstructed image u can be represented by the dictionary
4536
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 21, NO. 11, NOVEMBER 2012
D ∈ Rn×k by coefficients vector αi j ∈ Rk with a bounded error. The third part is a (weighted) sparsity constraint, which demands that the number of coefficients required to represent any patch is small. Note that the values μi j are patch-specific weights that are determined hiddenly by the optimization procedure. The minimization (7) with respect to the unknowns D, αi j , u yields the denoising algorithm. The choice of D is of high importance for the performance of the algorithm. In [4], it is shown that dictionary training can be done by the minimization problem (7). The algorithm proposed there is an iterative block-coordinate relaxation method which fixes all the unknowns except for the one to be updated. Thus, (7) is solved with a three-step iterative method: learning typical patterns (dictionary) from the noisy image via K-SVD; representing patches sparsely via the Orthogonal Matching Pursuit (OMP) algorithm [27]; averaging the denoised patches to finally obtain a global denoised image. III. P ROPOSED M ODEL Inspired by K-SVD [4] for dictionary learning and the “AA” model [9], here we propose a log-l 0 minimization approach, that is, an approach based on the sparse representation (l 0 ) of the log-image over a certain unknown dictionary of image patches. Denoting by√ 1 the √ constant 1 over the discrete image domain (a N × N grid) and by ·, · the usual scalar product in Euclidean spaces, the proposed model reads f ˆ αˆ i j , u} { D, ˆ = arg min λ log u + , 1 + γ loguT V u {D,αi j ,u>0} 1 + Dαi j − Ri j log u||2 2 (i, j )∈P + μi j αi j 0 (8) (i, j )∈P
where λ, γ are positive regularization parameters, P is given in (6) and u ∈ R N is the estimated image. The · T V term is defined, in the discrete setting, by summing over the image domain the norm of ∇u, the classical 2-neighbors discrete gradient estimate (see, e.g., [1]). As in (7), Ri, j ∈ Rn×N is the matrix corresponding to the extraction of the patch located in (i, j ), and αi, j ∈ Rk is the sparse vector of coefficients to represent the patch Ri, j log u with the dictionary D ∈ Rn×k . The hidden parameters (μi j )(i, j )∈P are determined by the optimization procedure described in [4]. The interpretation of the new model (8) is rather clear. The first to a constant, is the discrete version of the term term, up f (log u+ ) u in model (3); the second term is inspired by [11], [14], [15] where the non-convex issue of (3) is avoided; the last two terms assume that the small patches of the transformed image log u can be represented sparsely over a certain adaptive dictionary D. As TV is efficient to recover cartoon images while dictionaries are well adapted to textures, we can expect that Equation (8), which defines a kind of “cartoon + texture” model, will lead to good results. Similarly to [4], we propose a two-phases algorithm to approximate this log-l 0 minimization problem:
1) Learn the dictionary D from the log image and represent the log image patches via the learned dictionary to get Dαi, j . 2) Use a fast iterative method to solve a reduced convex model to get the recovered image. A. Learning a Dictionary from the Degraded Image and Obtaining its Sparse Representation In the first step of our method, a dictionary D which contains typical patterns of the observed image log f is obtained through the K-SVD method, and each patch (i, j ) of log f is represented sparsely with the coefficient vector αi, j in this dictionary [4], [11]. This is just the OMP and K-SVD steps in the additive K-SVD denoising algorithm, and the procedure is as follows. For the set of image patches Ri, j log f , we seek a dictionary D ∈ Rn×k and αi, j that minimize the cost function 1 {D, αi, j } = arg min Ri, j log f − Dαi, j 2 {D,αi, j } 2 (i, j )∈P + μi, j αi, j 0 , (9) (i, j )∈P
where μi, j is the hidden parameter which indicates how those two forces should be weighted. An iterative algorithm is designed to solve (9) effectively [4]. Here we note that in the proposed model, the sparse representation is performed in the logarithm domain of the image. Since the log function is strictly monotone, f and log f have the same edges, so we can expect that the efficiency of the K-SVD method remains on log f [11]. B. Dictionary and Total Variation Based Multiplicative Noise Removal Using the sparse dictionary representation Dαi, j obtained for each patch Ri, j log f via K-SVD in the first phase, we now minimize (8) with respect to u, that is, f arg min λ log u + , 1
u u>0 1 + Ri, j log u − Dαi, j 2 + γ loguT V . 2 (i, j )∈P
(10) As we shall see later, the experimental results show that the proposed model (10) is effective. However, the criterion is not convex with respect to u, and thus the global minimizer is difficult to obtain. Inspired by [11], [14], let us introduce the transformation u = exp(y). Then we can obtain a new equivalent model, which is now convex with respect to the variable y: arg min λ y + f e−y , 1
y
+
1 Ri, j y − Dαi, j 2 + γ yT V . (11) 2 (i, j )∈P
Using the notations RiTj Ri j , W = (i, j )∈P
M=
(i, j )∈P
RiTj Dαi, j
HUANG et al.: MULTIPLICATIVE NOISE REMOVAL VIA A LEARNED DICTIONARY
(note that W and M have the same dimensions as u and f ), Equation (11) can be rewritten as arg min λ y + f e y
−y
, 1
1 + W y, y − M, y + γ yT V . (12) 2 Let us denote by i, j the coefficient (i, j ) of any matrix . We can notice that since for any (i, j ) ∈ , the weight Wi, j counts how many times pixel (i, j ) is used to construct √ the √ image patches of size n × n, we have 1 ≤ Wi, j ≤ n. Proposition 1: Assume that inf f > 0, then the objective function in (12) is strictly convex. Proof: The function t ∈ R → t + f i, j e−t is convex, hence so is y → y + f e−y , 1 ; y → − M, y is also convex; y → W y, y is strictly convex since Wi, j ≥ 1; and classically y → yT V is convex. As a sum of convex and strictly convex functions, the objective function (12) is strictly convex. Thanks to the Proposition above, the minimization problem (12) can be handled by many typical optimization algorithms. Below we will design an efficient dual approach to solve it. Denote X = R N . The gradient ∇u is a vector in the vector space Y = X × X. And for p = ( p1 , p2 ), q = (q 1 , q 2 ) ∈ Y , we define a scalar product by pi,1 j qi,1 j + pi,2 j qi,2 j . p, q Y = √ 1≤i, j ≤ N
We also consider the discrete divergence operator div : Y → X, which is the opposite of the adjoint of the discrete gradient operator [1]. Indeed, we have − div = ∇ ∗ , that is, ∀u ∈ X, ∀ p ∈ Y, ∇u, p Y = − u, div p X . Notice that classically, we have yT V = ∇ y1 = max ∇ y, p
p∞ ≤1
where similarly to [28], we define p∞ , p1 as | pi, j | p∞ = max | pi, j |, p1 = i, j
4537
Algorithm 1 Chambolle-Pock Algorithm Solving (12) 1: Initialization: y¯ (0) , y (0), p(0) . 2: For k = 0, 1, 2, . . ., repeat until a stopping criterion is reached p(k+1) := argmin p∞ ≤1 {− p, ∇ y¯ (k) + β p − p(k) 2 }, y (k+1) := argmin y∈X {G(y) − div p(k+1) , y
+βy − y (k) 2 }, y¯ (k+1) := 2y (k+1) − y (k) , 3: Output: yˆapprox := y (k+1) .
formula, while the second step can be handled with Newton’s method efficiently. Theoretically, if the positive parameter β is chosen big enough, then the above algorithm converges to the solution of (13). √ Proposition 2: If β > 2, then Algorithm 1 converges to the unique minimizer of (12). Proof: The problem (12) can be rewritten as min G(y) + ∇ y1 . y∈X
This is exactly the nonlinear primal problem studied in [28] [see (3) there] and one can carefully verify that the Program 1 here is the Algorithm 1 proposed there. By Theorem 1 of [28], the algorithm converges as soon as 1 2 1 · L 2 since L ≤ 2 2 (see [1], [28]). The functional proposed in (10) is a dictionary + TV-based multiplicative noise removal model. The dictionary approach permits to restore textures efficiently, while TV is known to preserve edges well. As we shall see, the proposed model combines both merits together and manages to efficiently restore images corrupted by a multiplicative Gamma noise.
i, j
IV. BASIC P ROPERTIES
with | pi, j | = ( pi,1 j )2 + ( pi,2 j )2 . Now let us consider 1 1 G(y) = λ y + f e−y , 1 + W y, y − M, y . γ 2 We can easily prove that the dual problem of (12) is max min G(y) − div p, y .
p∞ ≤1 y∈X
(13)
Inspired by [28], [29], an iterative splitting method can be applied to solve (13). Starting from an initialization y¯ (0) , y (0), p(0) the sequence of iterates y (1), p(1) , y (2) , p(2) , · · · can be obtained from Program 1. Note that in this algorithm, the first step can be updated by means of a closed-form
Before we present numerical results, let us study the basic properties of the model (10) rewritten in the continuous settings, that is, f log u + λ inf log u∈B V (),u>0 u 1 + Rs log u − Dαs 2 2 s∈I
+ γ loguT V
(14)
where (s )s∈I is a finite set of small patches covering , and Rs is the restriction on s . We need a simple lemma. Lemma 3: Let a > 0 and M0 = min{0, 2 + 2 log( a2 )}. Then ∀x ∈ R, x + ae−x ≥ |x| + M0 .
4538
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 21, NO. 11, NOVEMBER 2012
Proof: Set g(x) = x + ae−x . Then g(x) − x ≥ 0 for any x ∈ R, and g(x) + x = 2x + ae−x attains its minimal value 2 + 2 log a2 in x = log a2 . We thus have ∀x ∈ R, g(x) − |x| = min(g(x) − x, g(x) + x) ≥ M0 as announced. Proposition 4: Let be a bounded open subset of R2 with Lipschitz boundary, and f a bounded function such that inf f > 0. For u satisfying log u ∈ BV (), u > 0, let f E(u) := log uT V + λ log u + . u Then E(u) has a unique minimizer. Proof: Let v = log u, then E(u) = F(v) with F(v) := vT V + λJ (v) and J (v) := (v + f e−v ). By Lemma 3, we have (15) J (v) ≥ (v + (inf f )e−v ) ≥ v1 + M0 ||
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
for some constant M0 (depending on f ). Thus, F is convex and bounded from below. Consider a minimization sequence (v n ) for F. Then v n T V and J (v n ) are bounded from above. By (15), so is v n 1 . Hence, (v n ) is bounded in BV (). As BV is compact in L 1 (), there is a v ∈ BV such that a subsequence (v nk ) converges to v in L 1 , and we may assume (up to another extraction) that v nk → v pointwise almost everywhere. As the BV norm is l.s.c., vT V ≤ lim inf v nk T V . Since v nk + f e−v nk is bounded from below, by Fatou’s Lemma, we have J (v) ≤ lim inf J (v nk ). Thus F(v) ≤ lim inf F(v nk ) = inf B V () F, and v minimizes F. Consequently, E admits at least one minimizer, u = ev . Readily, F is a strictly convex function. Hence F (and consequently E) admits a unique minimizer. Theorem 5: Let be a bounded, open subset of R2 with Lipschitz boundary. Let f be a positive, bounded function. For u > 0 satisfying log u ∈ BV (), let E(u) : = γ log uT V f 1 +λ Rs u log u + + u 2 s∈I
−Dαs
2
(16)
where γ > 0. Suppose that (s )s∈I is a finite local coverage 2 of and s∈I Dαs < +∞. Then E(u) has a unique minimizer. Proof: Consider a minimization sequence (u n ) of E. Let v n = log u n . Similarly to the previous proposition, one can take a subsequence (v nk ) which converges both in L 1 and BV . V. N UMERICAL R ESULTS In this section, the denoising results are presented to demonstrate the efficiency of the proposed method. Comparisons are also made with the most efficient methods that have been proposed recently for multiplicative noise removal. These five methods are: the “BS” algorithm [30]; the “AA” method [9]; the “DFN” algorithm [22]; the “DRSM” method [15] and the updating nonlocal method “NL” [16]. The parameters of all
Fig. 1. Comparative results for the Cameraman image (noise level: K = 10). (a) Original. (b) Noisy (K = 10). (c) AA (24.48). (d) BS (24.43). (e) DFN (26.08). (f) DRSM (25.93). (g) NL (26.55). (h) Proposed (27.32).
(a)
(b)
(c)
(d)
(e)
(f)
Fig. 2. Sample of the recovered camera stand part of Cameraman corresponding to Fig. 1 when K = 10. (a) Original. (b) BS (24.43). (c) DFN (26.08). (d) DRSM (25.93). (e) NL (26.55). (f) Proposed (27.32).
the methods are selected with respect to the best recovered PSNR and MAE values (see definitions below). Except for the experiments appearing in [22], the parameters used in the “AA,” “BS” and “DFN” methods are the same as in this paper. As in [22], two quantities are to assess the quality √ used √ of the recovered images. For a N -by- N pixels noise-free image u and its denoised version uˆ by a given algorithm, we define the peak signal-to-noise ratio (PSNR, in decibels) by PSNR = 10 log10
N · |max(u) − min(u)|2 (in dB) uˆ − u2
where |max(u)−min(u)| measures the gray-scale range of the original image. The Mean Absolute-deviation Error (MAE)
HUANG et al.: MULTIPLICATIVE NOISE REMOVAL VIA A LEARNED DICTIONARY
(a)
(b)
(d)
4539
(a)
(b)
(c)
(d)
(e)
(f)
(c)
(e)
(f)
Fig. 5. Sample of the recovered yard part of Nîmes corresponding to Fig. 4 when K = 10. (a) Original. (b) BS (27.0). (c) DFN (27.80). (d) DRSM (27.95). (e) NL (28.42). (f) Proposed (28.85).
Fig. 3. Sample of the recovered camera stand part of Cameraman when K = 4. (a) Original. (b) BS (22.36). (c) DFN (22.98). (d) DRSM (23.67). (e) NL (24.20). (f) Proposed (24.91).
(a)
(a)
(b)
(c)
(d)
(e)
(f)
(b) Fig. 6. Sample of the recovered yard part of Nîmes when K = 4. (a) Original. (b) BS (24.92). (c) DFN (25.84). (d) DRSM (25.75). (e) NL (26.07). (f) Proposed (26.34).
(c)
(d)
(e)
(a)
(f)
(g)
(b)
(h)
Fig. 4. Comparative results for the Nîmes image (noise level: K = 10). (a) Original. (b) Noisy (K = 10). (c) AA (25.26). (d) BS (27.0). (e) DFN (27.80). (f) DRSM (27.95). (g) NL (28.42). (h) Proposed (28.85).
(c)
(d)
(e)
(f)
(g)
(h)
is defined by 1 uˆ − u1 . N The images Cameraman with size 256-by-256, Nîmes (a French city) with size 512-by-512, Peppers with size 256by-256 and Barbara (head part) with size 256-by-256 are used in our experiments. The images are contaminated by multiplicative Gamma noise with two different levels K = 10 and K = 4 respectively. The original images are shown in (a) of Fig. 1, 4, 7, 10, respectively. The corresponding corrupted images are shown in (b) of Fig. 1, 4, 7, 10, respectively, and the corresponding recovered images by the “AA,” “BS,” “DFN,” “DRSM,” “NL” and the proposed method are shown MAE =
Fig. 7. Comparative results for the Peppers image (noise level: K = 10). (a) Original. (b) Noisy (K = 10). (c) AA (24.40). (d) BS (25.09). (e) DFN (26.40). (f) DRSM (26.52). (g) NL (27.07). (h) Proposed (27.30).
in subparts (c) to (h) of these figures respectively. Certain corresponding zoomed samples of the recovered images by the “BS,” “DFN,” “DRSM,” “NL” and the proposed method
4540
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 21, NO. 11, NOVEMBER 2012
(a)
(b)
(d)
(c)
(e)
(a)
(b)
(c)
(d)
(e)
(f)
(f)
Fig. 8. Sample of the recovered base part of Peppers corresponding to Fig. 7 when K = 10. (a) Original. (b) BS (25.09). (c) DFN (26.40). (d) DRSM (26.52). (e) NL (27.07). (f) Proposed (27.30).
(a)
(b)
(c)
(d)
(e)
(f)
Fig. 9. Sample of the recovered base part of Peppers when K = 4. (a) Original. (b) BS (22.82). (c) DFN (23.81). (d) DRSM (24.19). (e) NL (24.32). (f) Proposed (24.78).
(a)
Fig. 11. Sample of the recovered shoulder part of Barbara corresponding to Fig. 10 when K = 10. (a) Original. (b) BS (22.55). (c) DFN (23.80). (d) DRSM (23.06). (e) NL (24.60). (f) Proposed (25.47).
(a)
(b)
(c)
(d)
(e)
(f)
(b)
Fig. 12. Sample of the recovered shoulder part of Barbara when K = 4. (a) Original. (b) BS (20.97). (c) DFN (21.60). (d) DRSM (21.63). (e) NL (22.70). (f) Proposed (23.03).
(c)
(d)
(e)
(f)
(g)
(h)
Fig. 10. Comparative results for the Barbara image (noise level: K = 10). (a) Original. (b) Noisy (K = 10). (c) AA (22.56). (d) BS (22.55). (e) DFN (23.80). (f) DRSM (23.06). (g) NL (24.60). (h) Proposed (25.47).
are shown in subparts (b) to (f) of Fig. 2, 3, 5, 6, 8, 9, 11, 12, respectively. The recovered PSNR values are provided in the subpart captions of each figure.
Since the readability of the whole recovered images is rather poor, we only show them in the case K = 10 and essentially base our observations on close-ups. Moreover, due to the inferior performance of the “AA” method, we do not show the corresponding images for this method in the close-up figures. Tables I to IV report two performance assessments (the PSNR and MAE measurements, with the best results in boldface) for each of the six considered methods, along with the computation time (in seconds) reported in the “TIME” column. For the proposed method, the computing time is expressed as the sum of the dictionary training time and the remaining computation time. Both in terms of visual effect and PSNR/MAE measurements, it clearly appears that the proposed method outperforms the five others.
HUANG et al.: MULTIPLICATIVE NOISE REMOVAL VIA A LEARNED DICTIONARY
TABLE I P ERFORMANCE E VALUATION FOR THE Cameraman I MAGE , T IME FOR
P ROPOSED M ETHOD = D ICTIONARY T RAINING + R EMAINING
4541
TABLE III P ERFORMANCE E VALUATION FOR THE Peppers I MAGE , T IME FOR P ROPOSED M ETHOD = D ICTIONARY T RAINING + R EMAINING
Recovered K
10
4
Method
Recovered
PSNR
MAE
TIME(s)
AA BS DFN
24.48 24.43 26.08
9.90 9.44 7.40
19.35 1.00 58.10
DRSM NL Proposed
25.93 26.55 27.32
7.24 6.82 6.63
AA BS DFN
21.97 22.36 22.98
DRSM NL Proposed
23.67 24.20 24.91
K
Method
PSNR
MAE
TIME (s)
AA BS DFN
24.40 25.09 26.40
10.91 9.90 8.49
21.25 0.99 54.17
5.62 4.67 16.93 + 2.90
DRSM NL Proposed
26.52 27.07 27.30
8.13 7.62 7.57
5.96 5.74 13.60 + 2.93
14.26 12.14 10.61
19.58 0.99 53.74
AA BS DFN
22.68 22.82 23.81
13.67 13.08 11.52
20.41 1.00 52.67
9.23 8.97 8.97
5.65 4.57 11.55 + 2.94
DRSM NL Proposed
24.19 24.32 24.78
10.47 10.15 10.14
6.11 5.72 10.66 + 1.61
10
4
TABLE II P ERFORMANCE E VALUATION FOR THE Nîmes I MAGE , T IME FOR
TABLE IV P ERFORMANCE E VALUATION FOR THE Barbara I MAGE , T IME FOR
P ROPOSED M ETHOD = D ICTIONARY T RAINING + R EMAINING
P ROPOSED M ETHOD = D ICTIONARY T RAINING + R EMAINING
Recovered K
10
4
Method
Recovered
PSNR
MAE
TIME (s)
AA BS DFN
25.26 27.0 27.80
8.83 7.69 7.21
113.0 3.9 243.2
DRSM NL Proposed
27.95 28.42 28.85
6.96 6.59 6.33
AA BS DFN
24.55 24.92 25.84
DRSM NL Proposed
25.75 26.07 26.34
K
Method
PSNR
MAE
TIME (s)
AA BS DFN
22.56 22.55 23.80
12.82 12.38 10.86
15.45 1.06 56.30
29.46 19.25 30.12 + 9.60
DRSM NL Proposed
23.06 24.60 25.47
11.65 9.82 8.99
6.79 3.95 16.39 + 1.24
10.06 9.87 9.09
120.4 3.7 240.3
AA BS DFN
18.69 20.97 21.60
20.05 15.16 14.21
15.31 0.92 55.04
8.93 8.75 8.60
29.49 18.07 18.31 + 9.43
DRSM NL Proposed
21.63 22.70 23.03
13.87 12.34 12.04
6.90 2.37 13.17 + 1.71
For the Cameraman image with K = 10 (see Fig. 1, Fig. 2 and Table I), the estimated PSNR for “DRSM” and “NL” methods are 1.39 and 0.77 smaller than with the proposed method, while the MAE are 0.61 and 0.19 larger. In the case K = 4 (again, for Cameraman image), the best PSNR and MAE scores are obtained with the proposed method too. Comparing the visual appearance of the recovered images in Fig. 2 and Fig. 3, we find that the “NL” method yields a good result too, but the camera stand is better recovered and the artifacts are smaller with the proposed method. We also observe that the “AA” method over-smooths the image and sometimes leaves outliers in the recovered images; “BS” and “DFN” give rise to more artifacts in the recovered images; the recovered effect obtained by “DRSM” is not as good as the results of the “NL” and the proposed method. For the larger size picture Nîmes, Fig. 5 and Fig. 6 clearly show that the proposed method can preserve some details in the yard part of the images, while the same part is more obscure in the images recovered with the other methods.
10
4
For the picture Peppers, the recovered results shown in Fig. 7, Fig. 8 and Fig. 9 display the same recovered qualities as for Cameraman and Nîmes images: less artifacts appear in the recovered images and the textures are best preserved with the proposed method. Concerning Barbara image (zoom on head part, as in [15]), Fig. 10, Fig. 11 and Fig. 12 further show the superiority of the proposed method to recover the texture of the original image. Fig. 13 show the recovered effect when TV is not used in the proposed method [γ is zero in (8)] for the restoration of Barbara image in the case K = 4. Clearly, the result obtained with TV has less artifacts than without TV, and this is especially clear for the flat region. Therefore, the combination of TV regularization and dictionary learning improves the efficiency of the proposed method to a certain degree. All the recovered results show that the proposed method can effectively remove the multiplicative noise in the degraded images, and it manages to preserve well both edges and textures at the same time.
4542
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 21, NO. 11, NOVEMBER 2012
R EFERENCES
(a)
(b)
(c)
Fig. 13. Sample of the recovered hand part of Barbara when K = 4, with and without the TV term in (8). (a) Original. (b) Without TV (22.65). (c) Proposed (23.03).
Fig. 14.
Dictionary learned for the (log) image Cameraman when K = 10.
Concerning the results obtained with the method proposed in [11], as expected they are not sensitive to the initialization of the iterative algortihm used, but we observed that they were just a little better than those obtained by the “AA” method. Also, the method proposed in [14] always takes a long time to complete the restorations. Both methods are not as good as the methods “DFN,” “DRSM,” “NL” and the proposed method. Note that √ in all experiments, the parameter λ was set to be λ = 1/ (1, K ), where is the logarithmic derivative of the Gamma function, β was set to 10, and γ = 1. The proposed method is much less sensitive to the parameter values. To illustrate the mechanism of the proposed method, we show the dictionary learned from the logarithm of the Cameraman image (K = 10) in Fig. 14. VI. C ONCLUSION In this paper, we considered the multiplicative noise removal problem. Writing a MAP formulation for the log-transformed image with a combination of a TV-regularization and a sparse representation in an adaptive dictionary of image patches, we obtained a new variational formulation that combines the edge-preserving efficiency of TV-regularization with the ability of dictionary-based methods to recover textures well. Using the Orthogonal Matching Pursuit and an adaptation of Chambolle-Pock algorithm, we built an algorithm that permits to efficiently remove multiplicative noise and outperforms the most recent ones, both visually and in quantitative terms. ACKNOWLEDGMENT The authors would like to thank S. Durand, J. Fadili, M. Nikolova, and Dr. Teuber, for kindly providing the datum and original codes of their papers. The authors would also like to thank the anonymous reviewers for their extremely helpful comments.
[1] A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imag. Vision, vol. 20, nos. 1–2, pp. 89–97, 2004. [2] T. Chan and J. Shen, Image Processing and Analysis: Variational, PDE, Wavelet, and Stochastic Methods. Philadelphia, PA: SIAM, 2005. [3] L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Phys. D, vol. 60, nos. 1–4, pp. 259–268, 1992. [4] M. Elad and M. Aharon, “Image denoising via sparse and redundant representations over learned dictionaries,” IEEE Trans. Image Process., vol. 15, no. 12, pp. 3736–3745, Dec. 2006. [5] S. Lintner and F. Malgouyres, “Solving a variational image restoration model which involves L ∞ constraints,” Inverse Problems, vol. 20, no. 3, pp. 815–831, 2004. [6] D. Donoho and I. Johnstone, “Ideal spatial adaptation via wavelet shrinkage,” Biometrika, vol. 81, no. 3, pp. 425–455, 1994. [7] J. Starck, E. Candes, and D. Donoho, “The curvelet transform for image denoising,” IEEE Trans. Image Process., vol. 11, no. 6, pp. 670–684, Jun. 2002. [8] A. Buades, B. Coll, and J. Morel, “A review of image denoising algorithms, with a new one,” J. Multi. Model. Simul., vol. 4, no. 2, pp. 490–530, 2005. [9] G. Aubert and J. Aujol, “A variational approach to remove multiplicative noise,” SIAM J. Appl. Math., vol. 68, no. 4, pp. 925–946, 2008. [10] J. Bioucas-Dias and M. Figueiredo, “Multiplicative noise removal using variable splitting and constrained optimization,” IEEE Trans. Image Process., vol. 19, no. 7, pp. 1720–1730, Jul. 2010. [11] Y. Huang, M. Ng, and Y. Wen, “A new total variation method for multiplicative noise removal,” SIAM J. Imag. Sci., vol. 2, no. 1, pp. 20–40, 2009. [12] K. Krissian, C. Westin, R. Kikinis, and K. Vosburgh, “Oriented speckle reducing anisotropic diffusion,” IEEE Trans. Image Process., vol. 16, no. 5, pp. 1412–1424, May 2007. [13] L. Rudin, P. Lions, and S. Osher, “Multiplicative denoising and deblurring: Theory and algorithms,” in Geometric Level Sets in Imaging, Vision, and Graphics. New York: Springer-Verlag, 2003, pp. 103–119. [14] J. Shi and S. Osher, “A nonlinear inverse scale method for a convex multiplicative noise model,” SIAM J. Imag. Sci., vol. 1, no. 3, pp. 294– 321, 2008. [15] G. Steidl and T. Teuber, “Removing multiplicative noise by DouglasRachford splitting methods,” J. Math. Imag. Vision, vol. 36, no. 2, pp. 168–184, 2010. [16] T. Teuber and A. Lang, “Nonlocal filters for removing multiplicative noise,” in Scale Space and Variational Methods. New York: SpringerVerlag, 2012, pp. 50–61. [17] Y. Yu and S. T. Acton, “Speckle reducing anisotropic diffusion,” IEEE Trans. Image Process., vol. 11, no. 11, pp. 1260–1270, Nov. 2002. [18] V. Dutt and J. Greenleaf, “Adaptive speckle reduction filter for logcompressed B-scan images,” IEEE Trans. Med. Imag., vol. 15, no. 3, pp. 802–813, Dec. 1996. [19] M. Isana, R. Wagner, B. Garra, D. Brown, and T. Shawker, “Analysis of ultrasound image texture via generalized rician statistics,” Opt. Eng., vol. 25, no. 6, pp. 743–748, 1986. [20] R. Wagner, S. Smith, and J. Sandrik, “Statistics of speckle in ultrasound B-scans,” IEEE Trans. Sonics Ultrason., vol. 30, no. 3, pp. 156–163, May 1983. [21] A. Kryvanos, J. Hesser, and G. Steidl, “Nonlinear image restoration methods for marker extraction in 3d fluorescent microscopy,” Proc. SPIE, vol. 5674, pp. 432–443, Mar. 2005. [22] S. Durand, J. Fadili, and M. Nikolova, “Multiplicative noise removal using L 1 fidelity on frame coefficients,” J. Math. Imaging Vision, vol. 36, no. 3, pp. 201–226, 2010. [23] L. Denis, F. Tupin, J. Darbon, and M. Sigelle, “SAR image regularization with fast approximate discrete minimization,” IEEE Trans. Image Process., vol. 18, no. 7, pp. 1588–1600, Jul. 2009. [24] T. Teuber and A. Lang, “A new similarity measure for nonlocal filtering in the presence of multiplicative noise,” to be published. [25] C. A. Deledalle, L. Denis, and F. Tupin, “Iterative weighted maximum likelihood denoising with probabilistic patch-based weights,” IEEE Trans. Image Process., vol. 18, no. 12, pp. 2661–2672, Dec. 2009. [26] F. Malgouyres and T. Zeng, “A predual proximal point algorithm solving a non negative basis pursuit denoising model,” Int. J. Comput. Vision, vol. 83, no. 3, pp. 294–311, 2009. [27] Y. Pati, R. Rezaiifar, and P. Krishnaprasad, “Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition,” in Proc. 27th Asilomar Conf. Signals Syst. Comput., vol. 1. 1993, pp. 629–639.
HUANG et al.: MULTIPLICATIVE NOISE REMOVAL VIA A LEARNED DICTIONARY
[28] A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging,” J. Math. Imag. Vision, vol. 40, no. 1, pp. 120–145, 2011. [29] M. Zhu and T. Chan, “An efficient primal-dual hybrid gradient algorithm for total variation image restoration,” Dept. Math., UCLA, Los Angeles, CAM Rep. 08-34, 2008. [30] C. Chesneau, M. Fadili, and J. Starck, “Stein block thresholding for image denoising,” Appl. Comput. Harmon. Anal., vol. 28, no. 1, pp. 67–88, 2010.
Yu-Mei Huang received the M.Sc. degree in computer science from Lanzhou University, Lanzhou, China, in 2000, and the Ph.D. degree in computational mathematics from the Hong Kong Baptist University, Kowloon, Hong Kong. She is currently an Associate Professor with Lanzhou University. Her current research interests include signal and image processing, scientific computing, and numerical linear algebra.
Lionel Moisan received the Ph.D. degree from Paris-Dauphine University, Paris, France, and the HDR degree from Paris XI University, in 1997 and 2003, respectively. He studied mathematics and computer science from the École Normale Supérieure, Paris. He was a CNRS Researcher with École Normale Supérieure de Cachan, Bruz, France, from 1998 to 2003, and is currently a Professor of applied mathematics with Paris-Descartes University. His current research interests include variational models for image processing (restoration, smoothing) and a-contrario statistical models for image analysis (gestalt and structure detection).
4543
Michael K. Ng received the B.Sc. and M.Phil. degrees from the University of Hong Kong, Hong Kong, in 1990 and 1992, respectively, and the Ph.D. degree from the Chinese University of Hong Kong, Hong Kong, in 1995. He was a Research Fellow with the Computer Sciences Laboratory, Australian National University from 1995 to 1997, and an Assistant and Associate Professor with the University of Hong Kong from 1997 to 2005. He has been a Professor with the Department of Mathematics, Hong Kong Baptist University since 2005. His current research interests include bioinformatics, image processing, scientific computing, and data mining. Dr. Ng serves on the editorial boards of international journals.
Tieyong Zeng received the B.S. degree from Peking University, Beijing, China, the M.S. degree from École Polytechnique, Palaiseau, France, and the Ph.D. degree from the University of Paris XIII, Paris, France, in 2000, 2004, and 2007, respectively. He was a Post-Doctoral Researcher with the Centre de Mathématiques et de Leurs Applications, ENS de Cachan, France. He is currently an Assistant Professor with Hong Kong Baptist University, Hong Kong. His current research interests include image processing, statistical learning, and scientific computing.
