A wavelet frame coefficient total variational model for image restoration
Wei Wang Xiang-Gen Xia Shengli Zhang Chuanjiang He


Critical role of water in defect aggregation and chemical degradation of perovskite solar cells
    [
    Received: date / Accepted: date
===============================================================================================



[1]College of Information Engineering, Shenzhen
	University, Shenzhen, China. (mailto:aiyulunhui@sina.comaiyulunhui@sina.com, mailto:xxia@ee.udel.eduxxia@ee.udel.edu, mailto:zsl@szu.edu.cnzsl@szu.edu.cn
mailto:cjhe@cqu.edu.cncjhe@cqu.edu.cn) 
[2]Department of Electrical
	and Computer Engineering, University of Delaware, Newark, DE 19716, USA.
[3]College of Mathematics and Statistics,Chongqing University,Chongqing 401331, China.
	
This research is partially supported by the research grants from the Chinese NSF project (61372078), National Natural Science Foundation of China (No. 61701320, No. 61561019), Guangdong NSF project (2014A030313549), Key project of Guangdong (2016KZDXM006) and the Shenzhen NSF project (JCYJ20160226192223251).


definition
LemmaLemma
RemarkRemark 
TheoremTheorem
DefinitionDefinition
*keywordskeywords






   In this paper, we propose a vector total variation (VTV) of feature image model for image restoration.  The  VTV imposes different smoothing powers on different features (e.g. edges and cartoons) based on choosing various regularization parameters. Thus, the model can simultaneously preserve edges and remove noises. Next, the existence of  solution for the model is proved and the split Bregman algorithm is used to solve the model. At last, we use the wavelet filter banks to explicitly define the feature operator and 
   present some experimental results to show its  advantage over the related methods  in both quality and efficiency.




  Image restoration, wavelet, vectore total variation, variational method，






§ INTRODUCTION

Image restoration including image denoising, deblurring, inpainting etc. is a procedure of improving the quality of a given image that is degraded in various ways during the process of acquisition and commutation. Since many advanced applications in computer vision depend heavily on the input of high quality images, image restoration becomes an indispensable and preprocessing step of these applications. Therefore, image restoration is a basic but very important area in image processing and analysis. 
Image restoration can be modeled as a linear inverse problem:

    f = Au + η

where f  is the observed images, A  is a  linear operator and η  represents the additive noise whose type depends on its probability density function (PDF). The goal of image restoration is to find the unknown true image u  from the observed image  f. The problem is usually an ill-posed inverse problem and often solved by imposing a prior regularization  assumption. Among all regularization-based methods for image restoration, variational based <cit.><cit.><cit.><cit.><cit.><cit.><cit.> and wavelet frame based <cit.><cit.><cit.><cit.><cit.><cit.><cit.><cit.><cit.><cit.><cit.><cit.><cit.><cit.><cit.> methods have attracted a lot of attention in the past.

The common assumption of the regularization based methods is that images can be sparsely approximated in some transformed domains. Such transforms can be gradient operator, wavelet frame transform, Fourier transform, Gabor transform  etc. In order to utilize the sparsity, one solves (1) by finding a sparse solution in the corresponding transformed domain. Typically, the l_1  norm is used as a penalty of the sparsity. 
One of the well-known variational approaches is the Rudin-Osher-Fatemi (ROF) <cit.> model, which imposes the total variation (TV) regularization on the image  u:

    inf_u ∈ BV(Ω ){λ∫_Ω|∇ u|  + ∫_Ω(u - f)^2 dx}.

Here, Ω is the image domian, BV(Ω) is the bounded variational space,  ∫_Ω|∇ u|  and ∫_Ω(u - f)^2 dx are the TV (regularization) term and fidelity (fitting) term, respectively. The ROF model performs well for removing noise while preserving edges. However, it tends to produce piecewise constant results (called staircase effect in the literature). After the ROF model was proposed, the TV regularization has been extended to many other image restoration applications (such as deblurring, inpainting, superresolution, etc.) and has been modified in a variety of ways to improve its performance <cit.><cit.><cit.><cit.>.To solve the ROF model, many algorithms have been proposed <cit.><cit.><cit.><cit.><cit.> <cit.><cit.>and the split Bregman method <cit.> is the one of the most widely used ones . 

Wavelet frames represent images as a summarization of smooth components (i.e. cartoon), and local features (i.e. singularities). In wavelet frame domain, smooth image components are the coefficient images obtained from low-pass filters, while local features are those obtained from high-pass filters. Thus, we can impose different strengths of regularization on smooth components and local features separately when restoring images in the wavelet frame transform domain. There exist plenty of wavelet frame based image restoration models in the literature, such as, the synthesis based approach <cit.><cit.><cit.><cit.>, the analysis based approach <cit.><cit.><cit.><cit.>, and the balanced approach <cit.><cit.><cit.>. A typical analysis based model of l_1  norm regularization has a similar form with the ROF model:

    min_u λWu_1 + 1/2Au - f_2^2,

where W  represents the wavelet frame transform. Since the l_1  norm regularizations usually exist in the wavelet frame based models, the split Bregman algorithm is also a widely used method to solve them.

In this paper, we propose a vector total variation (VTV) of feature image model  for image restoration where the VTV is used as the regularization. The vector images are generated by mapping the original image into the feature spaces which represent different features of the image such as edges and cartoons. By choosing different regularization parameters, we can impose different smoothing powers on different features. Thus, the model can simultaneously preserve edges and remove noises. Next, the existence of  solution for the model is proved and the split Bregman algorithm is used to solve the model. At last, we use the wavelet filter banks to explicitly define the feature operator and  present some experimental results to show its  advantage over the related methods  in both quality and efficiency.

The rest of the paper is organized as follows. In Section 2.1, we describe the proposed model, prove the existence of its solution and use the convolution to define the feature operators. In Section 2.2, we use the split Bregamn algorithm to solve the model. In Section 3, we present some experimental results to show the advantage of our model. In Section 4, we conclude this paper.



§ THE WAVELET  COEFFICIENT TOTAL VARIATIONAL MODEL




 §.§ Proposed model

Before describing the model, we briefly introduce some mathematical definitions related to our paper.

Let Ω be a bounded open subset of R^2, f∈ L^2(Ω,R),  F=(F_1,F_2,...,F_m) be an operator such that Ff=(F_1f,F_2f,...,F_mf)∈ L^2(Ω,R^m) and F_if∈ L^2(Ω,R). In this paper, we call F_if a feature image of f and F a feature operator.

	The inner product in  L^2(Ω,R^m) is defined as 
	
    <u,v>_L^2(Ω,R^m)=∫_Ω∑_i=1^mu_iv_i,

	where u=(u_1,u_2,...,u_m),v=(v_1,v_2,...,v_m)∈ L^1(Ω,R^m).
	

	BV(Ω) is the subspace of L^1(Ω,R) such that the following quantity is finite:
	
    |f|_TV=sup_η∈ K<f,divη>_L^2(Ω,R),

	where divη=∂η_1/∂ x+∂η_2/∂ y, K={η=(η_1,η_2)∈ C_0^1(Ω,R^2):η_L^∞(Ω,R^2)≤ 1}, C_0^1(Ω,R^2) is the space of continuously differentiable functions with compact support in Ω, η_L^∞(Ω,R^2)=sup_x∈Ω√(η_1^2+η_2^2) and the quantity |f|_TV is called the total variation of f.
	


  The vector total variation of function g=(g_1,g_2,...,g_m)∈ L^1(Ω,R^m) is defined as 
  
    |g|_VTV=sup_ξ∈ P<g,divξ>_L^2(Ω,R^m),

  where ξ=(ξ_1,ξ_2,...,ξ_m), divξ=(divξ_1,divξ_2,...,divξ_m), P={ξ∈ C_0^1(Ω,R^2× m):ξ(x)_L^∞(Ω,R^2× m)≤ 1} and ξ(x)_L^∞(Ω,R^2× m)=max_i=1,2...,mξ_i(x)_L^∞(Ω,R^2).
	
For g smooth enough, we have <cit.><cit.>

    |g|_VTV= ∑_i=1^m|g_i|_TV= ∑_i=1^m∫_Ω|∇ g_i|.

For more details about the vector total variation, refer to <cit.><cit.><cit.><cit.><cit.>.

In this paper, we propose the following model for image recovery:

    inf_u ∈ BV(Ω ){ E(u) = | λ⃗Fu|_VTV +1/2Au - f_L^2(Ω,R)^2},

where A is a linear operator, f is an observed image, F is a feature operator, λ⃗=(λ_1,λ_2,...,λ_m)∈ (R^+)^m is the regularization parameter vector and λ⃗Fu=(λ_1F_1u,λ_2F_2u,...,λ_mF_mu) represents the feature images. The ideal case is that different F_mu can represent different features of the image u, such as edges, cartoons, etc. and thus we can impose different smooth strengths on them by choosing different λ_i. Clearly, if m=1 and F=I, then the proposed model is reduced to the ROF model.
   
In the following, we show the existencce of the solution for problem (<ref>).

	If F_i is linear, then the adjoint operator F^* of F is
	
    F^*g=∑_i=1^mF_i^*g_i,

	where  g=(g_1,g_2,...,g_m)∈ L^2(Ω,R^m) and F_i^* is the adjoint operator of F_i.


	Let f∈  L^2(Ω,R) and g ∈ L^2(Ω,R^m), then 
	
    <Ff,g>=   ∑_i=1^m<F_if,g_i>
    
    	=   ∑_i=1^m<f,F_i^*g_i>
    
    	=   <f,∑_i=1^mF_i^*g_i>

Thus, 	F^*g=∑_i=1^mF_i^*g_i.



	If F is linear and bounded and for any ξ∈ C_0^1(Ω,R^2× m), the operator div and F^* are commutative, i.e., 
	F^*divξ=div F^*ξ,
	then for any f∈ BV(Ω), we have
	
    |Ff|_VTV≤F|f|_TV.




    |Ff|_VTV=   sup_ξ∈ P<Ff,divξ>
    
    =   sup_ξ∈ P<f,F^*divξ>
    
    =   sup_ξ∈ PF^*ξ_L^∞(Ω,R^2)<f,div(F^*ξ/F^*ξ_L^∞(Ω,R^2))>
    ≤   sup_ξ∈ PF^*ξ_L^∞(Ω,R^2× m)<f,div(F^*ξ/F^*ξ_L^∞(Ω,R^2))>
    ≤   F^*|f|_TV
    
    =   F|f|_TV



	If F is linear and bounded,  invertible, and for any η∈ C_0^1(Ω,R^2), the operator div and (F^-1)^* are commutative, i.e., 
	(F^-1)^*divη=div (F^-1)^*η,
	then for any f∈ BV(Ω), we have
	
    |f|_TV≤F^-1|Ff|_VTV.




	
    |f|_TV=   sup_η∈ K<F^-1Ff,divη>
    
    	=   sup_η∈ K<Ff,(F^-1)^*divη>
    
    	=   sup_η∈ K(F^-1)^*η_L^∞(Ω,R^2× m)<Ff,div((F^-1)^*η/(F^-1)^*η_L^∞(Ω,R^2× m))>
    ≤   sup_η∈ K(F^-1)^*η_L^∞(Ω,R^2)<Ff,div((F^-1)^*η/(F^-1)^*η_L^∞(Ω,R^2× m))>
    ≤   (F^-1)^*sup_η∈ K<Ff,div((F^-1)^*η/(F^-1)^*η_L^∞(Ω,R^2× m))>
    ≤   (F^-1)|Ff|_VTV



	From the proof, we can see the condition we need is that there exists F_L^-1: L^2(Ω,R^m)→ L^2(Ω,R) such that F_L^-1F=I and (F_L^-1)^*divη=div(F_L^-1)^*η for any η∈ C_0^1(Ω,R^2).


    If F is  linear and bounded, and invertible, then  for any ξ∈ C_0^1(Ω,R^2× m) and η∈ C_0^1(Ω,R^2), the following two conditions are equivalent:
    
    	
  (a) F^*divξ=div F^*ξ,
    	
  (b) (F^-1)^*divη=div (F^-1)^*η.
    
 By 
	(a) ⇒ (b):
	Let ξ=(F^-1)^*η in (a), then
	
    F^*div(F^-1)^*η=divF^*(F^-1)^*η=div(F^-1F)^*η=divη.

	Multiplying (F^-1)^* on both sides of the above equation, we obtain
	
    div(F^-1)^*η=(F^-1)^*divη,

	
	(b) ⇒ (a): Let η=F^*ξ in (b), then we have
		
    (F^-1)^*divF^*ξ=div (F^-1)^*F^*ξ=divξ.

		Therefore, 
    divF^*ξ=F^*divξ.

 


	By 1=F^-1F≤F^-1F, we only have F^-1≥F^-1. If  F^-1=F^-1, then combining Lemma 2 and Lemma 3, we can obtain |Ff|_VTV=F|f|_TV.


	If the conditions of Remark 1 are satisfied, λ_min=min_i=1,2,...mλ_i>0 and  A⊗ 1=∫_Ω A(x)dx 0, then Problem (<ref>) admits a solution.


	Let u_n be a minimizing sequence of Problem (<ref>), then there exists a constant M>0 such that E(u_n)≤ M, i.e.

    | λ⃗Fu|_VTV + 1/2Au - f_L^2(Ω,R)^2≤ M.

Therefore, 

    | λ⃗Fu|_VTV≤ M


    Au - f_L^2(Ω,R)^2≤ 2M

From equation (<ref>), by Lemma 3, we have 

    |u|_TV   ≤F^-1|Fu|_VTV
       ≤1/λ_minF^-1|λ⃗Fu|_VTV
       ≤M/λ_minF^-1

As proved in in <cit.><cit.>, from equation (<ref>), we can obtain

    u_L^2(Ω,R)≤ M


Combining equations (<ref>) and (<ref>), we have that the sequence {u_n} is bounded in BV(Ω).  Thus, there exists u^*∈ BV(Ω) such that u_n  converges to u^*  weakly in L^2(Ω )  and strongly in  L^1(Ω ). Since the operator F  is linear and bounded, we have that Fu_n also converges to Fu^*  weakly in L^2(Ω )  and strongly in  L^1(Ω ). Then a standard process can show that u^*  is a minimizer of  E(u).


	If there exists a F_k such that F_k is invertible and F_kdivη=div F_kη for any η∈ C_0^1(Ω,R^2), λ_min=min_i=1,2,...mλ_i> 0 and  A⊗ 1=∫_Ω A(x)dx 0, then Problem (<ref>) admits a solution.


	Let u_n be a minimizing sequence of Problem (<ref>), then there exists a constant M>0 such that E(u_n)≤ M, i.e.
	
    | λ⃗Fu|_VTV + 1/2Au - f_L^2(Ω,R)^2≤ M.

	Therefore, 
	
    | λ⃗Fu|_VTV=∑_i^m|λ_iF_iu|_TV≤ M

	and so
		
    |F_ku|_TV≤M/λ_k

	
	Since F_k is invertible and F_kdivη=div F_kη for any η∈ C_0^1(Ω,R^2), like the proof in Lemma 2, we can obtain
	
    |u|_TV≤F_k^-1|F_ku|_TV≤M/λ_kF_k^-1

	
	The left part of the proof is the same as that in Theorem 1 and so is omitted.	


To assure the existence of the solution, we can intentionally design a F_k such that F_k satisfies the conditions in Theorem 2. The simplest F_k satisfying the conditions is F_k=I.




In the following, we give the definition of the feature operator F. In  signal processing, a filter of convolution is usually used to extract some feature of a signal. In this paper, we also use the convolution to define F_i.

Let {K_i(x)} be a family of functions with a compact support E=[-r,r]× [-r,r] such that K_i(x)=0 for x∈ R^2-E and K_i(x)∈ L^∞(R^2,R). For any f∈  L^2(Ω,R), we define F_if as

    F_if(x)=f⊗ K_i=∫_Ω f(y)K_i(x-y)dy.

Since 

    |∫_Ω f(y)K_i(x-y)dy|^2   ≤∫_Ω f^2(y)dy∫_Ω K_i^2(x-y)dy
       ≤K_i^2_L^∞(R^2,R)f^2_L^2(Ω,R)

we have F_if∈ L^2(Ω,R).




By the Fubini's theorem,  we have

    F_i^*f(x)=f⊗K̂_i=∫_Ω f(y)K_i(y-x)dy,

where K̂_i(x)=K_i(-x).

For any η=(η_1,η_2)∈ C_0^1(Ω,R^2), we have

    F_i^*divη   =∫_Ωdivη(y)K_i(y-x)dy
       =-∫_Ωη(y)∇_y (K_i(y-x))dy
       =∫_Ωη(y)∇_x (K_i(y-x))dy
       =∫_Ωη_1(y)∂K_i(y - x)/∂x_1+ η_2(y)∂K_i(y - x)/∂x_2dy
       =div(∫_Ωη_1(y)K_i(y-x)dx,∫_Ωη_2(y)K_i(y-x)dx)
       =div(F_i^*η_1,F_i^*η_2)
       =divF_i^*η

where we assume x=(x_1,x_2).
Combining the above equation with Lemma 1, we thus have 

    F^*divη=divF^*η.



 §.§ Algorithm

In this section, we utilize the split Bregman algorithm <cit.> to solve problem (<ref>). Before describing the algorithm, we need to give the discrete version of Problem (<ref>).
 
Let f be a  matrix, such as an image of size s× t, K_i a matrix convolution filter of size (2r+1)× (2r+1). Then the discrete convolution of f and K_i is defined by

    F_if(k,l)=(f⊗ K_i)(k,l)
       =∑_p=k-r,q=l-r^p=k+r,q=l+rf(k,l)K_i(k-p,k-q)
       =∑_p,q=-r^rf(k-p,l-q)K_i(p,q)

where we assume the center of K_i is the  origin of coordinates and extend f circularly. 
By equation (<ref>), we have 

    F_i^*f(k,l)=(f⊗K̂_i)(k,l)
       =∑_p=k-r,q=l-r^p=k+r,q=l+rf(k,l)K_i(k-p,l-q)
       =∑_p,q=-r^rf(k-p,l-q)K_i(-p,-q)


Let g=(g_1,g_2,...,g_r), g_i be a matrix of size s× t. Define the p-norm of g as

    g_1   =∑_i=1^rg_i_1
    g_2   =√(∑_i=1^rg_i_2^2)

Then by equation (<ref>),  the discrete version of Problem (<ref>) is 


    min_u{E(u)=∑_i=1^mλ_i∇  F_iu_1+1/2Au - f_2^2}


 Let  d = ∇ (Fu), i.e. d_i=∇ (F_iu), d^0 = b^0 = 0, then Problem (<ref>) is equivalent to the following iterations:

for j=0,1,2,…

    u^j + 1 = min1/2Au - f_2^2 + ∑_i=1^mγ_i/2∇ (F_iu) - d^j_i + b^j_i_2^2
       d^j + 1 =min∑_i=1^m(λ_i d_i _1 + γ_i /2d_i - ∇ ( F_iu^j + 1) - b_i^j_2^2)
       b_i^j + 1 = b_i^j + (∇ ( F_iu^j + 1) - d_i^j + 1)


The KKT condition for problem (<ref>) is 

    [ A^*(Au - f) + ∑_i=1^mγ_iF_i^*∇ ^* (∇ F_iu - d_i^j + b_i^j)
    = 0 ]
 Therefore, the solution for Problem (<ref>) is 

    u^j + 1 = FFT^ - 1(FFT(∑_i=1^mγ_iF_i^*∇ ^* (d_i^j - b_i^j) + A^*f)/FFT(A^*A +  ∑_i=1^mγ_iF_i^*∇ ^*∇ F_i))


The solution for problem (<ref>) can be explicitly solved:

    d_i^j + 1 = TH(∇ ( F_iu^j + 1) + b_i^j,λ_i/γ_i )

where  TH(x,T) =  sgn (x)max (| x | - T,0).

At last, the algorithm is summarized in Algorithm <ref>.


	To implement Algorithm 1, if we don't consider the existence of solution, we only need F_i^* but don't need F^-1 or F_i^-1.

The proof of the convergence of the split Bregman algorithm was given in <cit.>. The  convergence of our algorithm can be proved accordingly and so is not listed here.




§ EXPERIMENTAL RESULTS

In this section, we conduct some numerical experiments on image denoising and image deblurring using Algorithm <ref>. To implement Algorithm 1, we need to construct  K_i explicitly. In our experiments, the piecewise linear B-spline wavelet frame is used for constructing K_i.
The filter banks  of the B-spline wavelet frame are

    [h_1( - 1),h_1(0),h_1(1)] = 1/4[1,2,1]
  
    [h_2( - 1),h_2(0),h_2(1)] = √(2)/4[ 1,0,-1]
 
    [h_3( - 1),h_3(0),h_3(1)] = 1/4[ - 1,2, - 1].

Then the K_i used in our experiment are
K_1=h_1^† h_1 and 

    K_3(i-1)+j=h_i^† h_j,    i,j=1,2,3.

Clearly, the feature image generated by K_1 have more even region than those by the other K_i. So in the experiments, we set a larger λ_1 and γ_1 to impose a more smooth effect on it.

By the unitary extension principle (UEP) in <cit.>,  we have that 

    F^*F=I.

 Thus, F_L^-1=F^*.
By Lemma 1 and (<ref>), we also have  

    ∑_i=1^mF_i^*F_i=1


If we set the parameter vector γ  as γ_1>γ_2=γ_3=...=γ_m, then by equation (<ref>), we have

    ∑_i=1^mγ_iF_i^*∇ ^*∇ F_i   =(γ_1-γ_2)F_1^*F_1∇^*∇+γ_2∇ ^*∇∑_i=1^mF_i^* F_i
       =(γ_1-γ_2)F_1^*F_1∇^*∇+γ_2∇ ^*∇,

where we can interchange these operators since they are all convolutions. Thus equation (<ref>) in Algorithm 1 can be reduced to

    u^j + 1 = FFT^ - 1(FFT(∑_i=1^mλ_iF_i^*∇ ^* (d_i^j - b_i^j) + A^*f)/FFT(A^*A +  (γ_1-γ_2)F_1^*F_1∇^*∇+γ_2∇ ^*∇))




In Figure <ref>, we show the Fast Fourier transform (FFT) of ∇ ^*∇ and
F_1^*F_1, from which we can see that  the effect of the FFT of F_1^*F_1 is opposite to that of ∇ ^*∇. Since the effect of  ∇ ^*∇ is smoothing the image, the term F_1^*F_1 may slow down the convergence speed and cause artifacts in the denoised images. Thus in the numerical simulations, we neglect the term (γ_1-γ_2)F_1^*F_1∇^*∇  for computational efficiency. Equation (<ref>) is thus changed into

    u^j + 1 = FFT^ - 1(FFT(γ∑_i=1^mλ_iF_i^*∇ ^* (d_i^j - b_i^j) + A^*f)/FFT(A^*A + γ _1∇ ^*∇ ))

Equation (<ref>) can be seen as an analogy of the original split Bregman algorithm who solves u^j+1 by

    u^j + 1 = FFT^ - 1(FFT(λ∇ ^* (d^j - b^j) + A^*f)/FFT(A^*A + γ _1∇ ^*∇ ))

By experiments, we find that using equation (<ref>) can indeed get better results than using equation (<ref>). To demonstrate it, we give  both results of using equations (<ref>) and  (<ref>)  in the compared experiments.




In Section <ref>, we compare the experimental results of our model with the related methods (such as split Bregman algorithm (SB) <cit.>, dual tree complex wavelet transform (DTCWT) <cit.>, local contextual hidden markov model (LHMM) <cit.> and the  l_o  minimization in wavelet frame based model ( l_o-WF) <cit.>. In Section <ref>, we perform the comparison on the simulation of image debluring with our model,   the  l_o-WF model and the EDWF<cit.> model. For measuring the image quality quantitatively, we use the index peak signal to noise ratio (PSNR) defined by

    PSNR: = 10log _10(255*255*N/ref - ũ_2^2)

where ref  and ũ  are the true image and recovered image, respectively.
 


 §.§ Image denoising


In this subsection, we compare the performances of our model with  SB, DTCWT, LHMM and  l_o-WF on image denoising. Six images (all of size 256*256) shown in Figure <ref> are  tested. The noisy images are generated by the MATLAB command `imnoise' with `type=Guassian', m = 0  and variance v = 0.01. 



In the implementation of image denoising, the parameters of our algorithm using equation (<ref>) are set as:  λ _1 =λ _2==...=λ _9 = 0.2,  γ _1 = 8,  γ _2=γ _3=...=γ _9 = 4,  tol = 1e - 4; the parameters of our algorithm using equation (<ref>) are set as:  λ _1 = 2,  λ _2=λ _3=...=λ _9 = 1.5,  γ _1 = 12,  γ _2=γ _3=...=γ _9 = 4.5,  tol = 5e - 4. For  SB algorithm, the parameters are set as:  λ  = 19,  gamma  = 1,  tol = 5e - 4. For  DTCWT, the decomposition level  L = 4, threshold  T = 20. For l_o-WF, the piecewise linear B-spline wavelet frame is also used and the decomposition level  L = 1,  λ  = 350,  μ  = 1,  γ  = 0.2,  tol = 5e - 4. The LHMM model has no free parameters. For fair comparison, all the paramters are uniformly set for  the six tested images and we try our best to tune them to get the highest average PSNR.

In Figure <ref>, we show the feature (coefficient) images of a noisy image and its denoised version by our method using equation (<ref>). From Figure <ref>, we can see that the noises in every feature images are all removed.



In Figure <ref> - Figure <ref>, we show the visual comparisons of the results from the four compared methods. From the results of SB, we can see that the recovered images  have a lot of artifacts in the homogeneous regions. The reason is that the regularization parameter λ is chosen optimal for PSNR value. In order to preserve the  sharp edges in images, a relatively small λ should be chosen. However, it  leads to lack of smoothness in  homogeneous regions. This visual effect is especially clear in the `cameraman' image (see the sky and grass). From the results of DTCWT and LCHMM, we can observe that the restored images introduce too much artifacts. one of the reasons is that they both use the decimated wavelet transform. From the last three results, we can see that the denoised images by  l_o-WF and ours have similar visual effects, but our method   has a slightly higher PSNR and less runtime averagely which can be seen in Table <ref>.


Table <ref> lists the PSNR and  total run-time of the compared five methods for image denoising. Since the LCHMM is coded by C++ and the others by matlab, we don't list the runtime of LCHMM.  From Table <ref>, we can see that our model outperforms the other image restoration methods in terms of PSNR, averagely. In addition, from Table <ref>, we can see that the DTCWT is the fastest (It spends less than 1 second for the six test images). However, the quality of the DTCWT is the worst. 

To show the convergence of the our algorithm using equation (<ref>), we plot the relative error ||u^n+1-u^n||/||u^n|| versus the iterations of the six test images in Figure <ref>. From Figure <ref>, we can observe  the relative error is monotonically decreasing with the iterations, which numerically proves that the  algorithm is convergent.










	









 §.§ Image debluring


In this subsection, we apply the proposed model to image deblurring and compare the results with those of  l_o-WF and EDWF. We also test the same six images for image deblurring. The blurred images are generated by convolution with the blur kernel A  and added by a Gaussian noise  η, where A  and η  are generated by the MATLAB command `A = fspecial('motion', 9, 0) ' and `η  = 5*randn(size(u)) ', respectively. 

In the implementation of image deblurring, the parameters of our model using equation (<ref>) are set as:  λ _1 =0.006, λ _2==...=λ _9 = 0.004,  γ _1 = 0.4,  γ _2=γ _3=...=γ _9 = 0.1,  tol = 5e - 4; the parameters of our algorithm using equation (<ref>) are set as:  λ _1 = 0.004,  λ _2=λ _3=...=λ _9 = 0.002,  γ _1 = 0.1,  γ _2=γ _3=...=γ _9 = 0.4,  tol = 5e - 4. For  l_o-WF, the piecewise linear B-spline wavelet frame is also used and the decomposition level   L = 1,  λ  = 30,  μ  = 0.01,  γ  = 0.003,  tol = 5e - 4. For  EDWF, λ=1, γ=1.5, ρ=0.2, v=0.15, μ_1=1, μ_2=1.

In Figure <ref>, the visual results of the three models are presented and Table <ref>, summarizes the PSNR and  total run-time of the three methods. From Table <ref>, we can see that the PSNR of our model using equation (<ref>) is averagely the largest. 








§ CONCLUSIONS
 In this paper, we proposed a vector total variation (VTV) of feature image model for image restoration.  The  VTV imposes different smoothing powers on different features and thus can simultaneously preserve edges and remove noises. Next, the existence of  solution for the model was proved and the split Bregman algorithm was used to solve the model. At last, we used the wavelet filter banks to explicitly define the feature operator and 
presented some experimental results to show its  advantage over the related methods  in both quality and efficiency.





unsrt

