A Convex Variational Approach for Image Deblurring With Multiplicative Structured Noise

The restoration of images corrupted by blurring and structured noise has attracted growing attention in the domains of image processing and computer vision. However, many works only focus on the restoration of the images degraded by blurring and additive structured noise or multiplicative structured noise separately. It is still a challenge and an open problem to restore degraded images with blurring and multiplicative structured noise, simultaneously. In this paper, based on the total variation (TV), the statistical property of the Gamma noise and the maximum a posteriori (MAP) estimator, we obtain a convex variational model to recover blurred images with multiplicative structured noise. Especially, to get this convex model, we reformulate the prior assumption of the images degradation model by division instead of multiplication. For solving this convex model, an effective alternating direction method of multipliers (ADMM) is employed. Numerical experiments are presented to illustrate the effectiveness and efficiency of the proposed model.


I. INTRODUCTION
Image deblurring and denoising are fundamental tasks in image processing and computer version. In real applications, image blurring is unavoidable due to many factors such as long transmission distance, environmental electromagnetic interference, and air turbulence [1]. Moreover, the observed images can be further affected by both blurring and noise interference in the process of image acquisition and transmission [2], [3]. Generally, we can consider two kinds of random noises: additive and multiplicative noise. Under the additive noise scheme, numerous deblurring methods have been proposed to handle these problems. We refer the readers to see [4]- [15]. Indeed, various regularization approaches (e.g. the wavelet frame [5], [14], [15] and the total variation (TV) [4], [6], [7], [11]) were used to recover good numerical results from the additive noise. Under the multiplicative noise scheme, the blurred image Hu was further corrupted by some The associate editor coordinating the review of this manuscript and approving it for publication was Emanuele Crisostomi . multiplicative noise η. The degradation model is as follows, where u ∈ R m×n is the original image; H is the blurring operator and u 0 is the observed image; and η ∈ R m×n is a random vector that could follow one or more statistical distributions such as Gamma distribution and Gaussian distribution. Let E denote the vector space of 2D images defined on = {1, . . . , m} × {1, . . . , n}. For any 2D images, the total number of pixels is m × n. The pixels of the image are identified by a multi-index i = (i 1 , i 2 ) ∈ . The point-wise product between two elements Hu and η is denoted Hu η, Unlike the additive noise with independent properties, multiplicative noise is generally coupled with the original images. Though this standard multiplicative noise assumption can model some real applications such as synthetic aperture radar (SAR), ultrasound imaging, and laser images [16], it still lost many image information seriously, especially the structural information. In coherent imaging systems, the photographs that we observed have distinct structural features. For instance, we could observe the ripple effect in X-ray imaging [17], [18], the speckle noise in ultrasonic imaging [19]- [21], and the stripes effect in medical equipment imaging [22]. In this paper, we mainly consider the restoration of blurred images under multiplicative structured noise. To the best of our knowledge, this problem has not yet been studied. To take the structure information into account, similarly to [23], [24], we use a stationary noise η to replace the standard assumption describing the structural information of the noise. We can generate the structured noise η by convolving the noise ν with a known kernel, where the operator ψ is a convolution filter that depends on the noise structure, the symbol is the convolution product, and the convolution product η = ψ ν is defined by where the noise ν is a random vector. It can be seen that structured noise η depends on the structure ψ and the noise ν. Until the past decade, a few variational methods have been proposed to handle the problem of images degraded by blurring and multiplicative noise. When the blur operator H = Id in (1) (i.e. Id is an identity operator), these problems become the multiplicative noise removal problems. Rudin et al. [25] firstly proposed the variation model to solve this issue (called ''RLO model''), where σ 2 is the noise variance, and u TV is the total variation of u to preserve the shape edges (see Section II-A for details).
The RLO model only uses the basic statistical properties, the mean and the variance, which limits the image restoration effect. Based on the Gamma noise with the mean of 1, Aubert and Aujol [26] used a maximum a posteriori (MAP) estimator to obtain a multiplicative denoising model (called ''AA model''), where α is a parameter that balances the data fidelity term and the TV term. Although the AA model is non-convex, Aubert and Aujol proved that the objective equation has an optimal solution. The non-convexity of the RLO and AA model poses serious numerical troubles. Shi and Osher [27] proposed a globally convex total variational framework (called ''SO model''). They transformed the multiplicative noise into the additive noise through logarithmic transformation and used the spatial method to remove multiplicative noise. Huang et al. [28] maintained the framework of the AA model by using a transformation w = log u, this model is strictly convex and the alternating direction method of multipliers (ADMM) was proposed to find the optimal solution of the objective function. The numerical results show that the restoration effect of this method is better than the AA model. Steidl and Teuber [29] proposed a convex variational model based on I-divergences to recover images corrupted by multiplicative noise. Dong et al. [30] proposed the nonlocal total variational model to restore multiplicative noise and used the split Bregman method to deal with this minimization problem. Inspired by the Gaussian adaptive method for denoising, Chen and Cheng [31] developed a multiplier denoising model with local constraints of TV regularization term, which can automatically select regular parameters by using the statistical properties of random variables. Except for TV regularization methods, there are also many other ways to remove multiplicative noise well such as dictionary learning, tight-frame approach and so on. We refer the reader to see [32]- [39].
In recent years, many scholars have investigated the problem of image deblurring with multiplicative noise. The classical RLO and AA model can be extended to restore blurred images, min u u TV , and To overcome the drawback of the non-convex model in (7) and (8), Dong and Zeng [40] proposed a convex denoising and deblurring model by introducing a quadratic penalty term into the AA model. In [41], Wang and Ng transformed the multiplicative noise and blur removal problem into the additive denoising and deblurring problem by logarithmic transformation and used the ADMM algorithm to solve it. Dong and Zeng [42] used an I-divergence technique to build a convex model. Zhao et al. [43] proposed a model containing the total variational regular term, the data fidelity term, and the variance term, and obtained analytical solutions by ADMM algorithm. Shama et al. [44] proposed a convex model, which contains the data fidelity term based on MAP, the penalty term based on total generalized variational regularization and noise statistical characteristics, separately. There are also other models and algorithms focusing on this topic [7], [16], [45]- [52].
In order to get higher quality images, the problem of restoring images with structured noise has been getting more and more attention in the fields of space, biological sciences, and medical imaging. Münch et al. [53] combined the wavelet and Fourier transform to remove stripes and ring artifacts. Boas and Fleischmann [54] summarized the algorithms and restoration models for recovering computed tomography (CT) imaging. Henrik Fitschen et al. [55] used a variational model with directional first and second order differences to recover stripes noise. Roels et al. [56] analyzed the main factors of image degradation and proposed a non-local means image restoration algorithm to restore images under the 3D electron microscope. Chen et al. [57] proposed a method based on TV regularization and group sparsity constraint to deal with the fringe noise of remote sensing images. Fehrenbach et al. [24] (see Section II-B for details) proposed a variational stationary noise removal method to remove the additive structured noise in the scanning electron microscope. In [23], the authors proposed a variational model based on Gamma distribution to restore multiplicative structured noise (see Section II-C for details). This model is convex and provides good recovering results. Han et al. proposed a novel algorithm based on a new definition of similarity coefficient to remove PolSAR image speckle [58]. Yang et al. used the Schatten 1/2-norm regularization to the remote sensing images destriping [59]. However, all of the above methods only focus on the restoration of the images degraded by blurring and additive structured noise or multiplicative structured noise separately. As far as we know, it is still a challenge and an open problem to restore degraded images with blurring and multiplicative structured noise, simultaneously.
In this paper, inspired by the pioneer works [23], [24], we study the deblurring issues with multiplicative structured noise. Comparing with previous works, our contributions are as follows.
• First, we propose a convex variational model to restore the images degraded by blurring and multiplicative structured noise.
• Second, we adopt the standard ADMM algorithm to solve the proposed model, and this algorithm has a better effect and speed.
• Third, numerical experiments show that our model and algorithm have good numerical performance and are very robust in recovering images under different blur kernels, different types of structured noise and different Gamma noise levels. The rest of this paper is organized as follows. In Section 2, we review the total variation regularization and the related works of [23], [24]. We propose our model to restore blurred images with multiplicative structured noise in Section 3. In Section 4, the ADMM algorithm is employed to solve the proposed model. In Section 5, we demonstrate the effectiveness and efficiency of the proposed model by some numerical experiments. Conclusions are drawn in Section 6.

