Abstract
The variational models with nonlocal regularization offer superior image restoration quality over traditional method. But the processing speed remains a bottleneck due to the calculation quantity brought by the recent iterative algorithms. In this paper, a fast algorithm is proposed to restore the multichannel image in the presence of additive Gaussian noise by minimizing an energy function consisting of an -norm fidelity term and a nonlocal vectorial total variational regularization term. This algorithm is based on the variable splitting and penalty techniques in optimization. Following our previous work on the proof of the existence and the uniqueness of the solution of the model, we establish and prove the convergence properties of this algorithm, which are the finite convergence for some variables and the -linear convergence for the rest. Experiments show that this model has a fabulous texture-preserving property in restoring color images. Both the theoretical derivation of the computation complexity analysis and the experimental results show that the proposed algorithm performs favorably in comparison to the widely used fixed point algorithm.
1. Introduction
The restoration of multichannel image is one of the most important inverse problems in image processing. Among all the methods for restoration, the vectorial total variational (VTV) regularization is most popular for preserving discontinuities in restored images. VTV has been introduced by [1] into multichannel image restoration as a vectorial extension of total variation (TV) in gray scale image processing.
The famous total variational (TV) methods have drawbacks due to its local property while calculating the gradient of the image. Thus the images that contain many small structures or textures cannot get satisfying restored results. As a result, TV method based on nonlocal (NL) operator is generally used in gray scale image processing because of its structure preserving property. This method is generalized from the Yaroslavsky denoising and the block method. The main idea is using the similar blocks in the whole image to estimate current pixel value. It was introduced by Efros and Leung [2] while studying texture synthesis. Then Buades et al. [3] generalized nonlocal method and proposed the well-known neighbor denoising filter, namely, nonlocal means (NLM). After that, Gilboa and Osher [4] united nonlocal method into the regularization framework by defining the nonlocal formed TV model including NL-TV- model, NL-TV- model, and NL-TV-G model and proposed the algorithms separately. NL-TV model can both preserve the edges and the other small structures in images and has excellent application effect.
Our previews works extends NL-TV regularization to NL-VTV regularization for multichannel image restoration (see [5]). The corresponding algorithms are based on the variation calculus by solving the Euler-Lagrangian equations by fixed point iterations. This kind of methods is tested effectively but has a low convergence rate.
Recently, the optimization algorithms based on variable-splitting and penalty techniques have attracted much attention in image and signal processing, such as the dual methods, the Bregman methods, and the alternating direction method of multipliers (ADMM) methods (see [6–8]). Correspondingly, these methods have been extended for multichannel images. Bresson and Chan [6] has extended Chambolle’s dual algorithm for TV problem to vectorial images by defining the standard dual definition of VTV. Wang et al. [9, 10] has applied the alternating minimization algorithm to VTV model for multichannel image deblurring.
In this paper, we propose a fast alternating algorithm for NL-VTV model to restore color images from noisy observations. The following sections are organized as follows. Section 2 is the brief introduction of the NL-VTV model and classical fixed point algorithm. In Section 3, we propose the alternating algorithm for NL-VTV model. In Section 4, we do the convergence analysis of this algorithm. In Section 5, we applied this algorithm to color image denoising to present the texture-preserving effect of the NL-VTV model and the high convergent rate of the proposed algorithm.
2. Nonlocal Vectorial Total Variation Model
In this section, we will briefly introduce the nonlocal vectorial total variational model for multichannel image processing. The usually used iterative algorithm will be introduced as well, which will be compared to our algorithm in Section 5.
2.1. Vectorial Total Variation Model
Consider a vector valued function , such as the multichannel images, defined on a bounded open domain : where denotes the th channel of .
The regularized VTV model for multichannel image reconstruction could be written as where is the observed noisy multichannel image and is the regularized parameter. The vectorial gradient of is defined as We refer to [6] for more information about the well-posedness analysis of the VTV model.
2.2. Nonlocal Functional
The neighborhood filter NLM proposed by Buades et al. developed to a famous image denoising method [3]. It uses the pixels similar to the current pixel to estimate the value of the current pixel. Let the reference image be . The similarity of the two pixels in the image , could be measured as the following weight function : where is the Gaussian kernel with the standard deviation of : where is the filter parameter. Generally, corresponds to the noise level, which is usually set as the standard deviation of the noise. It is shown that the weight function above could be calculated by the image blocks centered around , . The value of the weight function reflects the similarity of the two blocks. The estimation of the value at pixel by the nonlocal means is where the normalized factor . It can be seen that this method is a kind of special filter algorithm, and the estimation of the value at pixel is weighted by the values of all the pixels in the image.
Let the 2-dimensional image domain be , pixels , and is a real function on . The derivation under the nonlocal frame is where is the nonnegative symmetrical weight function defined as (4), which is calculated by the reference image.
The nonlocal gradient is defined as a vector composed by all the partial derivations at position : , where :
Define the nonlocal divergence as which satisfies the conjugation with the nonlocal gradient; that is,
Thus the nonlocal Laplacian is defined as The NL-TV norm of based on nonlocal operators, that is, the weighted gradient , is defined by where function is defined as (4). The NL-TV model can be written as The well-posedness of NL-TV model is guaranteed in [11].
2.3. Nonlocal Vectorial Total Variation Model and the Fixed Point Iterative Algorithm
In this section, the nonlocal gradient will be introduced to reconstruct the variation norm in the vectorial total variation model to deduce the NL-VTV model.
The vectorial nonlocal gradient could be derived by generalizing the nonlocal gradient defined previously: where , , represents the gradient of each channel separately, which is calculated by the image of its own channel.
Use the vectorial nonlocal gradient to replace the gradient in VTV norm; we get the NL-VTV model: where is the VTV regularization based on the nonlocal gradient, which could be denoted by
The existence and the uniqueness of the solution to (15) have been proven in our previews work, in which we proposed a fixed point iterative algorithm to solve (15) as well. Here we will give out a brief introduction as follows.
The Euler-Lagrangian equation of the functional (15) is
Calculate the nonlocal curvature channel by channel; that is, where , . Then the Euler-Lagrange equation becomes
According to (19), when the weight function is fixed, the Euler-Lagrange equation is linear, which could be solved by the gradient descent kind of algorithm as follows.
First of all, we discretize (19). Denote to be the pixel value at position of the th channel image, where , (if the size of the image is , ). is the discrete value of the weight function , that is, the weight function value with respect to pixel and pixel . The pixel traverses the image channel where the pixel is in. Denote the neighbor set to be , where denotes a threshold value as a judgment of the similarity between the two blocks centered at pixel and pixel , respectively. Let be the discrete : And let be the discrete : .
As for functions, the discrete inner product is , while, as for vectors, the discrete dot product is , the discrete inner product is . The norm of vector is .
Then for channel , where , the discrete nonlocal curvature is where and . Then the Euler-Lagrange equation becomes
Denote . Then the Euler-Lagrange equation could be shortly written as Furthermore, it could be written as
According to the fixed point iteration, we use the multichannel image in step to update image in step . The weights are calculated by the original observation image. Then we get the iteration formula as where , .
The convergence property of the fixed point algorithm has been proved in our in-press work. However, when , the NLTV functional is nondifferentiable, which is similar to the classical TV model. Thus a regularization is involved in the calculation of to escape the denominator to be zero. Consequently, this algorithm takes many iteration steps to obtain the best restoration result. Thus we consider designing a fast alternating minimization algorithm instead.
3. The Alternating Minimization Algorithm
For the sake of conquering the nondifferentiability of the VTV functional, we design an alternating minimization algorithm in this section to solve the problem (15). The main idea of the alternating minimization is introducing an auxiliary variable. Then the original problem could be split into two subproblems of which one is differential with respect to the original variable, while the other is convex with respect to the auxiliary variable, respectively.
Yang et al. [10] have already studied the alternating minimization algorithm for a general model: Model (26) reduces to (15) when , , and , where is the identity matrix of order and represents the Kronecker product. Yang et al. pointed out in [10] that if is the finite difference operator , then, under some boundary conditions for and some other assumptions for the blurring kernel , the problem (26) can be fast transformed by fast Fourier transform, which diagonalizes the block circulant matrices: , , and their transpose, and then can be solved by the fast algorithm such as the Gaussian elimination method.
However, the nonlocal gradient operator in (15) is not block circulant. We will figure out another method to speed up the algorithm. We briefly propose the algorithm in this section. The convergence of this algorithm will be analysed in Section 4, as well as others properties of this algorithm.
3.1. Variable Splitting (Half-Quadratic Formulation of (15))
We introduce an auxiliary variable to approximate the nonlocal gradient of , that is, . Then the minimization problem (15) is transformed into an equivalent constrained problem:
In order to solve problem (27), we use the penalty method to reformulate it as where is the penalty parameter. It could be easily proved that, as , the solution of problem (28) is convergent to the solution of problem (15).
3.2. Alternating Minimization Algorithm
The problem of minimizing (28) with respect to for fixed becomes The functional in (29) is convex with respect to . According to Lemma 3.3 in [10], the problem (29) has the unique solutions given by where the multidimensional shrinkage operator .
On the other hand, minimizing (28) with respect to for fixed becomes The functional in (31) is differential with respect to . Thus the solution to (31) satisfies the first order optimized condition; that is, the Euler-Lagrange equation of the functional (31) is After some steps of simple calculation, we get Then according to (11), we get Finally, the multichannel image is calculated as follows:
Since the graph Laplacian is negative semidefinite, the operator is diagonally dominant. Therefore we can solve by a Gauss-Seidel algorithm. However, the subproblem (35) consists in solving the linear equations with the scale of for each image channel. As a result, the calculation complexity of the Gauss-Seidel algorithm is heavy. We are inspired by Chen and Cheng [12] to simplify, and, without loss of generality as well, we chose . Furthermore, using Taylor expansion to approximate, there is where is the weight matrix of similar blocks of each channel, which is normalized to a symmetrical matrix with the sum of every row equal to 1. Thus the solution to (35) is According to (9), the divergence at pixel is , .
To conclude this section, we rewrite the alternating minimization algorithm as in Table 1.
4. Convergence Analysis
Since it has been proved in our in-press work that the model (15) has the unique solution, in this section, we prove that the solution obtained by Algorithm 1 is the unique one. First, the convergence of the quadratic penalty method as the penalty parameter goes to infinity is well known [13]. That is, as , the solution of (31) converges to that of (15). We analyze the convergence properties of Algorithm 1 for a fixed . In this section, we show that, for any fixed , the algorithm of minimizing and alternately has finite convergence for some variables and -linear convergence for the rest.
| ||||||||||
Theorem 1. The iterative sequence generated by (30) and (37) from any starting point converges to the solution of the minimization problem (31).
Proof. According to (30) and (37), the functional in (28) satisfies
Therefore, the sequence is nonincreasing. Thus it is convergent in because of its nonnegativity as well. We denote the limit of the sequence by . We prove that followed, where and are the solution of (31).
Because is coercive and the sequence is convergent, the sequence is bounded. Thus there exists a subsequence , which converges to an accumulation point when . According to (30) and (37), for all , ,
and, for all , ,
We denote as the accumulation point of the sequence , because of that the functional is continuous, according to (38),
Because of that the threshold function, shrink, in (30), is continuous, we compute the limits of both sides of the equation and get
Thus . According to (41), . On the other hand, is strictly convex with respect to , then we have because of the uniqueness of the minimizer of . Similarly, we prove that .
We compute the limits of (39) and (40),
According to (43), . Thus the limit points of all of the subsequence of are the solutions to (28). As a result, converges to .
Let and let . We define set , where is a positive constant. Then we have the following conclusion.
Theorem 2 (-linear convergence). Assume that the sequence ; then there is a constant , such that for all , the sequence generated by (30) and (37) is q-linear convergence.
Proof. Combine (30) and (37); we have
We know that is a fixed point of (44). Let . Then
Since the operator is nonexpansive,
Because ,
We write . Then by uniting (46) and (47), we have
Let . Then for all , the sequence generated by (30) and (37) is -linear convergent to .
5. Experiments
In this section, we apply Algorithm 1 to color image denoising to test its efficiency. We compare this method to that in [10] to show the detail-preserving property of the nonlocal VTV model. Then we compare the convergent rate of this algorithm with the fixed point algorithm presented in Section 2.
5.1. Texture-Preserving Property of the NLVTV Model
We choose the color Barbara picture size and the color Baboon picture size for testing, both of which contain many textures, such as the woman’s scarf and the baboon’s beard. We add the Gaussian white noise with standard deviation 20 to the two original pictures. Then we use the VTV model and the NLVTV model to restore them, separately. In this experiment, the signal to noise rate (SNR) is used to estimate the restoration effect of the methods, which is defined as follows: where represents the variance of all the elements of matrix .
As is shown in Figure 1 and Table 1, the NLVTV model obtains similar SNR to VTV model. In the experiment of restoring the color Barbara picture, we choose the restored image by NLVTV model under the SNR (13.0678) lower than but close to that of the VTV model (SNR = 13.1212) instead of the restored image under the best SNR (13.2873). Figure 2 shows the results. Figure 3 shows the zoom part of the restored image, that is, the woman’s scarf. We can see that the NVTV model preserves the details better than the VTV model. Figures 4 and 5 show the results of the experiment of restoring the color Baboon picture. Although the NLVTV model gets lower SNR than the VTV model, the textures in this picture are better preserved by the NLVTV model rather than the VTV model.
(a) Original image
(b) Noisy image
(c) Denoised by VTV (SNR = 13.1212)
(d) Denoised by NLVTV (SNR = 13.0678)
(a) Original image
(b) Noisy image
(c) Denoised by VTV
(d) Denoised by NLVTV
(a) Original image
(b) Noisy image
(c) Denoised by VTV (SNR=13.1618)
(d) Denoised by NVTV (SNR=12.5936)
(a) Original image
(b) Noisy image
(c) Denoised by VTV
(d) Denoised by NVTV
5.2. Computation Complexity Analysis and the Convergent Rate Comparison
In this section, we compare the computation complexity and the convergent rate of this algorithm to the fixed point algorithm.
First, we analyze the computation complexities of the two algorithms, respectively. Assume the size of the multichannel image is . The size of the nonlocal searching window is . We do not consider the computation complexity of the calculation of the nonlocal weight, because that is executed exactly the same in both algorithms.
At each iterative step of Algorithm 1, namely, (30) and (37), the dominant calculation of (30) contains the calculation of the nonlocal gradient and the threshold computation of times. Thus the computation complexity is . The dominant calculation of (37) contains a step of nonlocal divergence computation with the computation complexity of and a matrix multiplication with the computation complexity of . As a result, the total computation complexity for each iteration is .
On the other hand, each iterative step of the fixed point iteration contains the nonlocal curvature with the computation complexity of and three matrix multiplications with the computation complexity of . As a result, the total computation complexity for each iteration is .
Consequently, the computation complexity of the proposed algorithm is much less than the fixed point algorithm.
The convergence rate can be represented by the relationship between the iterative step and the standard deviation of the restored image, intuitively. Figure 6 shows the comparison of the convergence rate between the fixed point iteration and the proposed alternating minimization algorithm. The horizontal axis represents the iterative step. And the vertical axis represents the standard deviation of the restored image in each step. The standard deviation is calculated as where and represent the restored images of step and step , respectively. When calculating , we reshape matrix into a vector whose elements are taken column-wise from matrix . Then Figure 6 shows that the alternating minimization algorithm converges faster.
6. Conclusion
In this paper, we propose a fast alternating minimization algorithm for multichannel image processing, based on the new nonlocal vectorial total variational model introduced recently. We use the variable splitting and penalty techniques to split the energy minimization problem into two optimization subproblems. One is a simple minimization problem, which is solved by the threshold method. The other is a simple minimization problem, which is solved by the variation calculus. As a result, this algorithm is faster than the classical algorithm because of avoiding dealing with the nondifferentiability of the nonlocal TV norm.
The main contribution of this paper is introducing the alternating minimization algorithm for the nonlocal vectorial total variational model, which is proved fast and can be applied to many multichannel image processing problems. In this paper, we apply this algorithm to the RGB color image denoising as an example. It can be easily extended to other vectorial problems, such as color image deblurring, and to other color image models, such as hue-saturation-value (HSV) model, as well.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This paper and the related works are supported by grants from the National Natural Science Foundation of China under Contracts 61072142, 61271437, and 61201337 and the Science Research Project of National University of Defense Technology under Contracts JC12-02-05 and JC13-02-03.