A Cartoon-Texture Decomposition Based Multiplicative Noise Removal Method

We propose a new frame for multiplicative noise removal. To improve the multiplicative denoising performance, we add the regularization of texture component in the denoising model, designing a multiscale multiplicative noise removal model. The proposed model is jointly convex and can be easily solved by optimization algorithms. We introduce Douglas-Rachford splitting method to solve the proposed model. In the algorithm, we make full use of some important proximity operators, which have closed expression or can be executed in one time iteration. In particular, the proximity of norm is deduced, which is just the Fourier domain filtering. In the process of simulation experiments, we first analyze and select the needed parameters and then test the experiments on several images using the designed algorithm and the given parameters. Finally, we compare the denoising performance of the proposed model with the existing models, in which the signal to noise ratio (SNR) and the peak signal to noise ratios (PSNRs) are applied to evaluate the noise suppressing effects. Experimental results demonstrate that the designed algorithms can solve the model perfectly and the recovery images of the proposed model have higher SNRs/PSNRs and better visual quality.


Introduction
Image denoising is a basic and important task in image processing.The relatively mature developed denoising model is the additive noise model in which case the noise is assumed to obey a Gaussian distribution, that is,  =  + ,  ∼  (0,  2 ) . (1) However, the noises involved in many applications do not conform to the characteristic of the additive one; they may corrupt an image in other forms.In this paper, we are concerned with the denoising problem under the assumption that the original image has been corrupted by some multiplicative noise.The corresponding examples include the uneven phenomena during the magnetic resonance imaging (MRI) and the speckle noise in the ultrasonic and in synthetic aperture radar (SAR) image.The degradation model of the multiplicative noise can be expressed as where ⋅ refers to the componentwise multiplication.To model the actual problems, we assume that the multiplicative noise obeys some random distribution; for example, the noise in SAR images is assumed to obey the Gamma distribution and the one in ultrasonic image obeys the Rayleigh distribution.

Multiplicative Noise Removal
where the mean of noise is assumed to be 1 and the variance is assumed to be  2 .Since (3) is a nonconvex problem, it is very difficult to be solved.The second classic model is called AA model [2] which is given by Aubert and Aujol under the hypothesis of Gamma distribution: The objective function in ( 4) is also nonconvex.However, Aubert and Aujol showed the existence of minimizers of the objective function and employed a gradient method to solve (4) numerically.There are several improved works based on (4) developed in recent years.
Recently, Shi and Osher [3] proposed considering a logarithm transformation on the noisy observation, log  = log+log , and derived the TV minimization model for multiplicative noise removal problems.Huang et al. also proposed the log-domain denoising model using the transformation  = exp() in (4).Consider which is known as EXP model [4].The objective function in ( 5) is strictly convex in  and the authors promoted an alternating minimization algorithm to solve the model and showed the convergence of the method simultaneously.Nevertheless, model ( 5) is convex in  rather than in  in the original image domain.At the same time, Durand et al. [5] proposed a method composed of several stages.They also used the log-image data and applied reasonable suboptimal hard thresholding on its curvelet transform; then they applied a variational method by minimizing a specialized hybrid criterion composed of an  1 data-fidelity term of the thresholded curvelet coefficients and a TV regularization term in the logimage domain.The restored image can be obtained by using an exponential of the minimizer, weighted in such a way that the mean of the original image is preserved.Their restored images combine the advantages of shrinkage and variational methods.Besides the above approaches, dictionary learning methods and nonlocal mean methods have also been proposed and developed for multiplicative denoising [6][7][8][9].In [10], Zhao et al. developed a convex optimization model for multiplicative noise removal.The main idea is to rewrite a multiplicative noise equation such that both the image variable and the noise variable are decoupled.That is to say, rewrite problem (2) as where  = diag() is a diagonal matrix of which the main diagonal entries are given by []  .According to (2), when there is no noise in the observed image, we obtain that  = , a vector of all ones.When there is a multiplicative noise in the observed image, we expect that {[]  0 ̸ = 0, ∀}, and moreover they are greater than zeros.So we can say that  is invertible, and ( 6) is equivalent to where  = diag() is the diagonalizable matrix of vector  and In (8), the first term is to measure the variance of , the second term is the data term, and the third term is the TV regularization term.If the first term is absent and  can be arbitrarily assigned, then minimizing min In [11], it was pointed out that minimizing (10) can extract the large scale component and leave out the small scale contents of image .In other words, we can extract cartoon component with different scales based on the selection of parameter  =  2 / 1 .We verify the above conclusion through the following experiment.The example is shown in Figure 1(a), an image  with four squares with size of 3 × 3, 5 × 5, 20 × 20, and 80 × 80 pixels, respectively.The experiment results are the following: when  = 0.01, all the four squares can be extracted, as shown in Figure 1(b); when  = 0.05, three squares can be extracted while the one of 3×3 is lost, as shown in Figure 1(c); when  = 0.2, the two big squares appear, and another two are lost, as shown in Figure 1(d); when  = 2, only the biggest square 80×80 remains, while all other ones are lost, as shown in Figure 1(e).The above results are consistent with the theoretical analysis in [11]; the bigger  is, the bigger the scale of extraction is.As a result, we have sound reasons to believe that the minimizer  of model (10) is almost the cartoon component of the whole restored image we expect.
The above analysis shows that, if we want to recover more contents from an image , we should choose a small , which bring much more small scales and details.However, for a noisy image, the small scale component may contain much more noise, which may hinder the denoising quality.For this reason, we will improve the denoising model through adding the prior information of the texture component and design a cartoon-texture decomposition based denoising model.Then, we design the numerical algorithm based on the Douglas-Rachford splitting algorithm and three useful proximity operators.Before conducting the experiment, we analyze the selection of the corresponding parameters.Finally, we verify the performance of our proposed model through comparing the SNRs/PSNRs indexes with some existing models.
The rest of the paper is organized as follows.We propose our new model in next section.In Section 3, the numerical method is presented, in which the Douglas-Rachford (DR) splitting algorithm and some proximal operators used in this paper are firstly reviewed, and the algorithm is then designed and described in detail.In Section 4, the denoising experiments on two types of multiplicative noise are implemented, and the experiment results and performance analyses are unfolded.Finally, we conclude our work in Section 5.

