Cauchy Noise Removal by Nonconvex ADMM with Convergence Guarantees

Image restoration is one of the essential tasks in image processing. In order to restore images from blurs and noise while also preserving their edges, one often applies total variation (TV) minimization. Cauchy noise, which frequently appears in engineering applications, is a kind of impulsive and non-Gaussian noise. Removing Cauchy noise can be achieved by solving a nonconvex TV minimization problem, which is difficult due to its nonconvexity and nonsmoothness. In this paper, we adapt recent results in the literature and develop a specific alternating direction method of multiplier to solve this problem. Theoretically, we establish the convergence of our method to a stationary point. Experimental results demonstrate that the proposed method is competitive with other methods in visual and quantitative measures. In particular, our method achieves higher PSNRs for 0.5 dB on average.


Introduction
In many imaging applications, images inevitably contain natural non-Gaussian noises, such as impulse noise, Poisson noise, multiplicative noise, and Cauchy noise. At the same time, the images may have been blurred by the point spread function (PSF) during their acquisition. Therefore, the image restoration problem is an essential task. Researchers have proposed many methods to deblur and denoise images; see [12,16,17,27,35,36,41,54] and references therein. In this paper, we focus on recovering images corrupted by blurring and Cauchy noise. Cauchy noise usually arises in echo of radar, in the presence of low-frequency atmospheric noise, and in underwater acoustic signals [26,31,40]. According to [44,45], it follows Cauchy distribution and is impulsive.
We assume that the original gray-scale image u is defined on a connected bounded domain ⊂ R 2 with a compacted Lipschitz boundary. The observed image with blurs and Cauchy noise is given as follows: where f ∈ L 2 ( ) denotes the observed image, K ∈ L(L 1 ( ), L 2 ( )) represents a known linear and continuous blurring (or convolution) operator, and η ∈ L 2 ( ) denotes Cauchy noise. Our goal is to recover u from the observed image f . In recent years, much attention has been given to Cauchy noise removal, and several methods have been proposed. In [13], the authors applied a recursive algorithm based on the Markov random field to reconstruct images and retain sharp edges. In 2005, Achim and Kuruoǧlu utilized a bivariate maximum a posteriori estimator (BMAP) to propose a new statistical model in the complex wavelet domain for removing Cauchy noise [1]. In [34], Loza et al. proposed a statistical approach based on non-Gaussian distributions in the wavelet domain for tackling the image fusion problem. Their method achieved a significant improvement in fusion quality and noise reduction. In [46], Wan et al. developed a novel segmentation method for RGB images that are corrupted by Cauchy noise. They combined statistical methods with denoising techniques and obtained a satisfactory performance. Since TV regularization is able to preserve edges effectively while still suppressing noise satisfactorily [21], Sciacchitano et al. proposed a convex TV-based variational method for recovering images corrupted by Cauchy noise in [42]. The variational model in this method is as follows: where γ > 0 is the scale parameter of Cauchy distribution, and BV ( ) is the space of functions of bounded variation. Here, u ∈ BV ( ) if u ∈ L 1 ( ) and its total variation (TV) is finite, where (C ∞ 0 ( )) 2 is the space of vector-valued functions with compact support in . The space BV ( ) endowed with the norm u BV ( ) = u L 1 ( ) + |Du| is a Banach space; see, e.g., [21]. In (2), λ denotes the positive regularization parameter, which controls the trade-off between TV regularization and the fitting to f andũ,ũ is the result obtained by the median filter, and α is a positive penalty parameter. Note that if 8αγ 2 ≥ 1, the objective functional in (2) is strictly convex and leads to a unique solution. Because of strict convexity, the model avoids the common issues of nonconvex optimization: the solutions depend on the numerical methods and how they are initialized. But the last term in (2) in fact pushes the solution close to the median filter result, and the median filter does not always provide satisfactory removals of Cauchy noise. Hence, in this paper we turn our focus back to a nonconvex model.
Recently, researchers have discovered some useful convergence properties of the optimization algorithms for solving nonconvex minimization problems [24,47,48,53]. In particular, the paper [48] established the global convergence (to a stationary point) of the alternating direction method of multipliers (ADMM) for nonconvex nonsmooth optimization with linear constraints. To take advantages of the recent results, in this paper we develop the ADMM algorithm to solve the following nonconvex variational model directly for denoising and deblurring simultaneously: We prove that our algorithm starting from any initialization is globally convergent to a stationary point under certain conditions. Furthermore, we compare our proposed method to the state-of-the-art method proposed in [42] and show the effectiveness of our method in terms of restoration quality and noise reduction. The outline of the paper is summarized as follows. In the next section, we analyse some fundamental properties of Gaussian distribution, Laplace distribution and Cauchy distribution. In Sect. 3, we illustrate the nonconvex variational model for denoising and deblurring, and prove the existence and uniqueness of the solution. In Sect. 4, we develop our algorithm for the proposed nonconvex model and present the convergence results. In Sect. 5, we demonstrate the performance of our algorithm by comparing with other existing algorithms. Finally, we conclude the paper with some remarks in Sect. 6.

