Keywords

1 Introduction

The main goal of an image denoising algorithm is to gain access though both noise suppression and the maintenance of important features such as contours and textures. In many real world applications, images are not usually contaminated with additive noise. This can be expressed in a word the presence of a noise called multiplicative, which appears in various image processing applications, for instance in Synthetic Aperture Radar (SAR) and Ultrasound imaging [10, 15]. For a mathematical description of such a problem we suppose that \(\varOmega \) is a connected open subset of . Let be the original image, from the degraded mechanism of multiplicative noise, the observation image f is obtained by

$$\begin{aligned} f = u\otimes \eta , \end{aligned}$$
(1)

where \(\eta \) is multiplicative noise. In this work, we are interested on the multiplicative noise that follows the Gamma distribution. So far, a many variational models have been considered to manipulate the image restoration problem with the multiplicative noise, while a large amount of works on this subject are total variation (TV)-based. The first approach with the TV regularization devoted to multiplicative noise removal (suited for Gaussian multiplicative noise) also have proposed by Rudin et al. [11] (RLO-model). By using property of the Maximum a Posteriori (MAP) estimator and TV as edge preserver, Aubert and Aujol [2] have introduced a multiplicative denoising variational model (AA-model) based on the Bayes rule and Gamma distribution, whose energy functional can be written as follows

$$\begin{aligned} \min _{u} \left\{ \int _{\varOmega }|\nabla u|+\lambda \int _{\varOmega }(\log u + \frac{f}{u}) \right\} , \end{aligned}$$
(2)

where the first term is the regularization term which imposes some prior constraints on the original image, and the second term is the fitting term which measures the violation of relation between u and the observation image f. Obviously, because of the non-convexity of (2), it does not have the global optimal solution. To make this problem treatable, many approaches [3, 4, 9, 13, 14] have been proposed in the literature. Among those approaches, Huang et al. [9] considered an exponential-transformation \(u\longrightarrow \exp (u)\) in the fitting term of (2) to obtain the following strictly convex model

$$\begin{aligned} \min _{u} \left\{ \lambda _{1}\int _{\varOmega }|\nabla u|+\lambda _{2}\int _{\varOmega }|u-z|^{2}+\int _{\varOmega }(z+f\exp (-z))\right\} , \end{aligned}$$
(3)

where \(\lambda _{1}\) and \(\lambda _{2}\) are positive regularization parameters, and \(z=\log u\) is an auxiliary variable. In a recent work, Huang et al. [8] proposed a significant modification of the regular term in (2). In deed, they provided a non-convex Bayesian type model for multiplicative noise removal, which includes TV and the Weberized TV as regularizers [12]. This gives an effective improvement over the previous models in terms of visual appearance. This model can be described as follows

$$\begin{aligned} \min _{u} \left\{ \alpha _{1}\int _{\varOmega }|\nabla u|+\alpha _{2}\int _{\varOmega }\frac{|\nabla u|}{u}+\int _{\varOmega }(\log u+\frac{f}{u})\right\} , \end{aligned}$$
(4)

where the first two terms are the regularization terms, while the third one is the nonconvex data fitting term of (2). \(\alpha _{1},\alpha _{2}\) are regularization parameters. The main idea of the second regular term in (4) is to discover the universal influence of ambient intensity level u on human’s sensitivity to the local intensity increment \(|\nabla u|\), or so called Just Noticeable Difference (JND) [16], in the perception of both sound and light. This local fluctuation can be formulated by

$$\begin{aligned} \frac{|\nabla u|}{u} = \text{ Const }, \end{aligned}$$
(5)

according to Weber’s law [12, 16], when the mean intensity of the background is increasing with a higher value, then the intensity increments \(|\nabla u|\) also has higher value. In the other hand, Dong et al. [4] proposed to replace the classical TV-norm in [9, 13] with the discrete nonlocal norm one and they applied the split Bregman iteration to solve this two nonlocal models, by introducing the two constraints \(z=\log u\) and \(z_{0}=\log f\), this two models can be described as follows