II. RELATED WORKS A. TOTAL VARIATION REGULARIZATION
In [60], Rudin et al. firstly introduced the total variation u TV in image restoration to preserve sharp edges. The discrete total variation of u can be defined by [16] where j = 1, · · · , m, k = 1, · · · , n, u j,k is the (j, k)th pixel location of the image u. And (∇u) x j,k and (∇u) y j,k denote the horizontal and vertical first-order difference, respectively,

B. ADDITIVE STRUCTURED NOISE REMOVE
Let us review the previous works in [23], [24] briefly about generating and solving structured noise. In [24], Fehrenbach et al. used a stationary noise assumption to replace the white noise assumption. Based on this new additive noise assumption, the additive degradation model can be rewritten as: The operator ψ i is a known elementary pattern, to represent the image structural information. The vector ν i follows white noise with known probability density functions. The model was proposed to restore images affected by additive structured noise as follows, where the first term is TV regularization term; the second term is data fitting term which depends on the probability distribution of ν i , and φ i (ν i ) is Tikhonov regularization [61]. They exploited primal-dual algorithms to solve this model, and obtained good numerical results for restoring image under additive structured noise.

C. MULTIPLICATIVE STRUCTURED NOISE REMOVE
In [23], Escande et al. proposed a convex model to restore images with multiplicative structured noise. Here, the degradation model can be rewritten as where δ = ψ λ, and λ is noise that follows Gamma distribution. The point-wise division between two elements u and δ is denoted Based on this degradation model, a convex variation model to recover images affected by multiplicative structured noise is proposed as follows, where c is the parameter to balance the TV term and the data fidelity term; the parameters (a and b) are from Gamma distribution. This model can handle the problem of images with multiplicative structured noise well. However, the blurring effect is unavoidable and not taken into account in [23]. Thereby, in this paper, we consider the deblurring issues under multiplicative structured noise.

