A variational image denoising model under mixed Cauchy and Gaussian noise

1.   Introduction

Different image acquisition and transmission factors cause the observed images corrupted by a mixture of noise statistics. Image denoising aims to retrieve a clean image from the observed noisy image, which is an essential problem in image processing. We here focus on the image denoising problem in the presence of mixed Cauchy and additive white Gaussian noise. In real circumstances, there are noises with a strong impulsive nature which the Gaussian model fails to describe, such as atmospheric noise caused by lighting, picture noise, radar noise and so on. The Cauchy noise is a type of alpha-stable noise that is impulsive in nature and has a heavy tail. It can occur in low-frequency atmospheric signals ^[1], underwater acoustic signals ^[2,3], radar clutter ^[4,5], multiple access interference in wireless communication systems ^[6], powerline communication channels ^[7], air turbulence ^[8], biomedical images, and synthetic aperture radar images ^[9,10]. Thus, the removal of Cauchy noise has been of great importance in many applications such as sonar, radar, image processing and communications ^[11,12,13]. On the other hand, the additive white Gaussian noise frequently appears due to the temperature of the sensor and the level of illumination in the environment that corrupts every pixels. These noises can appear simultaneously in practice. Indeed, the mixed Cauchy and Gaussian noise can occur in real world applications such as communication systems, where in the receiver one has the sum of the Gaussian noise due to electronic components and impulsive noise due to environmental effects ^[14], or astrophysical image processing ^[15], where the cosmic microwave background radiation is described as the sum of Gaussian distributed from the antenna beam and symmetric alpha-stable distributed random variables from galaxies. Thus, it is a necessary problem to suppress this mixed noise.

Let $\Omega\subset\mathbb{R}^2$ be a connected bounded image domain with a Lipschitz boundary. We consider a noisy image $f:\Omega\rightarrow \mathbb{R}$ given by

where $n$ is the mixture of Cauchy noise and Gaussian noise. The Gaussian noise is assumed to follow a Gaussian distribution, $\mathcal{N}(0, \sigma^2)$, with zero mean and standard deviation $\sigma$, and the Cauchy noise is assumed to follow a Cauchy distribution, $\mathcal{C}(0, \gamma)$. Specifically, the Cauchy noise, $w$, is a random variable following a Cauchy distribution, denoted by $\mathcal{C}(\delta, \gamma)$, with the probability density function (PDF) ^[16,17]

where $\delta\in\mathbb{R}$ is the parameter representing the location of the peak, and $\gamma > 0$ is the scale parameter that determines the level of noise. The Cauchy distribution looks similar to a Gaussian distribution with a bell-shaped curve, but it is a heavy-tailed distribution, as shown in Figure 1(a). Its tail heaviness is determined by the parameter $\gamma$, and it increases as the value of $\gamma$ increases, which can be seen in Figure 1(b). Owing to the tail heaviness, the Cauchy distribution is more prone to producing values that fall far from its mean. Hence, the noise generated from the Cauchy distribution is more impulsive than the Gaussian one. For instance, the Cauchy noise tends to contain much larger noise spikes than the Gaussian noise. In this work, we intend to recover a clean image $u$ from the noisy image $f$, which is an ill-posed inverse problem.

To solve the ill-posed inverse problem (1.1), one of the well-known approaches is to solve a minimization problem of the following form:

where $\Phi$ is a data-fidelity term that depends on the type of noise, and $R$ is a regularization term that controls the smoothness of $u$. The most popular regularization is the total variation (TV) ^[18], due to its convexity and edge preserving property. Many TV based variational models have been proposed for restoring images with various types of noise, such as Gaussian noise ^[18], impulse noise ^[19,20,21], multiplicative noise ^{[22,23,24,25]}, Poisson noise ^[26], Rician noise ^[27,28], and Cauchy noise ^[29,30,31]. The different noise distributions yield different data-fidelity terms, so one can attain suitable variational image denoising models depending on noise types. For instance, the $L^2$ data-fidelity term, $\Phi(u, f) = \int_{\Omega} (f-u)^2\, dx$, is typically used for the image denoising models under Gaussian noise ^[18], and $\Phi(u, f) = \int_{\Omega} |f-u|\, dx$ is more appropriate for the ones under impluse noise ^[19].

While there are numerous denoising methods for Gaussian noise, several approaches to eliminate Cauchy noise have been suggested. In addition to the Markov random field or wavelet-based denoising methods ^[32,33,34], a TV-based model was proposed in ^[29], with the nonconvex data-fidelity term derived from the Cauchy distribution (1.2). The same authors also suggested a convex model by inserting a quadratic penalty term that involves a pre-denoised image achieved by applying the median filtering to the noisy data. However, the median filtering does not always yield satisfactory denoising results. In ^[30], Mei et al. showed the effectiveness of the nonconvex model in ^[29] combined with the nonconvex alternating direction method of multipliers (ADMM) ^[35]. Moreover, a convex TV model was proposed in ^[31] for restoring images with $\alpha$-stable noise with $\alpha\in(0, 2)$ (Cauchy noise when $\alpha = 1$, Gaussian noise when $\alpha = 2$). The experimental results showed that the model outperformed the $L^1$-TV ^[36] and Cauchy ^[29] models in impulsive noisy environments (i.e., $\alpha\in(0, 1.5)$), while providing comparable performance in less impulsive noisy environments (i.e., $\alpha\in (1.5, 2)$). The aforementioned nonconvex or convex TV models in ^[29] have also been extended in various works ^{[37,38,39,40,41,42,43]} by adopting other regularization terms instead of TV.

