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 L 1 -norm terms in the proposed model make it di ﬃ cultly solved by using some classically numerical methods directly, we ﬁrst employ the alternating direction method of multipliers (cid:3) ADMM (cid:4) to solve a general form of the proposed model. Then, based on the ADMM and the Moreau-Yosida decomposition theory, a more e ﬃ cient method called the proximal point method (cid:3) PPM (cid:4) is proposed and the convergence of the proposed method is proved. Some numerical results demonstrate the viability and e ﬃ ciency of the proposed model and methods.


Introduction
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 4 to overcome it.The general form of regularization methods consists in minimizing an energy functional of the following form: in Banach space X, where λ is the regularization parameter, R, is a regularization term, u 0 is the observed image, and u is the image to be restored.The development of the energy functional 1.1 actually profits from the ROF model 5 which is of the following form: where z is the logarithmic transformation of u and |z| BV keeps the total variation property related to |e z | BV and α > 0.Here u satisfies f uη for the noise η.
However, as we all know the | • | BV 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 |u| BV 2 : Ω |∇ 2 u|dx for the Hessian operator ∇ 2 .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 BV and BV 2 was proposed by Lysaker et al. 13 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 g : g u 0 : 1/ 1 |∇ u 0 * G σ x | with G σ x 1/2πσ 2 exp − x 2 /2σ 2 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 16 .Obviously, | 1 − g u| BV tends to be predominant where edges most likely appear and |gu| BV 2 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 22 , the multigrid method 23 , operator splitting method 24-26 , the inverse scale method 27 , and so forth.However, different from the ROF model 1.2 and the LLT model 1.6 , the model 1.7 includes two L 1 -norm terms which make it solved more difficultly.More generally, the model 1.7 can be fell into the following framework: where h 1 , h 2 : X → R are proper, convex, and lower semicontinuous l.s.c., Λ 1 and Λ 2 are bounded linear operators, and θ is a parameter.A specific form of 1.8 was considered by Afonso and Bioucas-Dias in 28 where a Bregman iterative method was proposed to solve the model with the combination of the L 1 -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 35 .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.

Notations and Definitions
Let us describe some notations and definitions used in this paper.For simplifying, we use X R n .Usually, we can set n 2 for the gray-scale images.The related contents can be referred to 1, 27, 36-43 .

Definition 2.1. The operator
where x ∈ X and τ > 0.
Definition 2.7.The shrinkage operator S τ • : X → X is defined by where we use the convention 0/0 0.

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.

The Alternating Direction Method of Multipliers (ADMM)
We now consider the following constrained problem: where ζ i is the Lagrangian multiplier and μ i is the penalty parameter for i 1, 2. Then we can use the following augmented Lagrangian method ALM : with choosing the original values ζ 0 1 and ζ 0 2 to solve 3.2 .If we set d n i ζ n i /μ i for i 1, 2 and omit the terms which are independent of u n , v n , z n in 3.3a , the above strategy 3.3a -3.3c can be written as for the original values d 0 1 and d 0 2 .Then we can use the following ADMM to solve 3.4a -3.4c .
Algorithm 3.1 ADMM for solving 3.4a -3.4c . 1 Choose the original values: v 0 , z 0 , d 0 1 , and 3 If the stop criterion is not satisfied, set n : n 1 and go to step 2 .Since h u is differentiable and strictly convex, we can get the unique solution of 3.5a which satisfies where Λ * 1 and Λ * 2 are the adjoint operators of Λ 1 and Λ 2 , w n 1 ∈ ∂h 1 v n and w n 2 ∈ ∂h 2 z n , respectively.It follows that the solution u n in 3.6a can be directly obtained by the following explicit formulation:

3.7
However, when the operator M 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 Then u * is the solution of the minimization problem 1.8 .Furthermore, the sequence { u n , v n , z n } generated by Algorithm 3.1 converges to u * , v * , z * .
Notice that Algorithm 3.1 can be actually looked at as the split Bregman method 25 .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 v and z 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.

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:

3.10e
Abstract and Applied Analysis 9 with the first order optimality conditions given by it follows from 3.11a -3.11e and Moreau-Yosida decomposition Theorem 31.5 in 35 that

3.13
So we propose the following algorithm to solve 3.1 .
3.14 Theorem 3.5.Assume that h 1 x and h 2 u are convex and proper.If , here • : max{ Kx L 2 : x ∈ X with x L 2 ≤ 1} for a continuous linear operator K, then the sequence { u n , d n 1 , d n 2 } generated by Algorithm 3.3 converges to the limit point d 1 , d 2 , u .Furthermore, the limit point u is the solution of 1.8 .

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 M × M and set S R M×M , T S × S, and Z T × T as in 17 .The usual scalar product can be denoted as p 1 , p 2 T : M i 1 M j 1 p 1 i,j p 2 i,j for p 1 , p 2 ∈ T and p, q Z M i,j 1 p 11 i,j q 11 i,j p 12 i,j q 12 i,j p 21 i,j q 21 i,j p 22 i,j q 22 i,j for p, q ∈ Z.The L 1 norm of p p 1 , p 2 ∈ T is defined by |p| 1 p 2 1 p 2 2 and the 1 norm of q q 1 , q 2 ; q 3 , q 4 ∈ Z is defined by If u ∈ S, we use ∇ ∇ x , ∇ y ∈ T to denote the first order forward difference operator with x , ∇ − y to denote the first order backward difference operator with for j N, 4.2 for i, j 1, . . ., N. 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 ∇ y ∇ − y u i,j .Then we can give the first order and the second divergence operators as

4.4
Furthermore, if we set p ∈ T and q ∈ Z, it is easy to deduce that div p 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 u n , 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.

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 Λ 1 ∇ and Λ 2 ∇ 2 , then the algorithms are proposed as follows. where 3 If the stop criterion is not satisfied, set n : n 1 and go to step 2 .
Remark 4.3.For the first subproblem 4.6a , we can use the Gauss-Seidel method as shown in 25 to get the solution.However, in this paper, we use the following strategy: where some information of operators Δ and Δ 2 related to u is used.The formulas 4.6b and 4.6c of Algorithm 4.2 can be easily deduced from

4.8
Furthermore, following form Theorem 3.2, we can also deduce that the sequence {u n } generated by Algorithm 4.2 converges to the solution of 1.7 . 2 Compute 4.9 3 If the stop criterion is not satisfied, set n : n 1 and go to step 2 .
For Algorithm 4.4, based on the relations in 4.5 and Theorem 3.5, we have the following result.Theorem 4.5.If μ 1 ∈ 0, 1/8θ and μ 2 ∈ 0, 1/64θ , then the sequence {u n } generated by Algorithm 4.4 converges to the solution of 1.7 .Remark 4.6.For the above two algorithms, we can also set g 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 g 0 or g 1.At this time, the iteration strategy can be simplified as for these two models, respectively.Furthermore, when g ∈ 0, 1 , these two algorithms correspond to solve the model which is the convex combination of the ROF model and the LLT model.  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 σ 12.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 17 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.

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 h z ∈ C 1 Ω .It is easy to find that the minimization problem D u 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 G is 1/2 -averaged nonexpansive.Theorem 4.9.Assume that the functional 12 is convex, coercive, and bounded below, then the sequence { u n , z n } generated by 4.11a -4.11b converges to a point u * , z * .Furthermore, when u * z * , u * is the solution of Remark 4.10.It should be noticed that the following three models satisfy the conditions of Theorem 4.9, so we can use the strategy 4.11a -4.11b to solve these models.

Image Deblurring
Now we apply the hybrid model to the image deblurring problems with the following formula: Because of the compactness, the blurring operator K 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 48 , 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 K from the original problem 4.14 so that we can avoid the ill-posed operator K. Furthermore, it is obvious that the problem 4.15b satisfies 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 σ 0.02.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 .