III. THE PROPOSED MODEL A. IMAGING PIPELINE
The classical blurred images with multiplicative noise degradation model is like (1) u 0 = Hu η. VOLUME 8, 2020 Based on this degraded function, many approaches were proposed to deal with blurring issues. However, this equation leads to serious numerical troubles since the MAP estimator is non-convex, like the AA model. To avoid the numerical difficulties, similarly to [23], let u denote the original image and u 0 denote the observed image, λ is the multiplicative noise, and H is a blurred operator, we change the degradation model as follows, where δ = ψ λ, the vector λ follows Gamma distribution, and its density function is where τ (α) is the Gamma function. The parameter a > 0 is called the shape parameter, the parameter b > 0 is called the inverse scale parameter. The advantages of (14) are as follows:  • We rewrite multiplication into the division, it makes u and λ become a linear system to overcome numerical troubles because of using MAP estimator.
• In this degradation model, we consider the factors of the blurring, multiplicative noise, and structure. The model can more accurately reflect coherent imaging systems, such as SAR, Single Plane Illumination Microscope (SPIM), and Atomic Force Microscope (AFM) imaging.

B. THE PROPOSED MODEL VIA THE MAP ESTIMATOR
In this section, based on (14), we propose a convex model to restore blurred images with multiplicative structured noise. VOLUME 8, 2020  We assume that P(x) is the probability of variable x and P(x|y) is the probability of x on observation y. Our target is to maximize P(u, λ |u 0 ). Using Bayes rules, we have By taking the negative logarithmic transformation, this problem can be rewritten as follows, − log(P(u 0 |u, λ )) − log P(u, λ) + log P(u 0 ).
Because u and λ are independent, we can get P(u, λ) = P(u)P(λ). Moreover, the clean image u and the observed image u 0 satisfy the linear relation. It means that P(u 0 |u, λ) and P(u 0 ) are constants [24]. In this paper, we assume that images have a low total variation, which can be expressed as As λ follows Gamma distribution, we can get the data fitting term. By introducing the constraint condition: u 0 = Hu (ψ λ), we can get the restoration model finally, The proposed model (19) is convex if we can keep the parameters α > 1, b > 0. By introducing auxiliary variable y and z, model (19) can be rewritten as