Removing the mixture of noise is more challenging because of the unique nature of each of two types of noise. For the removal of mixed impulse and Gaussian noise, various efficient two-phase methods that integrate variational methods with adaptive median filters have been proposed ^[44,45,46]. These two-phase approaches strongly rely on precise detection of noisy pixels. A unified framework of joint detection of noisy pixels and reduction of noise components has been suggested in ^[47,48,49]. Moreover, the combination of data-fidelity terms has been considered. Specifically, a linear combination of the $L^1$ and $L^2$ data-fidelity terms was considered ^[50,51]. The removal of mixed Poisson and Gaussian (MPG) has also been extensively studied. Some early works are based on the noise parameter estimation ^[52,53], while others are mainly based on transform-domain procedures ^{[54,55,56,57,58]}, such as variance-stabilizing transformation ^[59] or Haar wavelet transform. The MAP approaches ^[60,61] lead to practical difficulties since the log-likelihood function involves infinite summation. The combination of data-fidelity terms has also been considered for MPG removal ^{[61,62,63,64]}. In ^[65,66], the authors proposed new TV-based denoising models under the mixed salt-and-pepper and Gaussian noise or MPG noise, by utilizing a data discrepancy which characterizes the mixed noise as an infimal convolution of two noise distributions. These data-fidelity terms provided better denoising performance than a combination of different data-fidelity terms corresponding to noise types. In ^[67], new operator-splitting algorithms were suggested for solving the MPG model in ^[65]. On the other hand, there are a few works for dealing with mixed Cauchy and Gaussian noise, and they are only for 1D signals ^{[14,15,68,69]}. In this work, we introduce a new model for restoring images under mixed Cauchy and Gaussian noise, by following the idea of the work ^[65,66].

To the best of our knowledge, there is no unified model for simultaneously removing both Cauchy and Gaussian noise from images, so our main contribution is to propose a novel model for denoising images with the mixed noise. In spite of the nonconvexity of the proposed model, the existence of a solution is proved. We also present an efficient iterative optimization algorithm. The rest of the paper is organized as follows. Section 2 recalls some variational models for restoring images with Cauchy noise. In Section 3, we propose a minimization problem for image denoising under mixed noise, along with an optimization algorithm for solving the proposed model. Section 4 presents experimental results with comparisons to other existing models, and Section 6 concludes our work with some comments.

2.   Variational models for Cauchy noise removal

Assuming that the Cauchy noise follows a zero-centered Cauchy law ($\delta = 0$), Sciacchitano et al. ^[29] derived a TV based model for the removal of Cauchy noise:

where $\text{BV}(\Omega)$ is the subspace of functions $u\in L^1(\Omega)$ such that the following BV semi-norm is finite ^[18]:

where the vector measure $Du$ represents the distributional or weak gradient of $u$, and $\|\cdot\|_\infty$ is the essential supremum norm. This is also called the total variation of $u$, denoted by $\text{TV}(u)$. If $u\in w^{1, 1}(\Omega)$, then $\int_\Omega |Du| = \int_\Omega |\nabla u|\, dx$. $\lambda > 0$ is a tuning parameter that determines the smoothness of the restored image $u$. Despite the nonconvexity of the model (2.1), the existence and uniqueness of a minimizer was proved under certian conditions. Mei et al. ^[30] also showed the efficiency of the model (2.1) associated with the alternating direction method of multipliers (ADMM) ^[35].

Recently, Yang et al. ^[31] proposed a convex TV model for restoring images degenerated by the $\alpha$-stable noise ($0 < \alpha < 2$), including Cauchy noise ($\alpha = 1$):

where $\mu > 0$ is a parameter, and $g$ is the pre-denoised image obtained by applying a median filter to $f$. The existence and uniqueness of a minimizer for the model (2.3) was shown, and the model was efficiently solved by employing the primal-dual algorithm ^[70]. The numerical results showed that the model (2.3) provided better denoising performance than the existing models ^[20,29], especially when the noise has more impulsive properties (i.e., $\alpha\in(0, 1.5)$).

3.   Description of the proposed model and algorithm

3.1.   Proposed model and convergence analysis

In this section, we introduce a new image denoising model in the presence of the mixture of Cauchy noise and Gaussian noise.

In general, both Cauchy noise and Gaussian noise are additive noise, so we assume that both noises occur simultaneously and independently in the entire domain. Thus, we consider a noisy image $f\in L^2(\Omega)$ given by

