Image Restoration Based on the Hybrid Total-Variation-Type Model
We propose a hybrid total-variation-type model for the image restoration problem based on combining advantages of the ROF model with the LLT model. Since two -norm terms in the proposed model make it difficultly solved by using some classically numerical methods directly, we first employ the alternating direction method of multipliers (ADMM) to solve a general form of the proposed model. Then, based on the ADMM and the Moreau-Yosida decomposition theory, a more efficient method called the proximal point method (PPM) is proposed and the convergence of the proposed method is proved. Some numerical results demonstrate the viability and efficiency of the proposed model and methods.
Image restoration is of momentous significance in coherent imaging systems and various image processing applications. The goal is to recover the real image from the deteriorated image, for example, image denoising, image deblurring, image inpainting, and so forth; see [1–3] for details.
For the additive noisy image, many denoising models have been proposed based on PDEs or variational methods over the last decades [1–3]. The essential idea for this class of models is to filter out the noise in an image while preserving significant features such as edges and textures. However, due to the ill-posedness of the restoration problem, we have to employ some regularization methods  to overcome it. The general form of regularization methods consists in minimizing an energy functional of the following form: in Banach space , where is the regularization parameter, , is a regularization term, is the observed image, and is the image to be restored. The development of the energy functional (1.1) actually profits from the ROF model  which is of the following form: where . In the case without confusion, for simplification we omit the open set with Lipschitz boundary for and . Due to edge-preserving property of the term , this model has been extended to other sorts of image processing problems such as for image deblurring  and for image inpainting [2, 7], where is a blurring operator and is the inpainting domain. Furthermore, this model was applied to restore the multiplicative noisy image which usually appears in various image processing applications such as in laser images, ultrasound images , synthetic aperture radar (SAR) , and medical ultrasonic images . One of the models based on was proposed by Huang et al. (HNW)  with the form where is the logarithmic transformation of and keeps the total variation property related to and . Here satisfies for the noise .
However, as we all know the term usually reduces the computational solution of the above models to be piecewise constant, which is also called the staircasing effect in smooth regions of the image. The staircase effect implies to produce new edges that do not exist in the true image so that the restored image is unsatisfactory to the eye. To overcome this drawback, some high-order models [12–15] have been proposed such as a model proposed by Lysaker, Lundervold, and Tai (the LLT model) with the following form: where for the Hessian operator . However, these classes of models can blur the edges of the image in the course of restoration. Therefore, it is a natural choice to combine advantages of the ROF model and the LLT model if we want to preserve edges while avoiding the staircase effect in smooth regions. One of convex combinations between the and was proposed by Lysaker et al.  to restore the image with additive noise. But their model is not quite intuitive due to lack of gradient information in the weighting function. Since the edge detector function with can depict the information of edges, we can employ it as a balance function; that is, we can apply the following model: to restore the noisy image . Obviously, tends to be predominant where edges most likely appear and tends to be predominant at the locations with smooth signals. Based on the advantages of the hybrid model (1.7), we also extend it to the image restoration models (1.3)–(1.5) in this paper.
Another topic for image restoration is to find some efficient methods to solve the above proposed models. In fact, there are many different methods based on PDE or convex optimization to solve the minimization problem (1.1) by means of the specific models (1.2)–(1.6). For example, for the purpose of solving the ROF model (1.2) or the LLT model (1.6), we can use the gradient descent method [5, 13], the Chambolle’s dual method [13, 17, 18], the primal and dual method [19–21], the second order cone programming method , the multigrid method , operator splitting method [24–26], the inverse scale method , and so forth. However, different from the ROF model (1.2) and the LLT model (1.6), the model (1.7) includes two -norm terms which make it solved more difficultly. More generally, the model (1.7) can be fell into the following framework: where are proper, convex, and lower semicontinuous (l.s.c.), and are bounded linear operators, and is a parameter. A specific form of (1.8) was considered by Afonso and Bioucas-Dias in  where a Bregman iterative method was proposed to solve the model with the combination of the -norm term and the total variation (TV) term. Actually this splitting Bregman method is formally equivalent to the alternating direction method of multipliers (ADMM) [24, 29–34]. However, the ADMM ineluctably tends to solve some subproblems which correspond to the related modified problems. Furthermore, these make us obtain the numerical results by requiring much more computational cost. In order to obtain an efficient numerical method, it is a natural choice to avoid solving the related subproblems. In this paper, we propose a proximal point method (PPM) which can be deduced from the ADMM. This deduction is based on the connection that the sum of the projection operator and the shrinkage operator is equal to the identity operator; it is known as the Moreau-Yosida decomposition Theorem 31.5 in . Then the PPM not only keeps the advantages of the ADMM but also requires much less computational cost. This implies that the PPM is much more fast and efficient, especially for the larger scale images. Furthermore, using the monotone operator theory, we give the convergence analysis of the proposed method. Moreover, we extend the PPM to solve the model (1.7) to image deblurring, image inpainting, and the multiplicative noisy image restoration. Experimental results show that the restored images generated by the proposed models and methods are desirable.
The paper is organized as follows. In Section 2, we recall some knowledge related to convex analysis. In Section 3, we first propose the ADMM to solve the problem (1.8) and then give the PPM to improve this method. In Section 4, we give some applications by using the proposed algorithms and also compare the related models and the proposed methods. Some concluding remarks are given in Section 5.
2. Notations and Definitions
Definition 2.1. The operator is monotone if it satisfies for all and , and is maximal if it is not strictly contained in any other monotone operator on .
Definition 2.2. Let be a convex function. The subdifferential of at a point is defined by where is called a subgradient. It reduces to the classical gradient if is differentiable.
Definition 2.4. Let be a convex proper function. The conjugate function of is defined by for all .
Definition 2.5. Let be convex and . The proximal mapping to is defined by for .
Definition 2.6. The projection operator onto the closed disc is defined by where and .
Definition 2.7. The shrinkage operator is defined by where we use the convention .
Remark 2.8. It is obvious that the function and its conjugate function satisfy the following relationship: for . Especially, the projection operator and the shrinkage operator satisfy for any . In fact, this corresponds to the classic Moreau-Yosida decomposition Theorem 31.5 in .
3. The Alternating Direction Method of Multipliers (ADMM) and the Proximal Point Method (PPM)
Variable splitting methods such as the ADMM [29–32, 44] and the operator splitting methods [42, 45, 46] have been recently used in the image, signal, and data processing community. The key of this class of methods is to transform the original problem into some subproblems so that we can easily solve these subproblems by employing some numerical methods. In this section we first consider to use the ADMM to solve the general minimization problem (1.8). However, the computational cost of the ADMM is tediously increased due to its looser form. In order to overcome this drawback, we thus change this method to a compacter form called the proximal method based on the relationship (2.9) in Remark 2.8.
3.1. The Alternating Direction Method of Multipliers (ADMM)
We now consider the following constrained problem: which is clearly equivalent to the constrained problem (1.8) in the feasible set . Throughout the following subsections, we always assume that and are a surjective map. It seems that the problem (3.1) including three variables looks more complex than the original unconstrained problem (1.8). In fact, this problem can be solved more easily under the condition that and are nondifferentiable. In the augmented Lagrangian framework [33, 34, 47, 48], the problem (3.1) equivalently solves the following minimization Lagrangian function: where is the Lagrangian multiplier and is the penalty parameter for , . Then we can use the following augmented Lagrangian method (ALM): with choosing the original values and to solve (3.2). If we set for , and omit the terms which are independent of in (3.3a), the above strategy (3.3a)–(3.3c) can be written as for the original values and . Then we can use the following ADMM to solve (3.4a)–(3.4c).
(3) If the stop criterion is not satisfied, set and go to step .
Since is differentiable and strictly convex, we can get the unique solution of (3.5a) which satisfies where and are the adjoint operators of and , and , respectively. It follows that the solution in (3.6a) can be directly obtained by the following explicit formulation: However, when the operator is ill-posed, the solution is unsuitable or unavailable. Hence we have to go back to (3.6a) and to employ some iteration strategy such as the Gauss-Seidel method to solve this equation. On the other hand, it is obvious that (3.5b) and (3.5c) can be looked at as the proximal mapping, so the solutions of the minimization problems (3.5b) and (3.5c) can be obviously written as
Notice that Algorithm 3.1 can be actually looked at as the split Bregman method . The based idea of this method is to introduce some intermediate variables so as to transform the original problem into some subproblems which are easily solved. The connection between the split Bregman method and the ADMM has been shown in [29, 49]. However, our algorithm considers the sum of three convex functions, which is more general than the related algorithms in [25, 49]. Furthermore, it must be noted that and are completely separated in (3.4a), so the two subproblems (3.5b) and (3.5c) are parallel. Therefore the convergence results of the ADMM can be applied here.
3.2. The Proximal Point Method
Though the ADMM in Algorithm 3.1 can effectively solve the original problem (3.1), we have to solve five subproblems. This actually makes this method suffer from a looser form as in [25, 45] so that it can badly affect its numerical computation efficiency. In this subsection, we propose a compacter form comparing to the ADMM. This formation called the PPM by using the relationship (2.9) in Remark 2.8 can reduce the original five subproblems of the ADMM in Algorithm 3.1 to solve three subproblems, thus it can improve computation cost of the ADMM. Now we have rewritten (3.5a)–(3.5e) with a little variation as the following form: with the first order optimality conditions given by If (3.11e) is replaced by it follows from (3.11a)–(3.11e) and Moreau-Yosida decomposition Theorem 31.5 in  that So we propose the following algorithm to solve (3.1).
Lemma 3.4. Set , and is a maximal monotone operator, then the operators and satisfy
Theorem 3.5. Assume that and are convex and proper. If and , here with for a continuous linear operator , then the sequence generated by Algorithm 3.3 converges to the limit point . Furthermore, the limit point is the solution of (1.8).
4. Some Applications in Image Restoration
In Section 4.1, we apply the ADMM and the above PPM to the image denoising problem. Here we also compare the proposed hybrid model with the ROF model and the LLT model. Then, based on the proposed hybrid model, we set the PPM as a basic method to solve image deblurring, image inpainting, and image denoising for the multiplicative noise in the last three subsections. For simplicity, we assume that the image region is squared with the size and set , , and as in . The usual scalar product can be denoted as for and for . The norm of is defined by and the norm of is defined by . If , we use to denote the first order forward difference operator with and use to denote the first order backward difference operator with for . Based on the first order difference operators, we can give the second order difference operator as follows: Using the same approach, we can define some other second order operators such as , , , and . Then we can give the first order and the second divergence operators as Furthermore, if we set and , it is easy to deduce that
Remark 4.1. The related examples in the following subsections are performed using Windows 7 and Matlab 2009(a) on a desktop with Intel Core i5 processor at 2.4 GHz and 4 GB memory. All of the parameters for related models are chosen by trial and empirically which can yield better restored images. On the other hand, we should notice that it is not very expensive when we use the ADMM and the PPM to get , but the total computational effort of one outer iteration requiring many inner steps can be very huge. In order to reduce the computational effort and keep fair comparison of these two methods, we so simplify the inner-outer iterative framework by performing only one-step in inner iteration. It is obvious that these sets are very efficient from the following numerical experiences.
4.1. Image Denoising for the Additive Noise
In this subsection, we consider to use the ADMM and the PPM to solve (1.7) for restoring the additive noisy image. If we set and , then the algorithms are proposed as follows.
Algorithm 4.2 (ADMM to solve (1.7)). (1) Choose the original and . Set , , and .
(2) Compute by where and .
(3) If the stop criterion is not satisfied, set and go to step .
Remark 4.3. For the first subproblem (4.6a), we can use the Gauss-Seidel method as shown in  to get the solution. However, in this paper, we use the following strategy: where some information of operators and related to is used. The formulas (4.6b) and (4.6c) of Algorithm 4.2 can be easily deduced from Furthermore, following form Theorem 3.2, we can also deduce that the sequence generated by Algorithm 4.2 converges to the solution of (1.7).
Algorithm 4.4 (PPM to solve (1.7)). (1) Choose the original , and . Set , , and .
(2) Compute by
(3) If the stop criterion is not satisfied, set and go to step .
Remark 4.6. For the above two algorithms, we can also set as a constant between 0 and 1. In fact, it is easy to find that the algorithms correspond to solving the ROF model or the LLT model, respectively, when or . At this time, the iteration strategy can be simplified as for these two models, respectively. Furthermore, when , these two algorithms correspond to solve the model which is the convex combination of the ROF model and the LLT model.
Example 4.7. In this example, we compare the ADMM with the PPM for solving the ROF model (1.4), the LLT model (1.6), and the hybrid model (1.7). The original images with three different sizes shown in Figure 1 are added to the Gaussian white noise with the standard deviation . Iterations were terminated when the stop conditions are met. It is easy to find the related results from Table 1 that the PPM is faster than the ADMM. Especially, the average CPU time of the PPM compared with that of the ADMM can save about 50% for the ROF model and the LLT model. It saves about 40% for the hybrid model.
Example 4.8. In this example, the noisy image is added to the Gaussian white noisy with the standard deviation . The algorithms will be stopped after 100 iterations. We compare the results generated by the ROF model, the LLT model, the convex combination of the ROF model and the hybrid model. As we can see from Table 2, the hybrid model get the lowest MSE and the highest SNR; these imply that the hybrid model can give the best restored image. On the other hand, it is easy to find that the ROF model makes staircasing effect appear and the LLT model leads to edge blurring. In fact, they are based on the fact that the restored model by the ROF model is piecewise constant on large areas and the LLT model as a higher model damps oscillations much faster in the region of edges. For the convex combined model and the hybrid model, they can efficiently suppress these two drawbacks. Furthermore, the hybrid model is more efficient than the convex combined model, because the hybrid model uses the edge detector function which can efficiently coordinate edge information. It should be noticed that here we use the Chambolle’s strategy  to solve the convex combined model so that it is slower than the strategy by using the PPM to solve the hybrid model. To focus on these facts, we present some zoomed-in local results and select a slice of the images which meets contours and the smooth regions shown in Figures 2 and 3.
4.2. Other Applications
In this subsection, we extend the hybrid model to other classes of restoration problems. As we can see in Section 4.1, the hybrid model has some advantages compared with the ROF model and the LLT model. Since the PPM is faster and more efficient than the ADMM, we only employ the PPM to solve the related image restoration model. Now we first consider the following iteration strategy: where . It is easy to find that the minimization problem in (4.11b) is coercive and strictly convex, so the subproblem (4.11b) has a unique solution. Based on [50, Lemma 2.4], we also deduce that the operator is -averaged nonexpansive.
4.2.1. Image Deblurring
Now we apply the hybrid model to the image deblurring problems with the following formula: Because of the compactness, the blurring operator is not continuously invertible, which implies that we cannot directly use the proximal method to solve minimization problem (4.14). However, based on the penalty method , we can transform the problem (4.14) to solve the following problem: The key of (4.15a)-(4.15b) is to separate the blurring operator from the original problem (4.14) so that we can avoid the ill-posed operator . Furthermore, it is obvious that the problem (4.15b) satisfie Theorem 4.9, so we can use the proximal method to solve it.
Example 4.11. In this experiment, we use the image Lena, which is blurred with a Gaussian kernel of “hsize = 3” and in which is added the Gaussian white noise with the standard deviation . The related images and SNRs are arranged in Figure 4. As we can see in Example 4.7, the proximal method tends to be stable before 100 iterations, so here we fix the algorithm to be stopped when the iteration attains 100. It is easy to find that the hybrid model can give a better image quality than the other two models. Especially, we also observe some staircasing effect in Figure 4(b) and edge blurring in Figure 4(c). However, all of the drawbacks can be efficiently suppressed in Figure 1(d).
(a) Original image
(b) ROF model
(c) LLT model
(d) Hybrid model
4.2.2. Image Inpainting
Here we consider to use the hybrid model for the image inpainting problem with the following form: where is the inpainting domain and satisfies If we introduce an auxiliary variable , based on the penalty method  again, then the solution of the problem (4.15a)-(4.15b) can be approximated by solving the following problem:
Example 4.12. In this example, we show the results of real image inpainting in Figures 6 and 7. We consider to use these three models to restore images shown in Figures 5(c) and 5(d), which are, respectively, contaminated by the mask image and the noise with the standard deviation . The parameters of these three models are chosen for getting better restored images and the algorithms will be stopped after 500 iterations (see also Table 3). It is obvious that these three models can efficiently restore the original deteriorated images. However, there are much more inappropriate information restored by the ROF model than the LLT model and the hybrid model. These advantages is based on the fact that the fourth-order linear diffusion (the LLT model and the hybrid model) damps oscillations much faster than second-order diffusion (the ROF model). Especially, these unsuitable effects can be easily seen from the local zooming images in Figures 6 and 7. Furthermore, the restored image by the hybrid model looks more natural than the LLT model.
(a) Original image
(b) Mask image
(c) Painting image
(d) Noisy image
4.2.3. Image Denoising for the Multiplicative Noise
Based on the hybrid model (1.7), the multiplicative noise model (1.5) can be naturally extended to the following minimization: It is obvious that the problem (4.19) can be approximated byFor the first subproblem (4.20a), its solution can be approximately obtained by using the Newton method to solve the following nonlinear equation:
Example 4.13. In this example, we restore the multiplicative noisy image. The noisy Lena image shown in Figure 7 is contaminated by the Gamma noise () with mean one, in which probability density function is given by where is an integer and is a Gamma function. Based on the iteration formula (4.20a)-(4.20b), we should notice that there are two interior iterations. For employing the Newton method to solve the problem (4.21), we set the stepsize and the Newton method will be stopped when . For solving the second subproblem (4.20b), we set the fixed iterations with 40. In fact for using the PPM to solve the ROF model, the LLT model and the hybrid model, the solutions of these models will tend to the steady. The outer iteration will be stopped after 200 iterations. The related restored images are shown in Figure 8, and it is easy to find that the hybrid model gives a better restored image than the two other models.
5. Concluding Remarks
In this paper, based on the edge detector function, we proposed a hybrid model to overcome some drawbacks of the ROF model and the LLT model. Following the augmented Lagrangian method, we can employ the ADMM of multipliers to solve this hybrid model. In this paper, we mainly proposed the PPM to solve this model due to the fact that the PPM unnecessarily solves a PDE compared with the ADMM so that it is more effective than the ADMM. The convergence of the proposed method was also given. However, the convergence rate of the proposed method is only , so our future work is to improve our method with convergence rate based on the same strategies in [51, 52].
A. Proof of Theorem 3.2
Proof. Following the assumption, we can find that the saddle point of (3.9) satisfies
with and . If we set and , then (A.1) can be written as Then the above equations (A.2a)–(A.2e) can be rewritten as a compact form
which implies that is the solution of the minimization problem (1.8).
Set denote the related errors. Then subtracting (3.6a)–(3.6e) by (A.2a)–(A.2e) and using the similar strategy as in , we successively obtain where and . Summing the above five equations, we get So it is easy to deduce that Following the fact that is strictly convex, and , we can get Combined with (A.7), it follows that Then we deduce that So we have Hence the proof is complete.
B. Proof of Lemma 3.4
Proof. Let be solutions of then we can obtain which implies that the assertion holds.
C. Proof of Theorem 3.5
Proof. Setting and , then operators and are maximal monotone. Following from (3.11a)–(3.11e) and the definition of the Yosida approximation, (3.13) can be equivalently written as Then, by (C.1a)–(C.1c) and Lemma 3.4, we can deduce that
Combine (C.2) and (C.3) with
which implies that
Since and , by (C.3) we eventually get
as long as and
Furthermore, by (C.4) we also get
Based on (C.1a) and (C.1b), using the definitions of resolvent and the Yosida approximation, it is easy to find that
Then, using Corollary 4 in page 199 of , we can deduce that and .
On the other hand, if is the limit of the sequence generated by Algorithm 3.3, it is easy to find that (3.13) can be rewritten as as . Using the relationship (2.9) in Remark 2.8, (C.1a)–(C.1c) can be then rewritten as Following from (C.10a) and (C.10b), we can find that Submitting (C.10a)–(C.10c) into (C.11), we have which implies that is the solution of (1.8).
D. Proof of Theorem 4.9
Proof. From the assumption, the functional exists in at least a minimized point denoted by , which implies that satisfies
That is to say that is a fixed point of the operator .
Now we show that the sequence generated by (4.11a)-(4.11b) converges to this fixed point . In fact, based on the nonexpansive operators and , for the sequence we can get This implies that the positive sequence is monotonically decreasing. So we can find a subsequence from converging to a limit point as . Form the continuity of and