The Proposed Model
Based on the analysis in Section 1.2, the recovery image obtained by solving (8) mainly contains the cartoon component of what we expect, while most of the significant texture components are lost.To optimize the denoising performance, we improve model ( 8) by adding the texture information and alter the regularization term ‖‖ TV to be ‖‖ TV , where  denotes the cartoon component of the restored image .
Taking into consideration the prior information of texture V =  − , we select ‖V‖ 2  −1 as the regularization term about V, for that minimizing ‖⋅‖ 2  −1 can extract the texture component very well.In a word, the new model is and the restored image is  * =  * +V * .In (11),  1 ,  2 , and  are the positive regularization parameters to control the balance among the terms in the objective function.
The proposed model possesses the following superiorities.Firstly, we add the term ‖V‖ 2  −1 in the regularization part, which can help to recover more image information and improve the quality of the restored image.Secondly,  2 / 1 and  can be adjusted according to the noise level, which makes the model have multiscale property.Thirdly, when  → +∞, then V → 0, and (11) degrades to model (8); that is, (11) is the modified and improved version of (8).Finally, the proposed model keeps the structure of (8), is jointly convex for (, , V), and can be easily solved by optimal methods.In the following section, we use alternative iteration method and Douglas-Rachford method to deal with it, which combine the application of several proximity operators.