where $w$ is the Cauchy noise following the Cauchy distribution $\mathcal{C}(0, \gamma)$, and $v$ is the Gaussian noise following the Gaussian distribution $\mathcal{N}(0, \sigma^2)$. Since $w$ and $v$ occur independently, the order of $w$ and $v$ does not matter.

To eliminate both Gaussian noise $v$ and Cauchy noise $w = f-u-v$ from the data $f$, we follow the idea in ^[65] which suggested infimal convolution-type data-fidelity terms for the dismissal of mixed SP-Gaussian noise or mixed Poisson-Gaussian noise. Now we define an infimal convolution-type data-fidelity term to remove mixed Cauchy and Gaussian noise as

where $\Phi_1$ and $\Phi_2$ are the Gaussian and Cauchy noise components, respectively, defined as

where $\Phi_2$ is acquired from the data-fidelity term in (2.1), which is nonconvex. $\lambda_1$ and $\lambda_2$ are positive parameters that balances the smoothing effect of the regularization as well as the fitting with respect to the intensity of each single noise distribution in $f$.

By integrating the data-fidelity term in (3.2) with the TV regularization, we propose the following minimization problem for restoring images with mixed Cauchy and Gaussian noise:

where $X = BV(\Omega)\cap L^2(\Omega)$, $g$ is the pre-denoised image by applying a median filter to the data $f$, and $\mu > 0$ is a parameter. The median filtering does not always bring adequate denoising results, but we add the last quadratic term mainly for the proof of the existence of a minimizer. Indeed, we in practice set the value of $\mu$ very small so that the pre-denoised image $g$ barely impacts on the denoising results. We note that the proposed model can be extended, by utilizing other regularization terms, instead of TV, such as the higher-order regularization ^[71,72], nonlocal TV ^[73], dictionary learning ^[74,75], or tight-frame approach ^[76], etc. This work mainly focuses on introducing a new data-fidelity term for the mixed Cauchy-Gaussian noise model, so the adoption of new regularization terms are beyond the scope of our work.

The energy functional $\mathcal{E}$ in (3.4) is nonconvex without any specific assumptions on the parameters, but it is still available to show the existence of a minimizer for the minimization problem (3.4). In the following theorem, we prove the existence of a minimizer for the problem (3.4).

Theorem 1. Let ${f\in L^2(\Omega)}$, and $g\in L^2(\Omega)$. Then, the minimization problem (3.4) has at least one solution $(u^\ast, v^\ast) \in X\times L^2(\Omega)$ with $X = BV(\Omega)\cap L^2(\Omega)$.

Proof. Since $\Phi_2$ has the minimum value $2|\Omega|\log \gamma$ when $u+v = f$, the functional $E$ in (3.4) is bounded from below. Then we can choose a minimizing sequence $(u_n, v_n)$ in $X\times L^2(\Omega)$ such that $\mathcal{E}(u_n, v_n)\leq C$ for a constant $C > 0$. So all terms in $\mathcal{E}(u_n, v_n)$ are bounded, i.e.,

Hence we can extract a subsequence $\{v_{n}\}$ (still denoted in the same way) weakly converging to $v^\ast$ in $L^2(\Omega)$ and $v_n\rightarrow v^\ast$ a.e. in $\Omega$. Then we have that

Since $g\in L^2(\Omega)$, $\{u_n\}$ is bounded in $L^2(\Omega)$ from (3.5) and thus bounded in $L^1(\Omega)$. Besides, $\{ \int_{\Omega} |D u_n|\}$ is bounded, so $\{u_n\}$ is bounded in $BV(\Omega)$. Hence, there is a subsequence $\{u_{n}\}$ (still denoted in the same way) and $u^\ast$ in $BV(\Omega)$ such that (i) $u_{n} \rightarrow u^\ast$ strongly in $L^1(\Omega)$, (ii) $u_{n}\rightharpoonup u^\ast$ weakly in $L^2(\Omega)$, (iii) $u_{n}\rightarrow u^\ast$ in a.e. $\Omega$, and (iv)

Lastly, from Fatou's Lemma, we can finally attain that

which indicates that $(u^\ast, v^\ast)$ is a minimzer of $\mathcal{E}(u, v)$.

Remark. The functional $\mathcal{E}$ in (3.4) is convex under certain conditions. Specifically, it can be easily proven that $\mathcal{E}(u, \cdot)$ is strictly convex if $\lambda_2\leq 4\gamma^2 \mu$, while $\mathcal{E}(\cdot, v)$ is strictly convex if $\lambda_2\leq 8\gamma^2 \lambda_1$. Furthermore, if $\lambda_2 < \frac{1}{2} \min\{8\gamma^2 \lambda_1, 4\gamma^2 \mu\}$, then $\mathcal{E}(u, v)$ is strictly convex (for the proof, see Appendix A), so it has a unique minimizer. However, we in practice do not impose this condition on the parameter $\lambda_2$ due to the following reason: we want to reduce the influence of the pre-denoised image $g$, so we set the value of $\mu$ to be very small. This leads to a small value of $\lambda_2$ with the aforementioned condition, but a small value of $\lambda_2$ is not appropriate for the removal of Cauchy noise especially when the level of Cauchy noise is high. Therefore, without the above constraint on $\lambda_2$, the functional $\mathcal{E}$ is nonconvex, hence we practically solve the nonconvex minimization problem.