$$\begin{aligned} \begin{aligned}&\min _{z}~~\left\{ \sum _{i}|\nabla _{NL}z|_{i}+\frac{\lambda }{2}\sum _{i}|z_{i}-z_{0_{i}}|^{2}\right\} ,\\&\min _{z}~~\left\{ \sum _{i}|\nabla _{NL}z|_{i}+\lambda \sum _{i}(z_{i}-\exp (z_{0_{i}}-z_{i}))\}\right\} . \end{aligned} \end{aligned}$$
(6)

Inspired by the idea of [8], we propose a nonlocal variational model for multiplicative noise removal

$$\begin{aligned} \min _{u} \left\{ \alpha _{1}\int _{\varOmega }|\nabla _{NL} u|+\alpha _{2}\int _{\varOmega }\frac{|\nabla _{NL} u|}{u}+\int _{\varOmega }(\log u+\frac{f}{u})\right\} . \end{aligned}$$
(7)

The rest of this paper proceeds as follows. In Sect. 2, first we describe the necessary definitions and tools of some nonlocal operators. After, we give a detailed implementation of the efficient computational method for obtaining the numerical solution of (7) using an improved split Bregman algorithm. Numerical experiments intended for the effectiveness of the proposed method are provided in Sect. 3. Finally, conclusions are made in Sect. 4.

2 Our Proposed Model

2.1 Nonlocal Differential Operators

First, we give a brief overview of the definitions of some nonlocal functionalities (for more details see [1, 5, 6]). Let , be a real function and \(\mu \) be a probability measure. Let a nonnegative symmetric function represents the weights between the current pixel and the pixels in \(\varOmega \). Given a pair of points \((x,y)\in \varOmega \times \varOmega \), then the nonlocal gradient is defined as

$$\begin{aligned} \nabla _{NL}u(x,y) := \sqrt{\mathcal {K}(x,y)}(u(y)-u(x)). \end{aligned}$$
(8)

The inner product between two vectors and the associated norm of \(q_{1}\) at \(x\in \varOmega \) are defined as

$$\begin{aligned} <q_{1},q_{2}> := \int _{\varOmega }q_{1}(x,y)q_{2}(x,y)d\mu (y),\quad |q_{1}| = \sqrt{\int _{\varOmega }q_{1}(x,y)^{2}d\mu (y)}. \end{aligned}$$
(9)

Hence, the norm of \(\nabla _{NL}u\) at \(x\in \varOmega \) is referred to as

$$\begin{aligned} |\nabla _{NL}u|(x) = \sqrt{\int _{\varOmega }\mathcal {K}(x,y)(u(y)-u(x))^{2}d\mu (y)}. \end{aligned}$$
(10)

The nonlocal divergence operator can be defined by the standard adjoint relation with the nonlocal gradient, , as follows

$$\begin{aligned}<\nabla _{NL}u,q> := -<u,\text{ div }_{NL}q>, \end{aligned}$$
(11)

which leads to the definition of nonlocal divergence of the vector q

$$\begin{aligned} \text{ div }_{NL}q(x) := \int _{\varOmega }\sqrt{\mathcal {K}(x,y)}(q(x,y)-q(y,x))d\mu (y). \end{aligned}$$
(12)

The general nonlocal p-Laplacian operator for \(1\le p <\infty \), which is defined by

$$\begin{aligned} \begin{aligned} \varDelta _{NL}^{p}u(x)&:=\frac{1}{2}\text{ div }_{NL}(|\nabla _{NL}u|^{p-2}\nabla _{NL}u),\\ {}&=\int _{\varOmega }\mathcal {K}(x,y)^{\frac{p}{2}}|u(y)-u(x)|^{p-2}(u(y)-u(x)) d\mu (y). \end{aligned} \end{aligned}$$
(13)

2.2 The Proposed Model

In this subsection, we describe a method to solve our model (7). We first take \(\psi (u)=\alpha _{1}+\alpha _{2}/u\), the minimization problem (7) can be rewritten as

$$\begin{aligned} \min _{u}\left\{ \int _{\varOmega }\psi (u)|\nabla _{NL} u|d\mu (x)+\int _{\varOmega }(\log (u)+\frac{f}{u})d\mu (x)\right\} . \end{aligned}$$
(14)

