Non-convex TV denoising corrupted by impulse noise

We propose a non-convex type total variation model for impulse noise removal by incorporating TV and the quasi-norm \begin{document}$\ell_q $\end{document} , \begin{document}$0 . Since the proposed model is non-convex and non-smooth, an iteratively reweighted algorithm is adapted and combined with a linearized ADMM. The convergence of the proposed algorithm is established and numerical results are given to illustrate the validity and efficiency of the proposed model.


1.
Introduction. We consider the image denoising problem degraded by impulse noise, which is caused by faulty pixels in camera or defective memory positions in hardware, for example [1,4]. Two types of impulse noise are widely considered; salt-and-pepper noise and random-valued noise. Salt-and-pepper noise takes only the maximum and the minimum values in its range at noisy pixels. In the case of random-valued noise, any random value can be taken [4].
To remove impulse noise, median type filters are widely used as a robust estimator, since noise statistics are characterized by probability densities having heavier tails than that of Gaussian [1]. To overcome the limitation of the median filter, diverse remedies and successors have been proposed such as the adaptive median filter [9] and the multistate median filter [6]. We refer the interested reader to the paper [4] and the references therein.
On the other hand, denoising models using the ℓ 1 norm as a data-fidelity term have been proposed to deal with outliers and impulse noise [18,4,14]. Actually the median and the ℓ 1 norm are closely related, since the median is a minimizer of the mean absolute error (MAE) and thus the sample median is the maximum likelihood (ML) estimator in the presence of Laplacian noise. To preserve edges and maintain sharp discontinuities, total variation (TV) is applied for the regularizer in addition to the ℓ 1 norm as fidelity term [3,18].
However, TV often causes staircase artifacts, favoring piecewise constant solutions even in smooth regions. To alleviate the drawback in smooth regions and to promote piecewise smooth solutions, higher-order TV is applied for Gaussian noise [5,13]. Since a heavy-tailed distribution is often observed in natural image gradients, a hyper-Laplacian distribution p(x) ∝ e −k|x| α with 0.5 ≤ α ≤ 0.8 is introduced [10]. Combining the idea of higher-order TV and hyper-Laplacian, we apply a non-convex type TV regularizer [8,15], which is composed of the ℓ q (0 ≤ q < 1) norm 1 and TV, to remove impulse noise.
Since the proposed model is non-convex and non-smooth, we adapt an iteratively reweighted algorithm, which is widely applied to the ℓ q norm minimization arising in compressed sensing [2,7]. We establish the convergence property of the proposed algorithm by showing that the values of the objective function are nonincreasing in iterates; see Theorem 3.1. Numerical experiments show the efficiency and stability of the proposed model in impulse noise reduction.
In our notation, for any x ∈ ℜ n , x j denotes the jth component of x, and x p = n j=1 |x j | p The identity matrix is denoted by I and the vector of zero entries is denoted by 0. Unless otherwise specified, {x (k) } denotes the sequence x (0) , x (1) , ....
where u is the underlying clean image and η denotes impulse noise.
As noticed in Introduction, to alleviate the drawback in smooth regions and cope with a heavy-tailed distribution in image gradients, we propose the nonconvex TV regularized ℓ 1 norm (NCTVℓ 1 ) model: where 0 < q 1 , q 2 < 1 and c 1 , c 2 ≥ 0. Here, we denote ∇u q1 Due to non-convexity and non-smoothness, it is a challenging minimization problem. Thus we apply the iteratively reweighted ℓ 1 algorithm, given in the next section.
3. Iteratively reweighted ℓ 1 algorithm for NCTVℓ 1 . In this section, we describe the iteratively reweighted ℓ 1 (IRℓ 1 ) algorithm for the NCTVℓ 1 denoising model (3). For simplicity of notation, we let b ∈ ℜ n + and u ∈ ℜ n + be vectorized versions of the two dimensional observed image and the reconstructed image of the n = M × N size, respectively.
To overcome the non-convexity of NCTVℓ 1 (3), the following convex approximation of the proposed model is solved at each iteration of IRℓ 1 : with η (k) ≥ 0. The formal procedure of IRℓ 1 is given in Algorithm 1.

Algorithm 1 IRℓ 1
Update u (k+1) and η (k+1) from u (k) and η (k) : Theorem 3.1 establishes the convergence property of IRℓ 1 . The proof is based on that the objective value decreases as the iteration of IRℓ 1 proceeds. In what follows,ũ is called a stationary point of  Proof. To clear up the proof, we define the function f as follows: and consider the following equivalent problem to (4): with t i = (∇u (k) ) i + η (k) , s i = (∇ 2 u (k) ) i + η (k) for all i and η = η (k) . By letting t for all i, we observe that Next, we consider the minimization problem with u = u (k+1) and η = η (k) . Then the solution is easy: The above inequalities and the equation (10) imply that The inequalities (6), (8), and (10) imply that the sequence {f (u (k) , t (k) , s (k) , η (k) ) + u (k) − b 1 } is non-increasing. By a similar procedure as in (10), we get This convergence property with the monotonicity of Also, from the inequality (6), we further have In addition, by the equivalence between (4) and (5) Taking limits on both sides of this inequality with K and (11) imply that Therefore, u * is a stationary point of the problem (3).
Although the problem (4) is convex, we further need an efficient algorithm, because of the nonlinearlity and non-smoothness of (∇u) i , (∇ 2 u) i , and u − b 1 . For this purpose, we apply linearized ADMM (LADMM); we introduce new variables v = ∇u, w = ∇ 2 u, z = u − b and then reformulate the subproblem (4) as a linearly constrained minimization problem: To describe the algorithmic framework of LADMM, we introduce the corresponding augmented Lagrangian function and linearized augmented Lagrangian function: where Finally, we describe LADMM using the augmented Lagrangian function (13) and the linearized augmented Lagrangian function (14) to solve the linearly constrained reformulation (12).
The first four minimzations are easy to solve, since we have the following closed form solutions.