The Numerical Method
3.1.Douglas-Rachford Algorithm and Some Proximity Operators.Before proposing the numerical scheme of ( 11), we will review some basic algorithms and operators.The first one is the Douglas-Rachford splitting algorithm.Definition 1.Let Φ() and Ψ() be the proper, l.s.c., and convex functions such that and Ψ() + Φ() → +∞ as ‖‖ → +∞; then the problem admits at least one solution and, for any  ∈ (0, +∞), its solutions are characterized by Algorithm 2.
Secondly, we review the definition of proximity operator and present several special examples that will be used in the numerical scheme.Definition 3. Let Φ be a proper, l.s.c., and convex function; for every  ∈   , the minimization problem admits a unique solution, which is denoted by prox Φ (), and the operator prox Φ () :   →   thus defined is called the proximity operator of Φ.
In the following, we list three examples of proximity operator which will be used in this paper: (1) Φ() = ‖  ‖ 1 ,  = ( 1 ,  2 , . . .,   )  ∈   .Set  = prox Φ (); then  can be componentwise expressed in the close form: which is known as the soft thresholding operator, and we denote it in short by  = ST  (). ( There are many strategies to solve it.One of the most classical methods is the semi-implicit gradient descent algorithm proposed by Chambolle in [12]; another is the Forward-Backward method developed in [5].Both of the above methods need inner loop in the processing of numerical scheme.In this paper, we adopt the variable splitting method and express the problem in discrete form: where ∇ is the gradient operator given by ∇ = (∇ (1) , ∇ (2) ) , [∇ ()]  = ([∇ (1) ()]  , [∇ (2) ()]  ) ,  = 1, 2, . . ., , We introduce an auxiliary variable  ∈  ×2 to transfer ∇: with a sufficiently large penalty parameter .We solve (19) by an alternating minimization scheme given in Algorithm 4.
Algorithm 4 (proximity of TV).Consider the following: (1) Fixing , compute  componentwise It is worth mentioning that the above alternating minimization scheme is embedded into the whole algorithm with just one time iteration, which makes the algorithm have no inner iteration.
(3) Consider Φ() = ‖‖ 2  −1 .We will deduce prox Φ () in the continuous situation, with the purpose of describing the derivation process clearly; this does not affect the processing of the numerical implementation.The proximity of Φ() is expressed as Denote the objective function in (20) by while minimizing (V) is equivalent to minimizing Ĥ(V) in Fourier domain: minimizing Ĥ(V) in V yields the unique solution V = Lx, where Taking inverse Fourier transform, we have It is indicated from (24) that the proximity of  −1 norm is just the Fourier domain filtering.

Solving the Proposed Model
Numerically.We solve (11) by alternatively updating , , V in turn.The framework is given as in Algorithm 1.
In Algorithm 1, updating  is to compute the proximity of ‖ ⋅ ‖ 1 , which is just the soft thresholding operator and can be realized by (15).Updating  is to solve a TV-L1 problem, and there are many methods to deal with it; to improve the algorithm efficiency, we realize the updating using Douglas-Rachford splitting algorithm, combining the computing of prox ‖⋅‖ 1 and prox ‖⋅‖ TV .Updating V is to solve an L1- −1 problem, which is also dealt with by Douglas-Rachford splitting algorithm, combining the soft thresholding and Fourier domain filtering.A detailed description of the whole process is shown as in Algorithm 2.
There are some issues that need to be illustrated.

Implementation Details and Experimental Results
This section is mainly devoted to numerical simulation of image restoration in the presence of multiplicative noise.We test the performance of our model under the corruption of two types of multiplicative noise: Gamma and Rayleigh.
Gamma Noise.The probability density function of Gamma noise  is given by where  > 0,  > 1.The mean and the variance of  are As multiplicative noise of an image, the mean of  is set to be 1, so  = 1/, and the variance of  is  2 = 1/.In model (10),  = {1/()}  =1 , and the value of its mean is estimated as [10] Rayleigh Noise.The probability density function of Rayleigh noise  is given by where  is a positive parameter.The mean value of  is equal to √/2, the variance of  is equal to (2 − /2)/ 2 , and the We test experiments on three images: Cameraman (size of 256 × 256), Lena (size of 512 × 512), and Aerial image of some city (size of 512 × 512); see Figure 2.
In this paper, we use the signal to noise ratio (SNR) and the peak signal to noise ratio (PSNR) as the evaluation indicators.These indicators are measured between clean and restored images.Let  new represent the image restored from the noisy image  noisy , and let  2  0 represent the average variance of the clean image  0 .Then, we define the SNR indicator by the formula as follows: ) , and the PSNR is defined as ) , where  denotes the length of the images.
4.1.Parameters Selection.Some necessary parameters  1 ,  2 , ,  1 ,  1 ,  1 , and  2 are needed to be given to start up the new algorithm;  1 ,  2 , and  are regularization parameters to keep balance among the terms in model (11).In particular,  =  2 / 1 is an important indicator to trade off between the data term and the regularization term. 1 and  1 are two relaxation parameters of the Douglas-Rachford splitting algorithms in Algorithm 2 and are empirically set to be 0.5. 1 and  2 are the parameters of the proximity operator of ‖ ⋅ ‖ TV and ‖ ⋅ ‖ 2  −1 , respectively, and are empirically set to be 10.
According to the experimental analysis in [10], the ratio  =  2 / 1 in ( 8) is almost a constant depending on the type of noise and was empirically estimated to be 5/6, 5/8 for Gamma  noise and Rayleigh noise, respectively.From a number of experiments, we find that the best value of  2 / 1 can be chosen from the range [0.6, 0.85], and the exact selection depends on the different noise levels and images.We give the values of  2 / 1 in Table 1.
For , it is a key parameter to promote the performance of model (11).Since a proper value of  is critical to the recovery of a satisfying image, it is necessary for us to emphatically discuss its impact on the denoising results.Here, we show its influences on image decomposition by using the image Cameraman.Fixing other parameters, we decompose the image into cartoon component and texture component with different  and put the experimental results in Figure 3.It is clear from Figure 3 that the different selections of  can bring completely different results.If  is selected too large ( ≥ 100), there are very few contents in the texture V, and many details are found in .With the decrease of , more and more features get involved in the texture component V.
However, when  is selected too little (e.g.,  ≤ 0.1), the block effect in cartoon component  is very serious, and some major features are included in the texture component V.Moreover, in the process of denoising, the too little value of  may bring much more details or small scale parts appearing in V, including some noise we do not expect.For this reason, we experimentally select  in the interval [20, 80] according to the different noise level, and the exact selection can be found in Table 1.