IV. NUMERICAL METHOD
Let's start a general convex minimization model with separable structures, where f (x 1 ), g(x 2 ), and q(x 3 ) are lower semicontinuous proper convex function; A i are given matrices; b is a known vector; X i are nonempty closed convex sets. Hence, the convex optimization model (20) falls into the above form with the following specifications: , z), and the abstract sets are X i := R n ; • the matrices A i and b are given by, respectively: where I is the identity operator. The augmented Lagrangian function of the model (20) is given by where ω 1 , ω 2 , and ω 3 are the Lagrangian multipliers, µ 1 , µ 2 , and µ 3 are the penalty parameters. The augmented Lagrangian method (ALM) for problem (21) is an iterative algorithm based on the iteration, ω k+1 We apply the alternating direction method of multipliers (ADMM) to the minimization of L (z, u, y, λ; ω 1 , ω 2 , ω 3 ) in (21), (23d)

A. Z SUBPROBLEM
For z subproblem, we can use the shrinkage operator to solve directly, where the shrinkage operator shrink ·, c µ 1 is defined by Especially, the element ξ |ξ | should be taken as 0 if |ξ | = 0.

B. U SUBPROBLEM
For fixed z, y, λ, the minimization of L with respect to u can be rewritten as follows, where Under the periodic boundary condition for u, we know that ∇ T ∇ and H T H are block circulation, so we can use the fast VOLUME 8, 2020  Fourier transform (FFT) to find the solution of this subproblem efficiently, whereM u = F(µ 1 ∇ T ∇ + µ 2 H T H )F .

C. λ SUBPROBLEM
For this subproblem, this problem can be rewritten as then this question can be solved by where

D. Y SUBPROBLEM
For y subproblem, this problem is equal to We use Preconditioned Conjugate Gradients Method (termed ''PCG'') to solve (32). This method is based on the conjugate gradient method, which can improve the speed of the algorithm by preprocessing. Assuming that the preconditioned matrix P is symmetric and positive definite, the linear equations AY = B can be transformed into If Q −1 A has smaller condition number than A, it can improve the speed of solving this system of linear equations. Using this method, the condition is that A is symmetric and positive definite, and B is known. In this subproblem, we have Because ψ T ψ as symmetric and positive definite, µ 3 I is a diagonal matrix, we know that A is symmetric and positive definite, and B is known. It means that the y subproblem can be solved by this method. We use MATLAB built-in ''PCG'' function to solve, We set TOL as 1e −5 ; MAXIT is the maximum number of iterations (i.e., we set MAXIT is 40). The overall algorithm is summarized in Algorithm 1.

