A new approach for Cauchy noise removal

1.   Introduction

In many imaging applications, images inevitably contain noise. Most of the literatures deal with the reconstruction of images corrupted by additive Gaussian noise, for instance ^{[1,2,3,4,5,6]}. However, in many engineering applications the noise has an impulsive characteristic, which is different from Gaussian noise and cannot be modeled by Gaussian noise. Based on ^[7], a type of impulsive degradation is given by Cauchy noise, which follows Cauchy distribution and appears frequently in radar and sonar applications, atmospheric and underwater acoustic images, wireless communication systems. For more details, we refer to ^[8,9]. Recently, much attention has been paid to dealing with Cauchy noise and several approaches have been proposed. In ^[10], Chang et al. employed recursive Markov random field models to reconstruct images corrupted by Cauchy noise. Based on non-Gaussian distributions, Loza et al. ^[11] proposed a statistical approach in the wavelet domain. By combining statistical methods with denoising techniques, Wan et al. ^[12] developed a segmentation approach for RGB images corrupted by Cauchy noise. Sciacchitano et al. ^[13] proposed a total variation (TV)-based variational method for reconstructing images corrupted by Cauchy noise. The variational model in ^[13] (called as SDZ model) is

where $\Omega$ is a bounded connected domain in $\mathbb{R}^2$, $\text{ BV}(\Omega)$ is the space of functions of bounded variation, $u\in \text{ BV}(\Omega)$ (for more details, see (2.1)) represents the restored image and $\gamma > 0$ is the scale parameter of Cauchy distribution. In (1.1), $\lambda$ is a positive number, which controls the trade-off between the TV regularization term and the fidelity term, $u_0$ is the image obtained by applying the median filter ^[14] to the noisy image $f$, and $\eta > 0$ is a penalty parameter. If $8\eta\gamma^2\geq1$, the objective functional in (1.1) is strictly convex and its solution is unique. The term $\eta||u-u_0||_2^2$ in (1.1) results in the solution being close to the median filter result, but the median filter does not always perform well as to Cauchy noise removal. In order to avoid this, in ^[15], the authors developed the the alternating direction method of multipliers (ADMM) to solve the following non-convex variational model (called as MDH model) directly

where $K$ represents a linear operator. As we know, solutions of variational problems with $\text{ TV}$ regularization have many desirable properties, such as the feature of preserving sharp edges. However, these solutions are always accompanied by blocking artifacts due to the property of $\text{ BV}$ space.

In order to overcome blocking artifacts, we will employ TGV as a regularization term. In ^[16], Bredies et al. proposed the concept of TGV, and they applied TGV to mathematical imaging problems to overcome blocking artifacts. For more details of TGV, we refer interested readers to ^[17,18]. In order to overcome the defect of the median filter result, based on the proximal algorithm idea, we will use the term $||u-z||_2^2$ to convexify the non-convex fidelity term $\int_\Omega log(\gamma^2+(u-f)^2)dx$. To simplify computing, for the TGV regularization term, we employ the proximal method. Based on these, we propose the following model

with the constraint $u = z$. Meanwhile, we compare the proposed model (1.3) with the following model

where $u_0$ is the image obtained by applying the median filter ^[14] to the noisy image. According to Table 1, the numerical results show that the proposed model (1.3) is better than the model (1.4). Compared with previous reports, the main novelty of our proposed approach has been condensed into the following points:

1. Compared with the BV regularization term, we employ the TGV regularization term which preserves the image structure and we prove that the proposed model admits a unique solution.

2. Different from the constraint by applying the median filter, we use the constraint by applying the proximal approach and experiment results show better performance.

3. The previous literature used ADMM algorithm but we employ non-expansive operator and the fixed point algorithm such that the convergence of the proposed algorithm is more efficiently proved.

The next part is organized as below. We propose a new model and show the model has a unique solution in Section 2. In Section 3, we employ a minimization scheme to deal with the new model. We show the convergence of the proposed algorithm in Section 4. The performance of the new method is demonstrated by numerical results in Section 5. Some remarks are concluded in Section 6.

2.   Variational model

Similar to ^[13], we propose a new non-convex TGV model for denoising Cauchy noise. For completeness, firstly, a review of BV space and TGV space is given. For more details on TGV models and Cauchy noise removal, we refer to ^{[19,20,21,22]}.