Experimental Results in the Assumption of Gamma Noise.
In the assumption of Gamma noise corruption, we test denoising performances of AA model, convex model, and the proposed model on the images of Figure 2    which is solved by the gradient descent method where Δ is step size and is set as 0.01 for all experiments and  is the regularization parameter and is adjusted according to the images and noise level; the exact value can be found in Table 1.The convex model is solved by the ADMM method; see [10].
In the test, the original images are corrupted by Gamma noise with  = 4,  = 10, and  = 15, respectively.In Table 1, we list values of related parameters and obtained SNRs/ PSNRs by taking average of ten noise cases under the same experimental setting.It is clear that the proposed model performs quite better in terms of PSNR values and SNR values than AA model and convex model.
We further show the denoised results obtained by three methods in Figures 4-6.It is observed that when  = 15 and  = 10, as shown in Figures 4 and 5, the restored images from the new model possess more clear contents and much better visual effects than the AA model and the convex model.However, when  = 4, that is, the noise level is higher, we can observe from Figure 6 that the new model possesses little advantage over the convex model in the visual effect.We also observe from Table 1 that, when  = 4, the SNRs/PSNRs of our model and the convex model are neck to neck.The main reason for this phenomenon is that the larger the variance of the noise is, the more the texture components are seriously destroyed and regarded as noise part, and the bigger value  should be adjusted.As described in the analysis (c) of the new model, there will be a little content in the V parts, and denoising results of the new model will approximate that of the convex model.

Experimental Results in the Assumption of Rayleigh Noise.
The Rayleigh noise is generated by √−2 log(1 − ), where  is a uniformly distributed random variable generated by "rand" in MATLAB.The mean of Rayleigh noise is set to be 1, and the variance is 0.2732.Since AA model is deduced under the assumption of Gamma noise, we just test the denoising performance of the convex model and the proposed model.Table 2 shows the SNRs/PSNRs results of the two models.we notice that the proposed model gets higher SNRs and PSNRs than the convex model.The experiment results are shown as in Figure 7.It can be observed that the restored images from the new model possess more clear contents and much better visual effects.

Further Analysis of the Restored Images.
In this section, we further analyze the cartoon part  and the texture part V of the restored image  of model (11).Take the case under Gamma noise as an example; we exhibit the experiment results when  = 10.
It can be observed from Figure 8 that there are rich contents in the texture parts V, which enhance the image quality in both the SNRs/PSNRs index and the visual effects.So, the proposed model can also be regarded as a cartoontexture decomposition method under the multiplicative noise corruption, on which we need further research.

Conclusions
A novel cartoon-texture decomposition based multiplicative noise removal model has been proposed in this paper.The main advantage of our work is embodied in the following two aspects: the one in which we take most of the a priori information of texture component into the denoising model, which makes the model possess more flexibility and efficiency.Another is that we dexterously solve the new convex model using optimal splitting algorithms and Fourier domain filtering method.The experiment results demonstrate that the proposed model can better remove the noise than the two other models, and the recovery images have higher SNRs/PSNRs and better visual effects.

Figure 5 :
Figure 5: Experiment results under the corruption of Gamma noise,  = 10.

Figure 6 :
Figure 6: Experiment results under the corruption of Gamma noise,  = 4.

Figure 7 :
Figure 7: Experiment results under the corruption of Rayleigh noise with the mean 1 and the variance 0.2732.

Figure 8 :
Figure 8: Experiment results of the proposed model.From (a-l), (a, e, i) are the noisy images corrupted by Gamma noise ( = 10), (b, f, j) are the restored cartoon parts, (c, g, k) are the restored texture parts, and (d, h, l) are the whole recovery images.
[1]els.There are many multiplicative noise removal models; one of the most classic models based on TV regularization is called RLO model proposed by Rudin et al.[1]: 2, . . ., ;

Table 1 :
PSNR and SNR results for the Gamma noise removal.

Table 2 :
PSNR and SNR results for the Rayleigh noise removal.