Statistical Properties for Cauchy Distribution
The Cauchy distribution is a special kind of the α-stable distributions with α = 1 and is important as a canonical example of the "pathological" case [3,15,29]. It is closed under linear fractional transformations with real coefficients [30]. However, different from most α-stable distributions, it possesses a probability density function that can be expressed analytically [19,28] as: where the parameter μ specifies the location of the peak and the parameter γ > 0 decides the half-width at half-maximum. Here, we let C(μ, γ ) denote the Cauchy distribution. Its mode and median are both μ while the mean, variance, and higher moments are undefined. In addition, the Cauchy distribution is infinitely divisible, that is, for every positive integer n, there exist n independent identically distributed (i.i.d.) random variables X n1 , X n2 , . . . X nn such that X n1 + X n2 + · · · + X nn follows the Cauchy distribution. Due to their infinite divisibility, random variables following the Cauchy distribution obey the generalized central limit theorem [37].
The Cauchy distribution is closely related to some other probability distributions. The Cauchy distribution is heavy-tailed, and its tail's heaviness is determined by the scale parameter γ . In particularly, if X and Y are two independent Gaussian random variables with mean 0 and variance 1, then the ratio X/Y follows the standard Cauchy distribution C(0, 1) [6,38]. In Sect. 5, we will apply this property to simulate images corrupted by Cauchy noise.
Further to show the statistical properties of the Cauchy distribution, we compare it with two most commonly used probability distributions: the Gaussian distribution (N (μ, σ 2 ) with mean μ and variance σ 2 ) and the Laplace distribution L(μ, b) with mean μ and variance 2b 2 . Since the Gaussian and Cauchy distributions are α-stable distributions with α = 2 and α = 1, respectively, they are both bell-shaped. Moreover, we can easily obtain the following relation between them at x = 0. Proposition 2.1 Let X 1 and X 2 be two independent random variables. Assume that X 1 ∼ N (0, 1) and X 2 ∼ C(0, 2 π ). Then the values of their probability density functions (PDFs) at x = 0 are equal.
In addition, both the Laplace and Cauchy distributions are heavy-tailed distributions. We demonstrate their relation by the tails of their distribution curves in the following proposition. Proposition 2.2 Let P G , P L and P C denote the PDFs for N (0, σ 2 ), L(0, b) and C(0, 2 π ), respectively. Then, the followings hold: 1. At x = σ = b = γ , the ratio of P G , P L and P C is 1 : π 2e : e 2π ; 2. At x = 3σ = 3b = 3γ , the ratio of P G , P L and P C is 1 : π 2 e Based on Proposition 2.2, we can see that the probability density value of the Gaussian distribution at a rather small x, saying x = σ = b = γ , is the largest, which shows that the additive Gaussian noise tends to mainly produce small perturbations. However, at larger x, saying x = 3σ = 3b = 3γ , the density of the Laplace distribution is more than 5 times that of the Gaussian distribution, and the density of the Cauchy distribution is even more than 7 times. Hence, the Laplace and Cauchy distributed additive noise tend to corrupt images with high perturbations. Figure 1 depicts the PDFs of the Gaussian, Laplace, and Cauchy distributions. From Fig. 1a, we see that these three distributions have different behaviours at the peaks and tails. See the details in the zoom-ins. Figure 1b depicts the portion around the peaks of the three distributions. The Gaussian distribution has the same peak as the Cauchy distribution while the density of the Gaussian distribution is slightly higher on both sides of the peak.
(c) Fig. 1 Comparison for probability density functions of N (0, 1), L(0, 2 π ) and C(0, 2 π ). a The plots of three distributions, b the zoomed-in portion of the curves around the peaks, c the zoomed-in portion of the curves around the tails distribution is closer to that of the Cauchy distribution than Gaussian distribution, but there still exists a big gap between the densities of the Laplace and Cauchy distributions. Therefore, the Cauchy distribution cannot be simply replaced with the Gaussian or Laplace distribution during image restoration.