3.2.   Optimization algorithm

To solve the proposed model (3.4), we consider a discretized image domain $\Omega = \{(i, j):\, i = 1, 2, ..., M, j = 1, 2, ..., N\}$, and let $u_s$ be the pixel value of an image $u$ at location $s\in\Omega$. Then we compute numerically the solution pair of the following minimization problem:

where $\langle\, , \, \rangle$ is the inner product, $\|\cdot\|_2^2 = \langle \cdot, \cdot\rangle$, and $\|\nabla u\|_1$ is the discrete version of the isotropic TV norm:

with $\nabla u = [\partial_{x_1} u, \partial_{x_2} u]^T$, denoting $\partial_{x_1} u$ and $\partial_{x_2} u$ by the finite difference operators that estimate the partial derivatives of the image $u$ along the $x_1$-axis and $x_2$-axis, respectively.

To solve the nonconvex problem (3.9), we first adopt the alternating minimization algorithm (AMA). The AMA minimizes a function of two variables, and its essential idea is to keep one variable fixed while minimizing the other variable and iterate this process. This approach has practically performed well, even though the objective function is nonconvex. The AMA applied to (3.9) yields the following iterative algorithm:

In the subsequent paragraphs, we solve the two subproblems in (3.10).

3.2.1.   Solving the $u$-subprolem in (3.10)

First, to solve the $u$-subproblem in (3.10), we adopt the alternating direction method of multipliers (ADMM) ^[35] applied to nonconvex minimization problems with linear constraints. By introducing an auxiliary variable $z$, we can rewrite the $u$-subproblem in (3.10) as the following equivalent constrained problem:

Then the augmented Lagrangian function (ALF) corresponding to (3.11) is given by

where $p\in \mathbb{R}^{M\times N}$ is the Lagrangian multiplier, and $\tau > 0$ is a penalty parameter.

The ADMM applied to (3.11) leads to the following iterative algorithm:

Following the Theorem 4.2 in ^[30], we can prove a convergence of the ADMM in (3.13) as follows:

Theorem 2. If $\tau > \frac{2\lambda_2}{\gamma^2}-\mu$, then the sequence $\{(u^\ell, z^\ell, p^\ell)\}$ generated by the algorithm (3.13) converges globally to a point $(u^\ast, z^\ast, p^\ast)$, which is a stationary point of $\mathcal{L}_\tau$.

The proof of this theorem is omitted since it is similar to the proof given in ^[30]. The stationary point $(u^\ast, z^\ast, p^\ast)$ satisfies the Karush–Kuhn–Tucker (KKT) conditions of problem (3.11). However, the minimization problem (3.11) is non-convex, so the KKT conditions are only necessary optimal conditions to (3.11). Hence, we cannot assure that $(u^\ast, z^\ast)$ is an optimal solution of (3.11).

In the following paragraphs, we solve the two subproblems in (3.13).

$u$-subproblem in (3.13). First, the $u$-subproblem in (3.13) is convex, so it can be efficiently solved by various convex optimization algorithms ^{[70,77,78,79,80]}. Among them, we again utilize the ADMM ^[79,80]. Specifically, to handle the nondifferential $L^1$ term, we introduce an auxiliary variable $\vec{d}$ to replace $\nabla u$, and then achieve the following constrained problem:

Analogously, the ALF for the problem (3.14) is

where $\vec{q}\in (\mathbb{R}^{M\times N})^2$ is the vector of Lagrangian multipliers, and $\eta > 0$ is a penalty parameter. Then the ADMM algorithm applied to (3.14) results in the following iterative algorithm:

The $u$-subproblem in (3.16) is a least squares problem, so the solution $u^{m+1}$ in (3.16) can be attained by solving the following normal equation:

where $\nabla^T = -{div}$ with a discrete divergence opterator such that $div(w_1, w_2) = \partial_{x_1}w_1 + \partial_{x_2}w_2$. In discrete setting, $\nabla^T \nabla$ is equal to $-\Delta$, where $\Delta$ is a discrete Laplacian operator. The discrete Laplacian operator can be regarded as the convolution of the kernel $[0\, \, 1\, \, 0;1\, \, -4\, \, 1;0\, \, 1 \, \, 0]$. Then, since $\Delta$ is a symmetric convolution operator, $\mu+\tau + \eta \nabla^T \nabla$ can be diagonalized by the 2-dimensional discrete cosine transform (DCT2) under the symmetric boundary condition ^[81]. Thus we solve the Eq (3.17) using the DCT2, denoted by $\mathcal{F}$, under the symmetric boundary condition. Then we can obtain an explicit formula for $u^{m+1}$:

where $\mathcal{F}^{-1}$ denotes the inverse DCT2.

Lastly, the solution $\vec{d}^{m+1}$ in (3.16) can be explicitly obtained as

where $shrink$ is the soft thresholding operator defined as

where $|t_s| = \sqrt{(t_{1, s})^2+(t_{2, s})^2}$ with $t = [t_1, t_2]^T$.

$z$-subproblem in (3.13). Next we solve the $z$-subproblem in (3.13), which can be rewritten as

The first-order optimality condition for $z^{\ell+1}$ is given by

where $F(z) = \lambda_2\log(\gamma^2+(z+v^k-f)^2)+\frac{\tau}{2} (z- u^{\ell+1}-p^\ell/\tau)^2$. This normal Eq (3.22) can be efficiently solved by using Newton's method as follows:

The $z$-subproblem (3.21) is strictly convex if $\lambda_2 \leq 4\gamma^2\tau$. This condition is satisfied in practice, so Eq (3.22) has a unique real root. Hence, with a good initial guess, Newton's method converges fast within a few number of iterations. Meanwhile, Eq (3.22) can be rewritten as a cubic equation by multiplying with the denominator. This cubic equation can be explicitly solved using Cardano's formula ^[82]. In our simulations, we decide to utilize Newton's method due to its efficiency. Indeed, despite the same denoising results, the total computational time when using Newton's method is shorter than when using Cardano's formula.

3.2.2.   Solving the $v$-subprolem in (3.10)

Lastly, we solve the $v$-subproblem in the AMA (3.10):

The necessary optimality condition for $v^{k+1}$ is

The problem (3.24) is convex if $\lambda_2 \leq 8\gamma^2\lambda_1$. In practice, this condition is not satisfied due to the choice of small values for $\lambda_1$. Thus, the problem (3.24) is not convex, so we utilize the Cardano's formula to solve Eq (3.25). Specifically, Eq (3.25) can be written as the following cubic equation:

with $a = \lambda_1$, $b = 2 \lambda_1(u^{k+1}-f)$, $c = \lambda_1 (\gamma^2+(u^{k+1}-f)^2)+\lambda_2$, $d = \lambda_2(u^{k+1}-f)$. The Cardano's formula for real roots of a cubic polynomial is given in the following proposition ^[82].

Proposition 1. For the cubic equation with real coefficients

define $Q = \frac{3ac-b^2}{9a^2}$, $R = \frac{9abc-27a^2d-2b^3}{54 a^3}$, and the discriminant $D = Q^3+R^2$. If $D > 0$, the cubic equation (3.27) has only one real root, which is given by

Otherwise, if $D \leq 0$, Eq (3.27) has three real roots (possibly equal), which are given by

where $\theta = \cos^{-1}(R/\sqrt{-Q^3})$.

If the Eq (3.25) has three real roots, we choose one which yields the minimum value of the objective function in (3.24).

Consequently, the whole algorithm for solving the problem (3.9) is given in Algorithm 1.

4.   Experimental results

This section presents the experimental results of the proposed model (3.4) and comparisons to other existing models. We compare our model with two Cauchy denoising models, nonconvex TV model (Cauchy-TV) ^[29,30] and Yang et al.'s model ^[31] (Yang's), given in (2.1) and (2.3), respectively. Furthermore, due to the impulsive characteristics of Cauchy noise, we also present the denoising results of the median filtering (MF) ^[83] and $L^1$-TV model ^[19]. For solving the $L^1$-TV and Cauchy-TV models, we employ the convex or nonconvex ADMM algorithms.

4.1.   Experimental setting

The test images are given in Figure 2, and the size of images are $256\times 256$ or $481\times 321$. The range of intensity values in original images is assumed to be $[0, 255]$. The Cauchy noise $w$ is generated by applying the following property: If $X$ and $Y$ are two independent Gaussian distributed random variables with mean $0$ and variance $1$, then the ratio $X/Y$ follows the standard Cauchy distribution, i.e., $\mathcal{C}(0, 1)$ ^[84,85]. Thus $w = \gamma\frac{n_1}{n_2}$, where $n_1$ and $n_2$ follow the standard normal distribution, $\mathcal{N}(0, 1)$, independently. In the experiments, we consider the four mixed noise cases: $(\gamma, \sigma) = (10, 20)$, $(10, 30)$, $(15, 10)$, $(15, 20)$. In the Yang's model (2.3), the parameter $\gamma$ is set to be the exact Cauchy noise level like our model.

The quality of restored images is measured by the Peak-Signal-to-Noise-Ratio (PSNR) value as

where $u$ and $u_\ast$ are the restored and original images respectively. We also compute the structural similarity index measure (SSIM) ^[86], which is a perception-based measure that uses information about the structure of the objects in the visual sense. Specifically, we compute the mean SSIM index value between two images $u$ and $u_\ast$ using