2.1.   Preliminaries

For convenience, we introduce the following notations. The function $u\in \text{ BV}(\Omega)$ iff $u \in {\text{L}^1}(\Omega)$ and its TV is finite, where TV of $u$ is

The space $\text{ BV}(\Omega)$ is a Banach space with the norm $||u||_{\text{ BV}(\Omega)} = ||u||_{{ \text{L}^1}(\Omega)}+\int_\Omega|Du|$ ^[23,24].

Throughout the paper, we denote the dimension by $d$, which is typically 2 or 3. For convenience, $C_c^k(\Omega, \rm Sym^k(\mathbb{R}^d))$ expresses the space of compactly supported symmetric tensor field, where $\rm Sym^k(\mathbb{R}^d)$ represents the symmetric tensors space on $\mathbb{R}^d$, which can be written as ^[16]

The TGV of order $k$ with positive weights $\alpha = (\alpha_0, \alpha_1, \cdot\cdot\cdot, \alpha_{k-1})$ is defined as ^[16]

When $k = 1, \alpha = 1$, $\rm Sym^1(\mathbb{R}^d) = \mathbb{R}^d$, $\rm TGV_\alpha^1(u) = TV(u).$ When $k = 2$, $\rm Sym^2(\mathbb{R}^d)$ represents all symmetric $S^{d\times d}$ matrices as follows, for $\xi\in \rm Sym^2(\mathbb{R}^d)$,

for more details, we refer to ^[16]. In the following part, we mainly use second-order TGV as

where

and

Following the notation in ^[25], we define the discretized grid as

for some positive $N_1, N_2\in \mathbb{N}$, where $h$ denotes the grid width and we take $h = 1$ for convenience. For convenience, we define $U, W, Z$ as

For simplicity, the $\rm TGV_\alpha^2$ functional will be discretized by finite differences with step-size 1. Based on ^[16], $\rm TGV_\alpha^2(u)$ can be reformulated as

where $w = (w_1, w_2)^T \in W$, $\varepsilon(w) = \frac{1}{2}(\nabla w+ \nabla^T w)$. The operators $\nabla$ and $\varepsilon$, respectly, denote

$\nabla: \ \ U\rightarrow W, \ \ \nabla u = \left(\begin{array}{c} \partial^+_x u \\ \partial^+_yu \end{array} \right)$,

$\varepsilon: \ \ W\rightarrow Z, \ \ \varepsilon(w) = \left(\begin{array}{cc} \partial^+_x w_{1} & \frac{1}{2}(\partial_y^+ w_{1}+\partial_x^+w_{2}) \\ \frac{1}{2}(\partial_y^+ w_{1}+\partial_x^+w_{2}) & \partial^+_y w_{2} \end{array} \right)$,

${\text{ div}:}\; W\rightarrow U, \; {\text{ div}}w = \partial_x^-w_1+\partial_y^-w_2$,

${\text{ div}}_h: \ \ Z\rightarrow W, \ \ {\text{ div}}_hz = \left(\begin{array}{cc} \partial^-_x z_{11}+\partial^-_y z_{12}\\ \partial^-_xz_{21}+ \partial^-_y z_{22} \end{array} \right)$.

For more details on the above discretion, we refer to ^[26].

By ^[27], the Cauchy distribution can be written as

where $x$ represents a random variable which obeys the Cauchy distribution, $\mu$ represents the peak location, $\gamma > 0$ is a scale parameter. The scale parameter is similar to the role of the variance. Here, we denote Cauchy distribution by $\mathcal{C}(\mu, \gamma)$.

2.2.   Proposed model

Following ^[13], we denote random variables by $\bf{f, u, v}$ and the respective instances by $f, u, v$. Denote the noisy image by $f = u+v$, where $v$ follows the Cauchy noise. We assume that the $v$ follows the Cauchy distribution with $\mu = 0$ and its density function is defined as follows

The MAP estimator of $u$ is obtained by maximizing the conditional probability of $\bf{u}$ with given $\bf{f}$. Based on Bayes' rule, we have

Equation (2.8) is equivalent to

where the term $logP(f(x)|u(x))$ presents the degradation process between $f$ and $u$, and $logP(u)$ denotes the prior information on $u$. For the Cauchy distribution $\mathcal{C}(0, \gamma)$ and each $x \in \Omega$, we have