v
where the shrinkage operator is adopted from [17] along with the convention 0 · (0/0) = 0. The convergence properties of LADMM for linearizing the quadratic penalty term can be found in [11].
Then the corresponding augmented Lagrangian function and linearized augmented Lagrangian function are as follows: The algorithmic framework to solve (16) is comparable to (15).
The first three steps have the following closed form solutions.

Numerical experiments.
In this section, we report numerical results on the proposed impulse noise removal model. The performance of NCTVℓ 1 (3) using IRℓ 1 is compared to these of TVℓ 1 (2), and the standard median filter. To show the efficiency of the model, we further compare it with the penalized model which is recently proposed in [12] and called AMM. This model uses the anisotropic TV and is solved by alternating minimization method. When the penalty parameters go to infinity, this model approaches (2) with replacing TV by the anisotropic TV. We use the relative difference for the stopping criterion of the proposed NCTVℓ 1 where C t > 0 is a stopping tolerance and set to C t = 5 × 10 −3 . For the inner loop, we choose the same criterion with C t = 1 × 10 −3 or we stop after 20 iterations. For quality evaluation, SNR and PSNR are employed, defined by SN R = 10 log 10 P s P n , P SN R = 10 log 10 255 2 P n , where P s = 1 n ||u|| 2 and P n = 1 n ||u −û|| 2 with the original and reconstructed images u andû, respectively.
To see the benefit of the NCTVℓ 1 model clearly, the proposed and two comparison models, the median filter and TVℓ 1 are applied to an one-dimensional piecewise smooth image under salt-and-pepper noise with various levels. The results are given in Figures 1 and 2. The parameters q 1 = q 2 = 0.9, c 1 = c 2 = 2.5, β = 1.5, γ = 0.01, δ = 2.5, and ρ = 8.0 are chosen for NCTVℓ 1 , and c = 2.5, β = 0.9, δ = 2.5, and ρ = 8.0 are for TVℓ 1 regardless of the noise levels. The NCTVℓ 1 model produces robust reconstruction in various noise levels, while two existing models suffer from the undesirable non-removed spikes and the number of those spikes increases as the noise level increases. The proposed model outperforms two models over 8 − 10dB for noise level 30% and about 10dB for 50%, in terms of SNR. Those models are applied to the four test images "lori", "boat", "cameraman", and "lighthouse" under salt-and-pepper noise and RVIN with various levels, given in Figures 3 -6. The model parameters q 1 , q 2 , c 1 , and c 2 are chosen by the best average PSNR on the chosen 3 training images ("babara", "boat", and "goldhill" of the size 512×512). The algorithm parameters β, γ, ρ, and δ are heuristically chosen and fixed for all experiments. The parameters are set to q 1 = q 2 = 0.9, c 1 = 0.3, c 2 = 0.4, β = 0.9, γ = 0.01, δ = 2.5, and ρ = 8.0 for NCTVℓ 1 , regardless of noise types and levels. For TVℓ 1 , the model parameter c is differently set to noise levels (c = 0.6 for 20%, c = 0.7 for 30%, c = 0.8 for 40%) for better performances and the others are fixed to β = 0.9, δ = 2.5, and ρ = 8.0. For the parameters of AMM, we refer to [12].
In Figure 3, TVℓ 1 yields non-smooth edges around hat and face line in the recovered image, while NCTVℓ 1 preserves smooth edges. For AMM, some unevenness is observed in the plain region. Figure 4 shows the reconstructed subimages of the "boat" image under the noise level 30%. The edges of poles on the boat and border lines of the boat body are non-smooth and rough in the reconstructions by TVℓ 1 . In Figure 5, TVℓ 1 and AMM provide jagged edges at the camera body and tripod regions. In contrast, the proposed NCTVℓ 1 shows well-preserved and smooth edges and consistent visual results in various noise levels. Figure 6 illustrates the staircase artifacts in the house and bottom regions of the recovered image by TVℓ 1 . The sky region is unevenly recovered in the AMM case.   Tables 1-3 show the performance comparison among three methods. They includes PSNR under various noise types and levels. In addition, we report the computational time, the number of iterations, and, if applicable, the total number of inner iterations. Overall, NCTVℓ 1 improves PSNR by 3 − 5dB comparing TVℓ 1 . Considering AMM, PSNR is 2 − 3dB higher. However, due to the complexity introduced by ℓ q norms and the Hessian term, NCTVℓ 1 is computationally slower than the TVℓ 1 model. 6. Conclusions. In this paper, we propose a non-convex variational model for impulse noise removal. The new model uses a non-convex regularizer which is composed of ℓ q norm and TV. To bypass non-convexity and non-smoothness, the iteratively reweighted ℓ 1 algorithm is applied. We show the convergence property of the proposed algorithm and demonstrate the effectiveness of the proposed model via numerical simulations.