Nonconvex Variational Model
This section describes our model of deblurring and denoising. In [42], a variational model for denoising was proposed. To make our exposition self-contained, we deduce a similar nonconvex variational model for deblurring and denoising based on the maximum a posteriori (MAP) estimator and Bayes' rule.

Nonconvex Variational Model Via MAP Estimator
We consider f (x) and u(x) as random variables for each x ∈ . The MAP estimator of u is the most likely value of u given f , i.e., u * = arg max u P(u| f ). Based on Bayes' rule and the independence of u(x) and f (x) for all x ∈ , we obtain arg max where the term log P( f (x)|u(x)) describes the degradation process that produces f from u based on (1), and log P(u) is the prior on u.
In addition, we use the prior P(u) = exp(− 2 λ |Du|). Then, we arrive at the variational model for deblurring and denoising: where λ > 0 is the regularization parameter. Although |Du| is convex, due to the logarithm in the data-fitting term, log γ 2 + (K u − f ) 2 dx is nonconvex. Therefore, the numerical solution of (5) depends on the numerical approach and how it is initialized.

Solution Existence and Uniqueness of the Model (5)
According to the properties of the total variation, we prove that there exists at least one solution for the nonconvex variational problem in the BV space.

Theorem 3.1 Assume that is a connected bounded set with compacted
Lipschitz boundary and f ∈ L 2 ( ). Suppose that K ∈ L(L 1 ( ), L 2 ( )) is nonnegative and linear with K 1 = 0. Then the model (5) has at least one solution u * ∈ BV ( ).
Obviously, E(u) is bounded from below. For a minimizing sequence {u k }, we know that E(u k ) is bounded, so both |Du k | and log γ 2 + (K u k − f ) 2 dx are bounded. Now we apply proof by contradiction to show that {K u k } is bounded in L 2 ( ) and therefore also bounded in L 1 ( ). Assume that K u k 2 = +∞, so there exists a set E ⊂ , whose measure is not zero, such that for any x ∈ E we have K u k (x) = +∞. Then, with f ∈ L 2 ( ) we will also have log

Based on
|Du k | being bounded, by the Poincaré inequality [2], we have where m (u k ) = 1 | | u k dx, C is a positive constant, and | | represents the measure of . As is bounded, Therefore, there exists a subsequence {u n k } in BV ( ) that converges strongly in L 1 ( ) to some u * ∈ BV ( ) as k → ∞, while {Du n k } converges weakly as a measure to Du * . Since K is linear and continuous, {K u n k } converges strongly to K u * in L 2 ( ). By the lower semicontinuity of total variation and Fatou's lemma, we conclude that u * is a solution of the model (5).
Although the objective function in (5) is nonconvex, we are still able to obtain a result on the uniqueness of the solution.
Theorem 3.2 Assume that f ∈ L 2 ( ) and K is injective. Then, the model (5) has a unique solution u * in U : Proof For each fixed x ∈ , we define a function g : R → R: Since the second order derivative of g: By the convexity of TV and linearity of K , the objective function of the model (5) is strictly convex in U . Hence, there exists a unique solution for the model (5) in U .
Note that Cauchy noise is so impulsive that, even with a small γ , many points in f are still heavily corrupted and thus some impulsive noise is still left in the images in U . If we also take the smoothing property of K into account, then the unique solution in U will not be satisfactory. In Sect. 5.1, we will demonstrate this point numerically.

Proposed ADMM Algorithm
Due to the nonconvexity of the variational model (5), different numerical algorithms and initializations may yield different solutions. Taking advantage of the recent result in [48], in this section we apply the ADMM algorithm to the minimization problem (5), which restores images degraded by blurring and Cauchy noise. Then, we prove that the proposed algorithm is globally convergent to a stationary point.

The ADMM Algorithm for Nonconvex and Nonsmooth Problem
We briefly review the ADMM algorithm and its recent convergence result under nonconvexity and nonsmoothness. Let We consider the minimization problem formulated as: min In general, F can be nonsmooth and nonconvex, and G can be nonconvex (but is differentiable as stated). By introducing a Lagrangian multiplier w ∈ R M for the linear constraint Ax + By = 0, we obtain the augmented Lagrangian: where β > 0 is a penalty parameter.
Extending from the classic ADMM [7,20], the multi-block ADMM generates the iterates where we use ∈ arg min when minimizers are not necessarily unique (in which case, any minimizer is fine). The general assumption is that all subproblems have minimizers. The convergence result of the ADMM algorithm under nonconvexity and nonsmoothness is summarized as follows [48]. We present the conditions that are simplified to fit our need yet more restrictive than those in [48].
Also, assume that A, B have full column rank 1