Generally, the nonlocal TV norm in (14) is difficult to be computed straightforwardly. To avoid this drawback, we adopt the split Bregman iteration, which was initially introduced by Goldstein and Osher [7]. This method is originally designed for the \(L_{1}\) regularization \(+ L_{2}\) minimization problem. Note that the second term contains a log-term, which will not allow to apply directly this method. To overtake this difficulty, we proceed as in [4] and introduce an auxiliary variable z, such that \(z=u\) and go on to solve the following unconstrained optimization problem

$$\begin{aligned} \min _{u,z}\left\{ \int _{\varOmega }\psi (u)|\nabla _{NL}u|d\mu (x)+\frac{\lambda }{2}\int _{\varOmega }|u-z|^{2}d\mu (x)+ \int _{\varOmega }(\log z+\frac{f}{z})d\mu (x)\right\} , \end{aligned}$$
(15)

where \(\lambda >0\) should be large enough so that z is sufficiently close to u in the sense of \(L_{2}\)-norm. Subsequently, the minimization problem (14) can be iteratively solved by minimizing the following two subproblems

$$\begin{aligned} \min _{u}= & {} \int _{\varOmega }\psi (u)|\nabla _{NL}u|d\mu (x)+\frac{\lambda }{2}\int _{\varOmega }|u-z|^{2}d\mu (x), \end{aligned}$$
(16)
$$\begin{aligned} \min _{z}= & {} \frac{\lambda }{2}\int _{\varOmega }|u-z|^{2}d\mu (x)+\int _{\varOmega }(\log z+\frac{f}{z})d\mu (x). \end{aligned}$$
(17)

The first subproblem (16) is a \(L_{1}\) regularization \(+L_{2}\) minimization problem, thus it can be efficiently solved by the classical split Bregman iteration. Let \(d=\nabla _{NL}u\), then (16) becomes

$$\begin{aligned} \begin{aligned} (u^{k+1},d^{k+1})=&\min _{u,d}\int _{\varOmega }\psi (u)|d|d\mu (x)+\frac{\beta }{2}\int _{\varOmega }|d-\nabla _{NL}u-b^{k}|^{2}d\mu (x)\\&\qquad \qquad +\frac{\lambda }{2}\int _{\varOmega }|u-z|^{2}d\mu (x),\\ b^{k+1}=&\ b^{k}+\nabla _{NL}u^{k+1}-d^{k+1}, \end{aligned} \end{aligned}$$

where \(\beta \) is a positive regularization parameter. To solve the problem (15) is equivalent to solve the following system

$$\begin{aligned} \begin{aligned} u^{k+1}&=\min _{u}\int _{\varOmega }\psi (u)|d^{k}|d\mu (x)+\frac{\beta }{2}\int _{\varOmega }|d^{k}-\nabla _{NL}u-b^{k}|^{2}+\frac{\lambda }{2}\int _{\varOmega }|u-z^{k}|^{2}d\mu (x),\\ d^{k+1}&=\min _{d}\int _{\varOmega }\psi (u^{k+1})|d|d\mu (x)+\frac{\beta }{2}\int _{\varOmega }|d-\nabla _{NL}u^{k+1}-b^{k}|^{2}d\mu (x),\\ b^{k+1}&=b^{k}+\nabla _{NL}u^{k+1}-d^{k+1},\\ z^{k+1}&=\min _{z}\frac{\lambda }{2}\int _{\varOmega }|u^{k+1}-z|^{2}d\mu (x)+\int _{\varOmega }(\log z+\frac{f}{z})d\mu (x). \end{aligned} \end{aligned}$$
(18)

Taking in the system (18), then (18) can be written as

(19)

The minimizers of first and last subproblems from system (18) are characterized by the optimality condition given by the Euler-Lagrange formulation. The second subproblem from system (18) can be implemented via the generalized shrinkage formula

$$\begin{aligned} shrink(t,\gamma )=sgn(t)\max \left\{ |t|-\gamma ,0\right\} . \end{aligned}$$
(20)