In order to overcome blocking artifacts, we employ the prior $P(u) = exp(-\frac{2}{\lambda}{\rm TGV}_{\alpha}^2(u))$. Then we obtain the TGV model for denoising as

where $\lambda > 0$ is the regularization parameter.

Next, we show that problem (2.11) admits at least one solution.

Theorem 2.1. The problem (2.11) has at least one solution in ${\text{ BGV}}_{\alpha}^2(\Omega)$, if $\gamma\geq 1, \lambda > 0$.

Proof. Clearly, if $\gamma\geq 1$, there exists a lower bound of the model (2.11). Assume that $\{ u^k \}_{k\in \mathbb{N}}$ is a minimizing sequence for problem (2.11).

By contradiction, we show that $\{ u^k \}$ is bounded in ${ \text{L}^2}(\Omega)$ and therefore bounded in ${\text{L}^1}(\Omega)$. Assume that $||u^k||_2 = +\infty$, so there exists a set $E \subset \Omega$ and measure($E$)$\neq 0$, such that for any $x \in E$, $u^k(x) = +\infty$. With $f\in {\text{L}^2}(\Omega)$, we have $log(\gamma^2+(u^k-f)^2) = +\infty$ for all $x\in E$, which contradicts to $\int_\Omega log(\gamma^2+(u-f)^2)dx < +\infty$.

Noting that $||\nabla u||_1$ and $||\varepsilon(w)||_1$ are both bounded, we obtain that $\{ u^k\}$ is a bounded sequence in ${\text{ BGV}}_{\alpha}^2(\Omega)$. According to Rellich-Kondrachov compactness theorem, there exists a function $u^* \in {\text{L}^1}(\Omega)$ such that $u^k\rightarrow u^*$. Because ${\rm TGV}_{\alpha}^2(u)$ is proper, semi-continuous and convex in ${\text{ BGV}}_{\alpha}^2(\Omega)$ ^[16], we obtain that ${\rm lim inf}_{k\rightarrow +\infty} {\rm TGV}_\alpha^2(u^k) \geq {\rm TGV}_\alpha^2(u^*).$ Meanwhile, according to the Fatou lemma, we can deduce that

which means that $u^*$ is the minimum point of (2.11), i.e., the problem (2.11) has at least one solution in ${\text{ BGV}}_{\alpha}^2(\Omega)$. Noting that the model (2.11) is strictly convex, based on the standard arguments in convex analysis^[28,29], we obtain that the minimum $u^*$ is unique.

3.   The algorithm for solving (2.11)

In order to obtain a convergent algorithm, we employ the alternating minimization algorithm for the variational model (2.11). The model (2.11) can be discretized as

Remark 3.1. The proximal operator ^[30] $prox_f:\mathbb{R}^n\rightarrow \mathbb{R}^n$ of $f$ is defined as

The definition indicates that $prox_f(v)$ is the trade-off between minimizing $f$ and being near to $v$. Based on this idea, we convexify the model (2.11) by adding a proximal term. The advantages are as follows:

(1) A strictly convex model is obtained due to the proximal term.

(2) The result of each iteration is near to the previous one.

In order to simplify the alternating minimization algorithm, we first introduce the following notations and definitions:

Now, we solve z-subproblem by (3.2). Based on ^[31], $\int_\Omega|Du|$ can be represented by

Equation (3.4) will make the calculation of the primal dual method very easy. Note that ${\rm TGV}_{\alpha}^2(u) = \mathop{\rm inf} \limits_{w\in W} \{ \alpha_0||\nabla u-w||_1+ \alpha_1||\varepsilon(w)||_1 \}$ (for more details, refer to ^[2,32]), where $\alpha_0, \alpha_1$ are positive constant parameters. Therefore the min-max problem of (3.2) can be reformulated as

where $p, q$ are the dual variables associated with the sets given by

Similar to ^[25], the solution of the min-max problem (3.5) can be solved as follows

where the projection can be computed as

$\sigma, \tau$ are positive parameters such that $\sigma\tau\leq 1/12$, and $k, l$ represent iteration numbers.

The optimality condition for (3.3) is

Based on the proximal-operator idea, we can take $v = u^k$ such that the result of each iteration is near to the previous one. Multiplying both sides of (3.11) by $\gamma^2+(u-f)^2$, one can obtain that (3.11) is equivalent to