where $\mathbf{x}$ and $\mathbf{y}$ are the spatial patches extracted from $u$ and $u_\ast$ respectively, $\mu_{\mathbf{x}}$ and $\sigma_{\mathbf{x}}^2$ represent the average and variance of $\mathbf{x}$ respectively, and $c_1$ and $c_2$ are some constants for stability. For a noisy data $f$, the intensity values of some noisy pixels are much larger (or smaller) than 255 (or 0) and they mainly affect the PSNR value of $f$, so the PSNR of $f$ cannot be properly computed. Thus, we compute the PSNR value of the cropped image, $\max(\min(f, 255), 0)$, instead of $f$, which are given in Tables 1–4. Moreover, we use this cropped image as an initial condition for $u$ from the experiments in ^[30].

The stopping criterion for our model is given by

where $g = u$ or $z$, $tol$ is a given tolerance, and $MaxIter$ is a maximum iteration number. For our model, we set $tol = 10^{-4}$. For the $L^1$-TV and Yang's models, we use the same stopping condition as ours, while for the Cauchy-TV model we use the stopping condition given in ^[30]. For our model, the maximum iteration number for Newton's method for $z$ is fixed as $N_z = 5$, and we fix $N_u = 10$ and $N_{uu} = 5$.

The selection of parameters for our model is as follows: The parameter $\mu$ is fixed as $10^{-7}$, and the penalty parameters $\tau$ and $\eta$ in the ADMM algorithms are fixed as $1$. The parameters $\lambda_1$ and $\lambda_2$ mainly influences on the quality of the restored images, and they are tuned to achieve the best denoising results. Indeed, the parameter $\lambda_1$ is set to be $0.05$ when $(\gamma, \sigma) = (15, 10)$, and $0.02$ for the other noise cases. The values of the parameter $\lambda_2$ are given in all figures.

4.2.   Denoising results

First, Figure 3 compares noisy images and signals degraded by the Cauchy noise with $\gamma = 10$, Gaussian noise with $\sigma = 20$, mixed Cauchy and Gaussian noise with $(\gamma, \sigma) = (10, 20)$, respectively. The cross-sectional lines in the right column show that the vertical scale for the noisy signal corrupted by the Cauchy noise goes from -200 to 900, which shows the impulsive feature of the Cauchy noise. Besides, although the noisy signal distorted by the mixed noise is mainly influenced by the Cauchy noise, it can be seen that the Gaussian noise also further distorts the signal.

In Figure 4, we compare the denoising results of all models when the noise level is $(\gamma, \sigma) = (10, 20)$. The data images are given in Figure 5. First we can see that the MF enables to remove both noise to some extent, but its restored images are noisier and blurrier than the ones using the other TV models. Among the TV models, the $L^1$-TV provides similar denoising results to the Cauchy-TV model. This is due to that the $L^1$ data-fidelity is suitable for eliminating both impulsive noise and Gaussian noise. On the other hand, the Yang's model supplies slightly better or similar denoising results than the $L^1$-TV model, which also explains the capability of the Yang's model to handle both noise. In fact, the Yang's model provides averagely higher PSNR and SSIM values than the $L^1$-TV model in all noise cases, as shown in Tables 1–4. Compared with these existing models, our model generates cleaner homogeneous regions, such as the sky areas in Boat and Cameraman, while better conserving textures and details, such as the textural part in Barbara, the ropes and iron pillars in Boat, the face and tripod parts in Cameraman. This can be more clearly seen in the zoomed images in Figure 6. Moreover, all these observations also match the highest PSNR and SSIM values of our model, as shown in Table 1. Overall, these examples justify the effectiveness of our data-fidelity term for getting rid of mixed Cauchy and Gaussian noise.

In Figure 7, we present the denoising results with higher Cauchy noise level and smaller Gaussian noise level, i.e., $(\gamma, \sigma) = (15, 10)$. Unlike the denoising results when $(\gamma, \sigma) = (10, 20)$, the Cauchy-TV model furnishes better denoising results than the $L^1$-TV and Yang's model. This implies that as the Cauchy noise level increases while the Gaussian noise level decreases, the performance of the Cauchy-TV model much enhances unlike the other models. In particular, the Yang's model retain some Cauchy noise in the restored images, which can be more obviously seen in the zoomed images in Figure 8. Indeed, this problem might be fixed by decreasing the regularization parameter $\lambda$ or increasing the size of a median filter applied to $g$, but these cause much smoother restored images with much less PSNR values. Thus we choose these given images as the best restored images of the Yang's model. Despite the improved performance of the Cauchy-TV model, our model provides better preserved details and cleaner smooth regions than the other models, which can be seen in the window parts in Castle and sky areas in Policemen. These also lead to the highest PSNR values of our model. These also show the outstanding performance of our model over the other models. Lastly, Figure 9 presents the final Gaussian components, $v$, corresponding to the final restored images, $u$, given in 5. This validates that our model properly separates the Gaussian and Cauchy noise components.