Which are outlined in the following system

$$\begin{aligned} \begin{aligned}&0 = \psi ^{'}(u)|d^{k}|+\lambda (u-z^{k})-\beta \text{ div }_{NL}(d^{k} -\nabla _{NL}u-b^{k}),\\&d^{k+1} = \frac{\nabla _{NL}u^{k+1}+b^{k}}{|\nabla _{NL}u^{k+1}+b^{k}|}\max \{|\nabla _{NL}u^{k+1}+b^{k}|-\frac{\psi (u^{k+1})}{\beta },0\},\\&0 = -\lambda (u^{k+1}-z)+\frac{z-f}{z^{2}}. \end{aligned} \end{aligned}$$
(21)

Using the Gauss-Seidel iterative scheme, then the numerical algorithm summarizes all of these elements

figure a

Let \(\psi ^{'}(u_{j}^{k+1,n}) = -\frac{\alpha _{2}}{(u_{j}^{k+1,n})^{2}}\), \(|d^{k+1,n}| = \sqrt{\sum _{j}(d_{i,j}^{k+1,n})^{2}+(d_{j,i}^{k+1,n})^{2}}\) and

\(d^{k+1,n=0}=d^{k}\), \(b^{k+1,n=0}=b^{k}\) and \(z^{k+1,l=0}=z^{k}\). We choose the weight function as

$$\begin{aligned} \mathcal {K}(x,y) = \exp \left\{ -\frac{d(x,y)}{h_{0}^{2}} \right\} , \end{aligned}$$
(22)

where \(d(x,y)=\int _{\varOmega }G_{\varrho }(s)|u(x+s)-u(y+s)|^{2}ds\) is the distance between patches located at x and y, \(G_{\varrho }\) is a Gaussian kernel of standard deviation \(\varrho \) and \(h_{0}\) is a filtering parameter. In all experiments, we fixed the number of outer iterations in the Split Bregman to be four (i.e., \(k=4\)) and number of inner iterations to be two (i.e., \(n=2\)). We use patches of \(5\times 5\) to compute the weight and a search window of size \(11\times 11\).

Fig. 1.
figure 1

(a) The original Barbara image (\(256\times 256\)); (b) the noisy image with level \(\sigma =0.08\); (c) the result of the AA-model; (d) the result of the HNW-model; (e) the result of the HXW-local-model [\(\alpha _{1}=50, \alpha _{2}=10^{-3}\)]; (f) the result of our model [\(\alpha _{1}=10^{-3}, \alpha _{2}=10^{-4}\)].

Fig. 2.
figure 2

(a) The original Cameraman image (\(256\times 256\)); (b) the noisy image with level \(\sigma =0.1\); (c) the result of the AA-model; (d) the result of the HNW-model; (e) the result of the HXW-local-model [\(\alpha _{1}=10^{-3}, \alpha _{2}=10^{-4}\)]; (f) the result of our model [\(\alpha _{1}=10^{-3}, \alpha _{2}=10^{-4}\)].

Fig. 3.
figure 3

(a) The original Lena image (\(256\times 256\)); (b) the noisy image with level \(\sigma =0.2\); (c) the result of the AA-model; (d) the result of the HNW-model; (e) the result of the HXW-local-model [\(\alpha _{1}=10^{-3}, \alpha _{2}=10^{-4}\)]; (f) the result of our model [\(\alpha _{1}=10^{-3}, \alpha _{2}=10^{-4}\)].

3 Numerical Results

In this section, we present numerical results so as to illustrate the performance of our proposed model. These results are compared to those obtained by the AA-model proposed by Aubert et al. (2), the HNW-model proposed by Huang et al. (3) and the HXW-model proposed by Huang et al. (4). Furthermore, we compared our model with the second nonlocal DZK-model proposed by Dong et al. (6). It’s worth anticipating that all the numerical simulations (exclude DZK-model) are implemented by the approach described in the previous Sect. 2.2. In all tests, each pixel of an original images (Figs. 1(a)–3(a)) are degraded by a noise which follows a Gamma distribution. The noise level is controlled by the value of \(\sigma \) in the experiments, the noisy images with different levels (\(\sigma =0.08,0.1,0.2\)) are shown in Figs. 1(b)–3(b). For measuring the quality of the restored image, three tools are considered. The first one is the structural similarity index (SSIM), the second one is the peak signal to noise ratio (PSNR), and the last one is the signal to noise ratio (SNR), thus these measures are given by

$$\begin{aligned} \begin{array}{c} \text{ SSIM } := \frac{(2\bar{u}\bar{u^{*}}+C_{1})(2\sigma _{uu^{*}}+C_{2})}{(\bar{u}^{2}+\bar{u^{*}}^{2}+C_{1})(\sigma _{u}^{2}+\sigma _{u^{*}}^{2}+C_{2})},\\ \text{ PSNR } := 10\log _{10}\left\{ \frac{M\times N}{||u-u^{*}||_{2}^{2}}\times \max \{u\}^{2}\right\} ,\, \text{ SNR } := 10\log _{10}\frac{||u-\bar{u}||_{2}^{2}}{||u-u^{*}||_{2}^{2}}, \end{array} \end{aligned}$$

where \(u^{*}\), u, and \(M\times N\) are respectively the restored image, the true image and the size of image. \(\bar{u}\) and \(\bar{u}^{*}\) represent the mean values of u and \(u^{*}\) in image domain \(\varOmega \). \(\sigma _{u}\) and \(\sigma _{u^{*}}\) denote the variances of u and \(u^{*}\). \(\sigma _{uu^{*}}\) is the covariance of u and \(u^{*}\). \(C_{1}\) and \(C_{2}\) are two variables to stabilize the division with a low denominator. Figures 1(c)–3(f) show the denoising results of the three noisy gray level images by different methods. In these experiments, it is clear that the restoration results obtained by the proposed method are visually better than those by the AA, HNW and HXW-models. Table 1 summarize the PSNRs and SSIMs results by different methods on the three noisy gray level images. From Table 1, we can see that PSNRs and SSIMs values of restored images using our method are wider than those restored by using the other three methods. In Table 2, we present the PSNRs and SSIMs results obtained by some specific cases of general model (7) on the three noisy gray level images. This can be included the models:

  • When \(\alpha _{1}=0\), this reduces to the nonlocal version of the model in [12] with the fitting term of (2).

  • When \(\alpha _{2}=0\), this reduces to the nonlocal version of the AA-model (2). We choose these tests to show that each model could be the best one among the three models. All the three nonlocal models obtain quite good results, than the other local models. From Table 2, It turns out that the results of the images restored, when noise level is small (\(\sigma =0.08,0.1\)) are better than those of other cases. But, when we increase \(\sigma \) (\(\sigma =0.2\)), we obtained almost the same results. Table 3 shows the SNRs results of our model compared with the nonlocal DZK-model on the two noisy gray level Barbara and Lena images under noise level \(\sigma =0.1\). From Table 3, we can see that our method gives better results than those obtained by DZK-model (6). The main reason for getting better results is that the nonlocal combination of TV and Weberized TV regularizers are dominated in our model. It can successfully remove the multiplicative noise following Gamma distribution, while preserving fine structures and textures as well.

Table 1. The PSNR and SSIM of noisy and restored images using fours methods
Table 2. Comparison results the different cases of general model (7) on the noisy gray level Figs. 1 and 3
Table 3. The SNR comparison of noisy and restored images using the nonlocal DZK-model with our model (\(\alpha _{1}=10^{-3}\), \( \alpha _{2}=10^{4}\)) under noise level \(\sigma =0.1\)

4 Conclusion

In this work, we have proposed a general nonlocal model for multiplicative noise removal problem that following Gamma distribution under the combination of two nonlocal regularizers the weberized TV and the TV operator. By combining these two regularizers, our nonlocal model outperform the classical AA-, HNW-, HXW- and DZK-models. It can both preserve more fine structures and remove the noise of this type. We have tested all of those models under different noise levels and their performances are evaluated and compared both visually and quantitatively.