and I m(A) ⊂ I m(B). Further assume that
F (x) is either restricted prox-regular 2 or piecewise linear, and G(y) is Lipschitz differentiable with the constant L ∇G > 0. Then, for any β larger than a certain constant β 0 and starting from any initialization (x 0 , y 0 , w 0 ), ADMM (8) produces a sequence of iterates that has a convergent subsequence, whose limit is a stationary point (x * , y * , w * ) of the augmented Lagrangian L β (x, y; w). If in addition L β satisfies the Kurdyka-Łojasiewicz (KL) inequality [4,8,32,33], then the result improves to global convergence to that stationary point.

The ADMM Algorithm for Solving (5)
Taking advantage of the ADMM convergence result, we apply it to solve the nonconvex variational model in (5) for simultaneous denoising and deblurring. Hereafter, we switch to the discrete form, but, for the sake of simplicity, we still use the same letters defined in the continuous context. We assume that the discrete image domain contains n × n pixels. The discrete minimization nonconvex model of (5) is formulated as follows: where f ∈ R n 2 is obtained by stacking the columns of the corresponding n × n gray-scale image, and K ∈ R n 2 ×n 2 . The TV regularization ∇u 1 is defined as: where ∇ x ∈ R n 2 ×n 2 and ∇ y ∈ R n 2 ×n 2 are the discrete first order forward differences in the x-and y-directions, respectively. The discrete gradient of u, ∇u, is defined as ∇u = [(∇ x u) , (∇ y u) ] ∈ R 2n 2 .
To derive the ADMM algorithm for our model, we introduce a new auxiliary variable v ∈ R n 2 and obtain the following constrained nonconvex minimization problem: min u,v∈R n 2 Let w ∈ R n 2 be the Lagrangian multiplier for the constraint K u = v. Then we have the corresponding augmented Lagrangian: where β > 0 is a penalty parameter. The whole algorithm for restoring the blurred images corrupted by Cauchy noise is given in Algorithm 1.
In Algorithm 1, the dominant computation is the steps to solve the two minimization subproblems in (11) and (12). The u-subproblem (11) can be efficiently solved by many 2 A function h : R N → R is restricted prox-regular, if for any sufficiently large M ∈ R + and any bounded set T ⊂ R N , there exists τ > 0 such that 3: If u k+1 satisfies the stopping criteria, return u k+1 and stop.
In addition, taking some specific properties of the variational model (9) into account, we provide a relatively simple proof.

Theorem 4.2
Let (u 0 , v 0 , w 0 ) be any initial point and {(u k , v k , w k )} be the sequence of iterates generated by Algorithm 1. Then, if β > λ γ 2 and K has full column rank, the sequence {(u k , v k , w k )} converges globally to a point (u * , v * , w * ), which is a stationary point of L β .
In order to prove Theorem 4.2, based on the model in (7), we define the following functions: The feasible set is F = {(u, v) ∈ R n 2 ⊗ R n 2 : K u − v = 0}. First, we give some useful lemmas that will be used in the main proof.

Lemma 4.1
The iterates of Algorithm 1 satisfy: Proof Substituting (13) on w k into the first-order optimality condition of the v-subproblem Since G is smooth, we can calculate its second derivative and thus L ∇G = λ γ 2 is a Lipschitz constant for ∇G. Consequently, we obtain the bound Proof According to the optimality condition of the u-subproblem (11), we define From (11) and the definition of subgradient F , it follows where the second equality follows from the cosine rule:

b − a and the last inequality follows from the convexity of F (u).
For the updates of v k+1 , w k+1 , by the cosine rule and Lemma 4.1, we have where we have applied the inequality In order to ensure C > 0, we need the penalty parameter β to satisfy: According to (17) and (18), we have This means that L β (u k , v k , w k ) is nonincreasing in k ∈ N.
As K has full column rank, there existsv such that K u k −v = 0. Therefore, we have Thus, we arrive at By the definition of the subgradient, we have Thus, according to the optimal condition 0 ∈ ∂F(u k+1 ) and combining (19), (20), (21), and Lemma 4.1, we arrive at Now we give the proof to our main convergence theorem.
Proof of Theorem 4.2 As K has full column rank, the feasible set F is nonempty. By Lemma 4.2, the iterative sequence {(u k , v k , w k )} is bounded, so there exists a convergent subsequence {(u n k , v n k , w n k )}, i.e., (u n k , v n k , w n k ) converges to (u * , v * , w * ) as k goes to infinity. Since L β (u k , v k , w k ) is nonincreasing and lower-bounded, we have K (u k −u k+1 ) → 0 and v k − v k+1 → 0 as k → ∞. According to Lemma 4.3, there exists p k ∈ ∂L β (u k , v k , w k ) such that p k → 0. Further, this leads to p n k → 0 as k → ∞. Based on the definition of the general subgradient [39], we obtain that 0 ∈ ∂L β (u * , v * , w * ), i.e., (u * , v * , w * ) is a stationary point.
Referring to [47,51], the function F (u) is semi-algebraic, and G(v) is a real analytic function. Thus, we conclude that L β satisfies the KL inequality [8]. Then, as in the proof of Theorem 2.9 in [5], we can deduce that the iterative sequence {(u k , v k , w k )} is globally convergent to (u * , v * , w * ).
Remark 1 In Theorem 4.2 we need K to have full column rank. Since K is a blurring matrix in our problem, this requirement is satisfied.

Numerical Experiments
In this section, we present the results of several numerical experiments to demonstrate the performance of the proposed method for restoring images corrupted by blurs and Cauchy noise. Here, we use ten 8-bit 256-by-256 gray-scale test images, see Fig. 2. All numerical results are performed under Windows 10 and Matlab Version 7.10 (R2012a) running on a Lenovo laptop with a 1.7 GHz Intel Core CPU and 4GB RAM.
We utilize the peak signal-to-noise ratio (PSNR) and the structural similarity index (SSIM) [49] as performance measures, which are respectively defined as PSNR = 20 log 10 255n , whereũ is the restored image, u is the original image, μũ and μ u denote their respective means, σ 2 u and σ 2 u represent their respective variances, σ is the covariance ofũ and u,  . 3 Comparison of different initial values for removing Cauchy noise in the image "Parrot", with γ = 5 (in the 1st row) and 10 (in the 2nd row). a Noisy images, b restored images of (I), c restored images of (II), d restored images of (III) and c 1 , c 2 > 0 are constants. PSNR is a good measure of the human subjective sensation, and a higher PSNR implies better quality of the restored image. SSIM conforms with the quality perception of the human visual system (HVS). If the SSIM value is closer to 1, the characteristic (edges and textures) of the restored image is more similar to the original image. In our method, we set the stopping condition based on the following relative improvement inequality: where E is the objective function in (9) and = 5×10 −5 . In addition, since the regularization parameter λ balances the trade-off between fitting f and TV, we manually tune it in order to obtain the highest PSNRs of the restored images. The selection method of λ is out of the scope in this paper. The parameter β in Algorithm 1 affects the convergent speed. Based on Theorem 4.2, we round β > λ γ 2 up to the nearest value with two digits after the decimal point as β. In addition, we set the iteration number for the Newton method while solving the  The largest values are given in bold v-subproblem as 3. The iteration number for solving the u-subproblem equals 5 in denoising and 10 in simultaneous deblurring and denoising.

Different Initializations
Since our model (9) is nonconvex, though we are able to prove that the ADMM algorithm converges globally to a stationary point from any given starting point (u 0 , v 0 , w 0 ), the local minimizers that we obtained may still depend on the initial points. To study the influence of initializations and obtain better restorations, in this section we test three different choices of u 0 in denoising:  In Table 1, we list PSNRs and SSIMs for different initial points for the test images "Parrot" and "Cameraman" at the noise levels γ = 5 and 10. The noisy images are obtained via f = u + γ η 1 η 2 , where η 1 and η 2 are independent random variables following the Gaussian distribution with mean 0 and variance 1. It is obvious that both PSNRs and SSIMs are highest in the case (I), and are lowest in the case (III), which shows that the unique solution in U is not a satisfactory local minimizer. Figure 3 depicts the restored "Parrot" images in order to compare the visual performance due to different initial points. Figure 3d shows the unique solution in U , and we can see that there is still some noise left in the restored images. The reason is that Cauchy noise is so impulsive that by corrections in a small range, [−γ, γ ], it is not enough to remove all noise. Compared with the results from (II), the ones from (I) include clearer features and less noise, especially in the region around the eye and black stripes of "Parrot". Hence, we choose (I) as initialization in our remaining numerical experiments.
Theorem 4.2 demonstrate that with any given initial points, Algorithm 1 converges globally to a stationary point. Figure 4 depicts the plots of the objective function values in (9) versus the number of iteration in order to observe the convergence of our method. It is clear that the objective function value keeps decreasing over the iterations. Furthermore, our method converges very fast except in case (III), which does not provide good restorations.