where

In order to solve (3.12), we need the following proposition.

Proposition 3.2. ^[33] A generic cubic equation with real coefficients

has at least one solution among the real numbers. Let

If there exists a unique real solution of (3.17), the discriminant, $\triangle = q^3+r^2$ has to be positive. Furthermore, if $\triangle \geq 0$, the only real root of (3.17) is given by

Since the problem (3.3) is strictly convex with respect to $u$, then there exists a unique real solution for (3.3) and it can be obtained by (3.19). Instead of the method presented above, the $u$ subproblem (3.3) can be solved by the Newton method because the objective function in (3.3) is twice continuously differentiable.

The alternating minimization algorithm for Cauchy noise removal is given in Algorithm 1, where $K$ represents the maximum iteration number.

4.   Convergence

In the following section, we prove the convergence of the proposed Algorithm 1.

Definition 4.1. (^[34]). An operator $Q:\mathbb{R}^N \rightarrow \mathbb{R}^N$ is non-expansive, if for $\forall y_1, y_2 \in \mathbb{R}^N$, there holds $||Q(y_1)-Q(y_2)||_2\leq ||y_1-y_2||_2$.

Clearly, the identity map $I(x) = x$ for all $x$ is non-expansive. One can easily check that the product and the sum of two non-expansive operators are also non-expansive respectively. For any fixed $v\in \mathbb{R}^N$, the maps $Q(y) = y+v$ and $Q(y) = y-v$ are non-expansive.

Definition 4.2. (^[34]). Given a non-expansive operator $P$, $T = (1-\beta)I+\beta P$, for some $\beta \in (0, 1)$, is said to be $\beta$-averaged non-expansive.

Definition 4.3. (^[34]). An operator $G:\mathbb{R}^N \rightarrow \mathbb{R}^N$ is called firmly non-expansive, if for any $x_1, x_2\in \mathbb{R}^N$, there holds

Remark 4.4. An operator $G$ is firmly non-expansive if and only if it is $\frac{1}{2}$-averaged non-expansive.

Lemma 4.5. (^[35]). Let $\varphi$ be convex and lower semi-continuous, and $\beta > 0$. Suppose $\hat{x}$ is defined as follows

Then $S$ is $\frac{1}{2}$-averaged non-expansive.

Since ${\rm TGV}_{\alpha}^2(u)$ is convex and lower semi-continuous, based on Lemma 4.5, it is obvious that $\mathcal{S}(u)$ is $\frac{1}{2}$-averaged non-expansive. Note that

Let $\varphi(u) = \sum_{i = 1}^{N} log(\gamma^2+(u_i-f)^2)+\frac{1}{2}||u-z^{k}||_2^2$, we have

Noting that $\varphi(u)$ is convex and by Lemma 4.5, we have that $\mathcal{L}(z)$ is $\frac{1}{2}$-averaged non-expansive.

Lemma 4.6. (^[36]) Let $P_1$ and $P_2$ be $\beta_1$-averaged and $\beta_2$-averaged non-expansive operators respectively.

By Lemma 4.5, we obtain that $\mathcal{L}\circ \mathcal{S}$ is $\frac{3}{4}$-averaged non-expansive.

Definition 4.7. (^[37]). A function $\phi: R^n\rightarrow R$ is proper over a set $X\subset R^n$ if $\phi(x) < +\infty$ for at least one $x\in X$ and $\phi(x) > -\infty$ for all $x\in X$.

Definition 4.8. (^[37]). A function $\phi:R^n\rightarrow R$ is coercive over a set $X\subset R^n$ if for every sequence $\{x_k\}\in X$ such that $||x_k||\rightarrow \infty$, we have $\mathop{\rm lim}\limits_{k\rightarrow \infty}\phi(x_k) = \infty$.

The following Lemma 4.9 can be shown easily, and we omit its proof here.

Lemma 4.9. The functional $E(z, u)$ in (3.1) is coercive.

Lemma 4.10. (^[28]). Let $\phi : R^N\rightarrow R$ be a closed, proper and coercive function. Then the set of the minimizers of $\phi$ over $R^N$ is nonempty and compact.

Lemma 4.11. The set of the fixed points of $\mathcal{L}\circ \mathcal{S}$ is non-empty.