Image Inpainting
Here we consider to use the hybrid model for the image inpainting problem with the following form: where D is the inpainting domain and λ satisfies 4.17  If we introduce an auxiliary variable z, based on the penalty method 48 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 σ 12.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.

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 by For the first subproblem 4.20a , its solution u n 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 L 33 with mean one, in which probability density function p u is given by  where L 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 t 1 and the Newton method will be stopped when u n 1 − u n 2 / u n 2 ≤ 1.0 × 10 −5 .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.

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 O 1/k , so our future work is to improve our method with convergence rate O 1/k 2 based on the same strategies in 51, 52 .

Appendices
A. Proof of Theorem 3.2 Proof.Following the assumption, we can find that the saddle point u * , v * , z * , w * 1 , w * 2 of 3.9 satisfies Then the above equations A.2a -A.2e can be rewritten as a compact form which implies that u * 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 24 , we successively obtain where d n 1e ∈ 1/μ 1 ∂h 1 v n 1e and d n 2e ∈ 1/μ 2 ∂h 2 z n 1e .Summing the above five equations, we get A.6 So it is easy to deduce that Following the fact that h x is strictly convex,

C. Proof of Theorem 3.5
Proof.Setting A 1 x : ∂h 1 x and A 2 x : ∂h 2 x , then operators A 1 and A 2 are maximal monotone.Following from 3.11a -3.11e and the definition of the Yosida approximation, 3.13 can be equivalently written as C.1c  Then, by C.1a -C.1c and Lemma 3.4, we can deduce that

D. Proof of Theorem 4.9
Proof.From the assumption, the functional M u, z exists in at least a minimized point denoted by u, z , which implies that u, z satisfies z F u , u G z .

D.1
That is to say that u is a fixed point of the operator G • F.

2 2 ≤ 8 p 2 11 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.

Figure 1 : 7 . 4 . 7 .
Figure 1: The original images and the noisy images with three different sizes in Example 4.7.

Figure 2 :
Figure 2: The original and the noisy image in Example 4.8.

Figure 3 :
Figure 3: The related restored images in Example 4.8.

Figure 5 :
Figure 5: The related images in Example 4.12.

Figure 6 :
Figure 6: The related restored images corresponding to the painting image in Example 4.12.

Figure 7 :
Figure 7: The related restored images corresponding to the noisy image in Example 4.12.

1 2L 2 , C. 4 Since μ 1 ≤ 1 /θ Λ 1 2 and μ 2 ≤ 1 /θ Λ 2 2 , 1 L 2 , d 2 − d n 2 L 2 ≤ d 2 −
by C.3 we eventually getd 1 − d n 1 L 2 ≤ d 1 − d n−1 1.2 where |u| BV : Ω |∇u|dx.In the case without confusion, for simplification we omit the open set Ω with Lipschitz boundary for • L 2 Ω and |•| L 1 Ω .Due to edge-preserving property of the term |u| BV , this model has been extended to other sorts of image processing problems such as for image inpainting 2, 7 , where K is a blurring operator and D 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 8 , synthetic aperture radar SAR 9 , and medical ultrasonic images 10 .One of the models based on BV was proposed byHuanget al. HNW 11 with the form min z α Ω z fe −z dx |z| BV , 1.5 1 for all y 1 ∈ A x 1 and y 2 ∈ A x 2 , and A is maximal if it is not strictly contained in any other monotone operator on X.The projection operator P B τ • : X → X onto the closed disc B τ : {x ∈ X : |x| L 1 ≤ 1/τ} is defined by If the stop criterion is not satisfied, set n : n 1 and go to step 2 .Set x 1 , x 2 ∈ X and A is a maximal monotone operator, then the operators H cA and J cA satisfy 3.3 PPM for solving 3.1 . 1 Choose the original values: d

Table 1 :
The related results in Example 4.7.

Table 2 :
The related data in Example 4.8 here R.P. is the regularization parameter .

Table 3 :
The related results in Example 4.12.