Figure 10 compares the denoising results when $(\gamma, \sigma) = (10, 20)$, $(10, 30)$, $(15, 20)$. First, it can be observed that the Cauchy-TV model does not perform well in the presence of high-level of Gaussian noise, although it properly removes the Cauchy noise. Moreover, the Yang's model seems to deal with both noise to a certain degree, but it keeps some Cauchy noise under a high level of Cauchy noise, as shown in the denoised image when $(\gamma, \sigma) = (15, 20)$. The $L^1$-TV model moderately performs for deleting both noise, but its performance does not exceed the Yang's one with respect to the PSNR and SSIM values. In contrast, our model well eliminates both noise in all noise cases, by sufficiently denoising homogenous regions while conserving fine structures and edges than the other models.

In Figure 11, we present zoomed denoised images when $(\gamma, \sigma) = (10, 30)$ and $(\gamma, \sigma) = (15, 20)$, respectively. The zoomed data images are given in the first row. When $(\gamma, \sigma) = (10, 30)$, the denoised images of our model and Yang's model seem to be visually indistinguishable, but our model brings cleaner homogeneous regions, such as the sky areas in both examples. On the other hand, the Yang's model fails to adequately remove the Cauchy noise, whose trace can be detected in some pixels in the Pirate image. But our model sufficiently denoise smooth regions with better keeping delicate features than the other models, leading to the highest PSNR values. These can be observed in the eye area in Parrot and face region in Pirate. Hence, our model attains more satisfactory denoising results even when the noise level of Cauchy or Gaussian noise is high.

In Tables 1–4, we report the PSNR and SSIM values of the restored images of all models. The proposed model yields the highest PSNR and SSIM values in almost all cases. Overall, our model leads to the best denoising results with regard to these image quality assessments. These also confirm the superior performance of our model over the existing models.

Figure 12 presents the plots of the energy values $\mathcal{E}(u^k, v^k)$ in (3.4) and PSNR values of $u^k$ via the outer iteration numbers $k$. In all cases, as the outer iteration increases, the energy values decrease and the PSNR values increase. Moreover, in Figure 13, we present the plots of relative errors of $u$ and $v$. We can see that these relative errors generally decrease as the iteration increases, despite slight fluctuations in the relative errors of $v$. All of these illustrate the convergence behavior of the proposed iterative algorithm even though it is not theoretically proven.

In Figures 14 and 15, we compare the denoising results of the proposed model with different values for the parameter $\lambda_2$ in (3.2). In Figure 14, we present the denoising results when $(\gamma, \sigma) = (10, 20)$ and $(10, 30)$, while in Figure 15, we present the denoising resuls when $(\gamma, \sigma) = (10, 20)$ and $(15, 20)$. First, it can be observed that as $\lambda_2$ increases, the restored image $u$ includes more details, but it also tends to retain noise. Thus, we try to select a sufficiently denoised image with PSNR and SSIM values as high as possible for the best restored image. The optimal values of $\lambda_2$ may be different depending on images, as shown in the 1st rows in both figures, because different images have different structures and characteristics. Moreover, the optimal value of $\lambda_2$ depends on the noise level, as $\lambda_2$ is a regularization parameter that needs to be adjusted to obtain an adequately denoised image. Specifically, the optimal value of $\lambda_2$ tends to decrease as the Gaussian noise level $\sigma$ increases, while it tends to increase as the Cauchy noise level $\gamma$ increases, which can be seen in the 2nd rows in both figures. However, we note that the optimal values of $\lambda_2$ are restricted to some values; for all the test images except Bird, the optimal $\lambda_2$ is chosen from $\{21, 22, 23, 24\}$ when $(\gamma, \sigma) = (10, 20)$, $\{20, 21, 22, 23\}$ when $(\gamma, \sigma) = (10, 30)$ and $\{29, 30, 31, 32, 33, 34\}$ when $(\gamma, \sigma) = (15, 20)$, while $\lambda_1$ is fixed at $0.02$. Furthermore, the restored images obtained by using the values in the same set do not change significantly and their PSNR and SSIM values change slightly, as shown in these examples. Hence, these examples show that the proposed model is not sensitive to the choice of parameter $\lambda_2$.

Figure 16 presents the denoising results of the proposed model with two different values of $\lambda_1$ when $(\gamma, \sigma) = (15, 10)$. We recall that for the other noise cases such as $(\gamma, \sigma) = (10, 20)$, $(10, 30)$ and $(15, 20)$, $\lambda_1$ is fixed at $0.02$. However, when $(\gamma, \sigma) = (15, 10)$, the restored images obtained with $\lambda_1 = 0.02$ tend to retain some Cauchy noise despite the parameter $\lambda_2$ being adjusted, as seen in the top row. This is due to that the Gaussian noise level is smaller than the other noise cases, which indicates that the parameter $\lambda_1$ is influenced by the Gaussian noise level. By increasing the value of $\lambda_1$ to $0.05$, we could obtain satisfactory restored images, as shown in the bottom row. We note that when $(\gamma, \sigma) = (15, 10)$, the optimal $\lambda_2$ is chosen from $\{27, 28, 29, 30, 31\}$ for all the test images. Thus, although the parameter $\lambda_1$ depends on the noise level, two values of $\lambda_1$ are enough to attain satisfactory denoising results. These imply that the choice of parameter $\lambda_1$ is not tricky in the proposed model.