Proof. By Lemma 4.9, the objective function $E(z, u)$ is coercive. Based on Lemma 4.10, the set of minimizers of $E(z, u)$ is non-empty. Set $(\hat{z}, \hat{u})$ is a minimizer of $E(z, u)$. Therefore we have

It indicates that

Thus we have $\hat{u} = \mathcal{L}\circ\mathcal{S}(\hat{u})$.

According to the Krasnoselskii-Mann (KM) theorem ^[38], noting that $\mathcal{L}\circ\mathcal{S}$ is non-expansive and the set of the fixed points of $\mathcal{L}\circ\mathcal{S}$ is nonempty, one has that the sequence $\{u_i\}$ converges weakly to a fixed point of $\mathcal{L}\circ\mathcal{S}$, for any initial point $u_0$. Since $E(z, u)$ is strictly convex and differentiable with $u$, one has that the minimizer of $E(z, u)$ is unique. Clearly, the fixed points of $\mathcal{L}\circ\mathcal{S}$ are just the minimizers of $E(z, u)$. Thus the sequence $\{u_k \}$ converges to the unique minimizer of $E(z, u)$. Therefore, we have the following theorem.

Theorem 4.12. The sequence $\{u_k \}$ converges to the unique minimizer of $E(z, u)$ as $k\rightarrow \infty$, for any initial point $u_0$.

5.   Numerical simulations

In this section we provide numerical results to show the performance of the proposed method for image restoration problems under Cauchy noise. Here, we compare our method with existing models as follows: ROF model^[1], the median filter^[39], SDZ model^[13] and MDH model^[15]. For ROF model, we use the primal dual method proposed in ^[31]. For SDZ model and MDH model, we use the source codes of ^[13] and ^[15] respectively.

Considering the quality of the restoration results, we measure them by different evaluation metrics: The peak-signal-to-noise ratio (PSNR) value and the structural similarity (SSIM) value, which are defined as

where $u_0$ denotes the original signal with mean $\bar{u}_0$ and $u$ is the denoised signal, $M\times N$ is the image size,

where $\mu_u, \mu_{u_0}, \sigma_u^2, \sigma_{u_0}^2, \sigma_{uu_0}$ denote, respectively, mean, variance, co-variance of the image $u$ and $u_0$, $c_1$ and $c_2$ are small positive constants. To compare with different approaches easily, we use the same stopping criterion for all the algorithms, that is

or the maximum number of iterations.

Combined with relevant reports^{[13,15,16,17,25,26]}, we adjust each parameter one by one. For each image in Figure 1, we try our best to tune the parameters of the compared algorithms to obtain the highest PSNR and SSIM. Based on hundreds of experiments, we observe that $\tau$ is the key parameter to control the restoration quality and convergence speed. For the proposed model, ROF model, MDH model and the median filter, the grey level range is [0,255]. For SDZ model, the grey level range is normalized to [0, 1]. For the proposed model, the range of $\tau$ is [0.3, 0.7], $\sigma = \tau/12$, the range of $\lambda$ is ^[15,30], and the range of $\eta$ is [0.9, 3]. For MDH model, the range of $\lambda$ is ^[25,50]. For SDZ model, the range of $\lambda$ is ^[2,9] and the range of $\eta$ is [0.3, 3]. For ROF model, the range of $\lambda$ is ^[1,8].

According to the images in Figures 2–4, images visual quality of our method is better than others. Compared with TV-based methods, block-effects can be more significantly reduced by our method. The reason is that the solution of kernel space of second-order $TGV$ is first-order polynomial, but not the piecewise constant function in $BV$ space. In Tables 2 and 3, we list PSNR, SSIM values of numerical results. Clearly, PSNR, SSIM values of our method are better than others. According to the zoomed images in Figures 2–4, our method enhances the image quality and reduces noise more significantly while there is much more noise residual in images of other methods. According to the images in Figure 4, the structure around the eye is better preserved in the proposed method than others.

6.   Conclusions

Based on the Moreau envelop^[30] idea and TGV regularization, we propose a new approach to Cauchy noise removal. We show that the solution of the proposed model is unique. In order to solve the new model, an alternating minimization method is employed and its convergence is proved. Numerical results demonstrate that the images quality of our method is better than that of some earlier restoration methods.

Acknowledgments

We really appreciate authors of ^[13] and ^[15] for their codes.

Conflict of interest

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.