V. NUMERICAL EXPERIMENTS
In this section, we use four images which are given in Figure 1   Core (TM) i7-8565U CPU at 1.80 GHz with 8 GB of RAM. Quantitatively, the quality of the recovered images is evaluated by the peak signal-to-noise ratio (PSNR), the relative error (ReErr), and structural similarity (SSIM) [62]. The PSNR and ReErr are defined as, where u 0 and u are the restored image and the original image, respectively. The size of these images is m×n. And the SSIM [62] can be defined as follows, where µ x , µ y are the averages of x, y; σ 2 x , σ 2 y are the variances; σ xy is the covariance of x, y, and θ 1 , θ 2 are two variables to stabilize the division with a weak denominator. The overall SSIM (MSSIM [63]) is the mean of local similarity indexes defined by where x i , y i are corresponding windows of clean/restored image indexed by i and N . Note that the larger PSNR, VOLUME 8, 2020  larger SSIM, and smaller ReErr means the better restoration effect. To illustrate the effectiveness of the proposed method, we compare the recovering results of our approach with those of the following image deblurring algorithms, • AA model is proposed by Aubert and Aujol [26], and solved by the gradient method; • HNZ model is proposed by Huang et al. [16], and solved by an alternating minimization algorithm; • DZ model is proposed by Dong and Zeng [40], and solved by primal-dual algorithm; • WZN model is proposed by Wang et al. [36], and solved by the ADMM algorithm; • LY model is proposed by Lu et al. [50], and solved by an alternating minimization algorithm.
The HNZ model and DZ model are extensions of the AA model. The WZN model is a framelet-based convex optimization model for multiplicative noise and blur removal problems. Wang, Zhao, and Ng proposed a new method to select the regularization parameter by using the l 1 norm-based L-curve method for these framelet-based models. The LY model is a variational model and this model utilizes the   favorable properties of framelet regularization. These models can be used to solve the blurred images with multiplicative noise, and make competitive results for restoring blurred images with multiplicative Gamma noise. Without loss of generality, like the DZ model, we also assume that the mean of λ equals 1 under the multiplicative noise scheme. Let K = α, 1/θ = b, then we have that K θ = 1 and its variance is 1/K . The density function (15) can be rewritten as In our degradation model, we mainly consider three factors: Gamma noise levels, blur kernels, and structured noise types. Firstly, we use the MATLAB built-in function to generate the Gamma noise, • random (''gamrnd'', K, 1/K), where the level of the Gamma noise is determined by the parameter K. Note that the smaller K is, the more serious the Gamma noise is. In all experiments, we choose different K numbers from the set of K ∈ {0.05, 0.15, 0.5, 1.5, 3, 5, 33}. For blur kernels, we consider three blur kernels (''Gaussian VOLUME 8, 2020   For GB, we test three different sizes (5 × 2, 5 × 5, 7 × 2); as for MB, we test three experiments. We test len = 3, 5 when the theta is the default number (i.e. theta = 0); when the theta = 30, the len is 5; about AB, we test hsize = 3, 5. Sure, in addition to these three common blur kernels, there are other blur kernels, such as ''disk blur'', ''laplacian blur'', ''Sobel blur'' and so on.
For structured noise types, we choose four types from six common structured noise types (see Figure 2), point, speckle, stripes, and grid noise, respectively. Speckle, stripes, and grid noise are common structured noise, point noise is multiplicative Gamma noise with no structure. As shown in Figure 3, the first row is the original images; the second row is the images only affected by multiplicative structured noise; the third row is the images only affected by blurring; the last row is the  In this paper, we use the ADMM algorithm to solve our model. About the regularization parameter c and the ADMM penalty parameters µ 1 , µ 2 , µ 3 , we choose the best one from the set of {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 50}. And it is based on different types of noise. After testing, we found that when c = 2, µ 1 = µ 2 = µ 3 = 8, our algorithm has a good recovery effect on the selected test set. In all the experiments, the stop criterion is either −4 or the maximum iterative number exceeds 500.
For the other methods (the AA, HNZ, DZ, WZN, and LY model), we choose the corresponding parameters such that the MSSIM result is highest. VOLUME 8, 2020

B. EXPERIMENT 1
In the first experiment, the original images (including ''peppers'', ''f16'', and ''zelda'') are blurred by three different blur kernels, and further degraded by multiplicative structured noise. We choose stripes noise as structured noise in this experiment. And the K value (Gamma noise level) we set is 1 or 0.15. Tables 1-3 display the value of ReErr, MSSIM, and PSNR for three blurring kernels with different parameters. Table 4 lists the CPU-time of these six models. In Figure 4, we show the restored results of the Gaussian blur (GB(5,2)/K = 1). Figure 5 shows the partial restored effect of motion blur (MB(5)/K = 3). And we display the restoration results for average blur (AB(5)/K = 0.5) in Figure 6. All these results demonstrate our model can remove different blur kernels with structured noise effectively and also confirm that our proposed model gets better results than the other five models.
C. EXPERIMENT 2 In this experiment, we show the recovering results of the images degraded by different types of structured noise. The original image is corrupted by three different noise types

Algorithm 1 ADMM for the Proposed Model
Input: the observed image u 0 , psi, and the blur operation H . Parameters: ω 1 , ω 2 , ω 3 , µ 1 , µ 2 , µ 3 , and stopping parameters tol. including speckle noise, grid noise, and point noise (point noise without structure), and then this image is also blurred by MB (5). The degraded images and the recovered images are shown in Figures 7-9, respectively. The values of ReErr, MSSIM, and PSNR are shown in Table 5. These recovering results demonstrate that the proposed model is very effective to recover blurred images with different types of structured noise.

D. EXPERIMENT 3
In this experiment, we consider showing our model to restore images with different Gamma noise levels. We let ''peppers'' as our test image. Except for the noise levels, the image ''peppers'' is also corrupted by GB (5, 2) and stripes noise in this experiment. For Gamma noise, the smaller the K value is, the image degradation effect is more serious. We choose three different K numbers (K = 0.1, 0.5, 1.5). Table 6 tells that the proposed model improves the quality of restoration in terms of ReErr, MSSIM, and PSNR compared with the other models, and the speed is also better than the WZN, LY, DZ, AA, and HNZ models. Figures 10-12 show the degraded images with different K values and the corresponding recovered images, respectively. E. EXPERIMENT 4 In Experiment 4, we consider a simplified model. When H = Id, it means that our model can be used to recover images that are only affected by multiplicative structured noise. Then, the ADMM algorithm in this paper can be compared with the primal-dual algorithm proposed by literature [23] (called the ''EWZ'' model), and the AA, DZ, WZN, and LY model. We choose the stripes noise as the structured noise, the K value is 1. The results are as follows. From Figure 13 and Table 7, although the DZ model takes the shortest time, our restoration effect is better significantly. In general, our model and algorithm are better than others.
F. EXPERIMENT 5 In the aforementioned four experiments, we mainly show the restoring effect on synthetic images. In this experiment, we show our model to recover the real image. Due to the complexity of the real image and the numerous influencing factors, we can only conduct qualitative analysis and compare them through visual effects. For the real image, it is important to estimate the blur kernel, and then use the no-blind method to estimate latent image. In this paper, we mainly propose a non-blind approach to restore the blurred image with multiplicative structured noise. In this experiment, we use the general Gaussian kernel to approach this real image blur kernel. There are two reasons to do this estimation. First, the Gaussian blur is widely used in the estimation of static VOLUME 8, 2020 blur in natural imaging [64]- [66]; the second reason is that we use a blurring kernel estimation method proposed by [45] to estimate roughly. The kernel is shown as Figure 14(b), which is very much like the Gaussian blur height. The more accurate blur kernel estimation will be studied in our future work. Based on this kernel estimation, the restoration result is shown in Figure 14. It can be found that the recovery effect of our model is better.

VI. CONCLUSION
In this paper, we discuss the image restoration problems degraded by blurring and multiplicative structured noise, simultaneously. The most relevant conclusions are summarized as follows, • A convex model based on the total variation (TV), the statistical property of the Gamma noise and the maximum a posteriori (MAP) estimator was proposed to restore degraded images with blurring and multiplicative structured noise, simultaneously.
• We employed an effective ADMM algorithm to handle this minimization problem, which is based on FFT, soft thresholding formula, and ''PCG'' method.
• We conduct many experiments on different blurring kernels, different types of structured noise, and different Gamma noise levels, and illustrate that the proposed model is robust and outperforms the state-of-the-art in image deblurring with multiplicative structured noise removing.
• Meanwhile, numerical results indicate that our algorithm can generate competitive results in terms of the CPU-time comparing to other algorithms and models.