Comparisons of Image Deblurring and Denoising
In order to demonstrate the superior performance of our proposed method, we compare it with two other well-known methods: the median filter (matlab function 'medfilt2') with window size 3 and the convex variational method in [42] ("conRE" for short). For fair comparison, we use the same stopping rule in the convex variational method and adjust the two parameters in the model for highest PSNRs. First, we compare the three methods for Cauchy noise removal, i.e., by setting K as the identity matrix. Tables 2 and 3 list the PSNRs and SSIMs of the restored images at the noise levels γ = 5 and γ = 10, respectively. Obviously, comparing to the two variational methods, the median filter provides the worst PSNRs and SSIMs. Our method always yields the highest PSNRs. Especially at the lower noise level (γ = 5), our PSNRs are about 1dB higher than the convex method [42]. Furthermore, in most cases, our SSIMs are also higher than others.
In Figs. 5 and 6, we present the results from different methods for removing Cauchy noise from the images "Parrot", "Cameraman", "Baboon", "Boat" and "Plane". Although the median filter effectively removes Cauchy noise, it also oversmooths the edges and destroys many details. It is obvious that two variational methods outperform the median filter. Comparing to the convex method, our nonconvex method can provide better balance between preserving detail and removing noise. To further illustrate the performance of our method, we show the zoomed regions of the restored images "Parrot", "Baboon" and "Boat" in Figs. 7 and 8, where we can clearly see the difference among the results from the three methods, Fig. 8 Zoomed version of the restored images in Fig. 6. a Original images; b the median filter; c the "conRe" model; d our method e.g., the stripes around the eye in "Parrot", the nose and whiskers of "Baboon", and the ropes and iron pillars of "Boat".
In the following experiments, we compare the three methods on recovering images corrupted by blurs and Cauchy noise. Here, we consider the Gaussian blur with size 7 and standard deviation 3, and the out-of-focus blur with size 5. Further, Cauchy noise with γ = 5 is added into the blurry images. Tables 4 and 5 list the PSNRs and SSIMs by applying different methods to the images "Parrot", "Cameraman", "Plane" and "Test". Figures 9 and 10 show the restored images. The largest values are given in bold The largest values are given in bold Fig. 9 Comparison of the restored results by applying different methods for deblurring and denoising the images degraded by a Gaussian blur (G, 7, 3) and Cauchy noise (γ = 5). a Degraded images; b the median filter; c the "conRe" model; d our method From Tables 4 to 5, we find that our method provides the highest PSNRs and SSIMs. Comparing with the convex method, our method can improve by at least 0.36dB on PSNR. In Figs. 9 and 10, it is easy to see that the restored images by the median filter are oversmoothing as the median filter does not deblur. The convex method can recover edges and textures, but some noise is not removed. However, our method not only preserves the fine features but also effectively removes Cauchy noise, which can be clearly seen in the zoomed regions in Fig. 11.  Fig. 10 Comparison of the restored results by applying different methods for deblurring and denoising the images degraded by the out-of-focus blur (A, 5) and Cauchy noise (γ = 5). a Degraded images; b the median filter; c the "conRe" model; d our method Fig. 11 Zoomed version of the restored results for the image "Parrot" degraded by the Gaussian blur (in the 1st row) and the out-of-focus blur (in the 2nd row), respectively. a Degraded images; b the median filter; c the "conRe" model; d our method

Conclusion
In this paper, we have reviewed and analyzed the statistic properties of the Cauchy distribution by comparing it with the Gaussian and Laplace distributions. Based on the MAP estimator, we have developed a nonconvex variational model for restoring images degraded by blurs and Cauchy noise. Taking advantage of a recent result in [48], the alternating direction method of multiplier (ADMM) algorithm is applied to solve the nonconvex variational optimization problem with a convergence guarantee. Numerical experiments show that the proposed method outperforms two well-known methods in both qualitative and quantitative comparisons.