In Figure 17, we present some denoising results of the proposed model with different values for the parameter $\gamma$ in (3.2). Throughout the experiments, for the Yang's model and our model, the parameter $\gamma$ is assumed to be the exact Cauchy noise level, denoted by $\gamma_\ast$. So we here present the denoising results when $\gamma$ is not the accurate Cauchy noise level. Two values for $\gamma$ are roughly chosen; one is smaller than $\gamma_\ast$ while the other one is larger than $\gamma_\ast$. For the Lena image, we test with $\gamma = 5$, $15$ when $\gamma_\ast = 10$, whereas for the Cameraman image, we test with $\gamma = 10$, $20$ when $\gamma_\ast = 15$. It can be observed that even with the values for $\gamma$ that are different from the precise noise level $\gamma_\ast$, we can attain similar denoising results to the ones obtained with $\gamma_\ast$, especially when $\gamma > \gamma_\ast$. These indicate that although $\gamma = \gamma_\ast$ generates the best denoising results, the parameter $\gamma$ does not significantly affect the denoising results. The estimation for the parameter $\gamma$ in the presence of mixed Cauchy-Gaussian noise is another effortful task, so we leave this problem as a future work.

Lastly, in Table 5, the computational cost of the variational models is reported. The computational time of our model is the slowest among the models. Thus, there is a tradeoff between computing time and restoration performance. Despite the efficiency of our optimization algorithm, reducing the computational time of our model remains an issue.

5.   Conclusions

In this paper, we introduced a novel image denoising model under mixed Cauchy and Gaussian noise. The model is composed of a nonconvex data-fidelity term, expressed as an infimal convolution combination of two data-fidelity terms associated with two noise distributions, and total variation regularization. This new data-fidelity term enabled the separation of Cauchy noise and Gaussian noise components. It facilitated simultaneous removal of both noise. Total variation regularization assisted in the sufficient elimination of the mixed noise in homogeneous regions, while keeping structural edges and fine features. Regardless of the nonconvexity of the model, we proved the existence of a minimizer. To solve the proposed model, we utilized an alternating minimization approach and the alternating direction method of multipliers. This contributed to an iterative algorithm, and its convergence was shown experimentally. Numerical results demonstrated the effectiveness of the proposed model, comparing to other existing models, regarding both visual aspect and image quality assessments. However, theoretical convergence analysis of the proposed algorithm remains an issue. Blind denoising with unknown Cauchy noise level $\gamma$ is another demanding problem to be investigated in the future.

Acknowledgments

The author was supported by the Hankuk University of Foreign Studies Research Fund and the National Research Foundation of Korea (2021R1F1A1048111).

Conflict of interest

The authors declare no conflicts of interest in this paper.

Appendix A

For each fixed $x\in\Omega$, we define the function $h:\mathbb{R}\times \mathbb{R}\rightarrow \mathbb{R}$ as

The first-order and second-order partial derivatives of $h$ are obtained as

where $C$ is defined as

It can be easily proven that if $4\gamma^2\mu \geq \lambda_2$, then $\frac{\partial^2 h}{\partial s^2}(s, \cdot)\geq 0$, so $h$ is convex with respect to the variable $s$. Then, $\frac{\partial h}{\partial s}(s, \cdot) = 0$ is equivalent to a cubic equation that has either two real roots or one real root. Thus, the function $h(s, \cdot)$ has only one minimizer, so $h(s, \cdot)$ is strictly convex if $4\gamma^2\mu \geq \lambda_2$. Similarly, if $8\gamma^2\lambda_1 \geq \lambda_2$, then $\frac{\partial^2 h}{\partial t^2}(\cdot, t)\geq 0$, so $h$ is strictly convex with respect to $t$. Moreover, for the convexity of $h$ on $\mathbb{R}\times \mathbb{R}$, the Hessian matrix of $h$ needs to be positive semi-definite. That is, in addition to the aforementioned constraints, the determinant of the Hessian matrix of $h$ should be nonnegative:

This inequality holds if $\lambda_2(\mu+2\lambda_1)\leq 8\gamma^2\mu\lambda_1$, which is satisfied when $\lambda_2 \leq \frac{1}{2} \min \{ 8\gamma^2\lambda_1, 4\gamma^2\mu \}$. Thus, if $\lambda_2 \leq \frac{1}{2} \min \{ 8\gamma^2\lambda_1, 4\gamma^2\mu \}$, then the Hessian matrix of $h$ is positive semi-definte, so $h$ is convex. In addition, $h$ is strictly convex if $\lambda_2 < \frac{1}{2} \min \{ 8\gamma^2\lambda_1, 4\gamma^2\mu \}$.