# TV+TV2 Regularization with Nonconvex Sparseness-Inducing Penalty for Image Restoration

Authors: Chengwu Lu; Hua Huang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

In order to restore the high quality image, we propose a compound regularization method which combines a new higher-order extension of total variation (TV+TV2) and a nonconvex sparseness-inducing penalty. Considering the presence of varying directional features in images, we employ the shearlet transform to preserve the abundant geometrical information of the image. The nonconvex sparseness-inducing penalty approach increases robustness to noise and image nonsparsity. In what follows, we present the numerical solution of the proposed model by employing the split Bregman iteration and a novel p-shrinkage operator. And finally, we perform numerical experiments for image denoising, image deblurring, and image reconstructing from incomplete spectral samples. The experimental results demonstrate the efficiency of the proposed restoration method for preserving the structure details and the sharp edges of image.


---

## Body

## 1. Introduction

Image restoration is one of the most classical problems in image processing. During the formation, transmission, or recording process, images may be degraded by noise and blur. Thus, an important problem in image processing is the restoration of the original true image from an observed image. Mathematically, it amounts to estimate an imageu∈ℝN from a measurement
(1)f=𝒜u+w
obtained through a noninvertible linear degradation operator 𝒜:ℝN→ℝM (i.e., an M×N matrix) and contaminated by an additive noise w∈ℝM. Typical degradation operators include identity as in image denoising, incomplete linear projector as in compressive sensing, convolution operator as in image deconvolution problems, 0-1 mask as in image inpainting, and so on.The recovery ofu from the observed image f∈ℝM can be considered as an inverse problem which is ill-posed in general. The information provided by f is not sufficient to ensure existence, uniqueness, and stability of a solution u. It is therefore necessary to regularize the problem by adding a prior constraint on the solution. One typical regularization approach takes the following form:
(2)min⁡u∈ℝN12∥𝒜u-f∥22+λℛ(u),
where the first term is the fidelity function (or fitting function) that forces closeness to data according to (1), ℛ:ℝN→ℝ- is the so-called regularization function (or penalty function) that embodies the priors about the unknown image u with certain characteristics (such as piecewise smoothness or sparseness), and λ>0 is a parameter that controls the trade-off between these two functions. The intuitive meaning of the linear inverse problem (2) is clear: its minimizers reach a compromise between lack of fitness to the observed image (as measured by the fidelity function) and a degree of “undesirability” (as quantified by ℛ(·)).The choice of the regularization functionℛ(·) significantly affects the quality of the restored image. The majority of the existing regularization functions favor simpler images with piecewise constant intensities. Details and fine features in images are not necessarily more complicated. The existing state-of-the-art image restoration models usually are designed for single convex regularization function such as total variation (TV) seminorm or sparseness-inducing penalty. The success and widespread use of the TV seminorm in the past two decades can be mainly attributed to its ability to produce well-preserved and sharp edges. Moreover, its convexity permits the design of efficient algorithms. Although the TV seminorm is extremely popular in a variety of applications, it also gives rise to some undesired effects. In particular, it leads to the well-known staircase effect when applied to images that are not necessarily piecewise constant [1]. To attenuate the staircase effect, there is a growing interest in the literature for replacing TV seminorm by combining a second-order regularization function with the TV seminorm [2–5] or employing a ℓ1 norm of some higher-order differential operator in a standalone way [6–8]. The motivations behind these attempts are to restore potentially a wider range of images, which comprise more than merely piecewise constant regions. More recently, Papafitsoros and Schönlieb constitute a straightforward higher-order extension of the TV seminorm by adding a nonsmooth second-order regularization function; it is referred to as TV+TV2 regularization [9]. TV+TV2 regularization can compete with total generalized variation (TGV) [10], infimal convolution [11], and Euler’s elastic [12], three other state-of-the-art higher-order regularizations. But most of all, TV+TV2 is simple and efficiently numerically solvable.On the other hand, problems of takingℛ(·) as sparseness-inducing penalty have become familiar over the past three decades, particularly in statistical and signal processing contexts. One typical sparseness-inducing penalty problem takes the following form:
(3)min⁡u∈ℝN12∥𝒜u-f∥22+λ∥Ψu∥1,
where ∥·∥1 denotes the ℓ1 norm and Ψ is a sparsity operator, such as a wavelet transform. When Ψ is a discrete gradient operator, the second term leads to TV regularization. In the 1990s, seminal work on the use of ℓ1 sparseness-inducing penalty appeared in the statistics literature, such as the famous basis pursuit denoising (BPDN) criterion [13] and the least absolute shrinkage and selection operator (LASSO, [14]). For brief historical accounts on the use of the ℓ1 penalty in statistics and signal processing, see [15]. Another intriguing new application for the ℓ1 sparseness-inducing penalty is compressive sensing (CS) [16–18] when 𝒜 is the measurement or sampling operator in the optimization problem (3). In addition, it is worth mentioning the numerical results in [19, 20] and theoretical works [21–23]. These results have justified that the robustness of restored images to noise and the nonsparsity can be increased by replacing ℓ1 norm in (3) with the nonconvex nonsmooth ℓp(0<p<1) quasinorm. Moreover, recent researches show that ℓp(0<p<1) quasinorm offers much richer possibilities to restore high quality images with neat edges.To restore high quality image, in the regularization approach (2), one may be desired to encourage the solution to exhibit characteristics that are piecewise smooth and sparse simultaneously. The resulting optimization problem cannot be dealt with by single regularization. It is most naturally expressed by the combination of more than one regularization function, such as compounding total variation and sparseness-inducing penalty [24].In this paper, we consider a new compound regularization, that is, linear combination of TV+TV2 regularization and sparseness-inducing ℓp(0<p<1) quasinorm. Considering the presence of varying directional features in image, we employ the shearlet transform to preserve the abundant geometrical information of the image. As mentioned above, TV+TV2 regularization can effectively attenuate the staircase effect of TV regularization. The nonconvex sparseness-inducing penalty approach increases robustness to noise and image nonsparsity.


## 2. Preliminary

In this section, we will briefly review TV+TV2 regularization and shearlet transform to make this paper self-contained.


### 2.1. TV+TV2 Regularization

In image restoration model (2), the most frequently chosen regularization term is TV seminorm
(4)ℛ(u)=TV(u)=∫Ω|∇u|dxdy,
where |∇u|=:|∂xu|2+|∂yu|2. It is well-known that TV seminorm is tuned towards the preservation of edges and performs very well if the reconstructed image is piecewise constant. However, as an image does not only consist of flat regions and jumps but also possess slanted regions, that is, piecewise linear parts, such regularization procedure will cause undesired staircase effect. One way to reduce this staircase effect is in fact to “round off" TV seminorm by using Huber regularization [25]. However, such procedure can only eliminate the staircase effect in areas with small gradient. Another way of improving TV seminorm is the introduction of higher-order derivatives in the regularization function. Especially in recent years, higher-order versions of nonsmooth image enhancing methods have been considered.To reduce the staircase effect of TV regularization, Chambolle and Lions [4] understood the image u as an addition of two image layers, that is, u=u1+u2 (u1 is the piecewise constant part of u and u2 is the piecewise smooth part), and proposed the inf-convolution regularization function:
(5)ℛ(u)=∫Ω(α|∇u1|+β|∇2u2|)dxdy,
where |∇u1|=:|∂xu1|2+|∂yu1|2 and |∇2u2|=:|∂xxu2|2+|∂yyu2|2+2|∂xyu2|2 is the distributional Hessian of u2. We denote TV(u1)=∫Ω|∇u1|dxdy and TV2(u2)=∫Ω|∇2u2|dxdy. The former is TV regularization and the latter is TV2 regularization, which is total variation of the gradient. In [3, 7, 26–28] the authors also consider using TV2 regularization.Following the idea of the inf-convolution regularization, Chan et al. proposed the following regularization [29]:
(6)ℛ(u)=∫Ω(α|∇u1|+β|Δu2|2)dxdy,p=1,2,
where the Hessian ∇2 used in the inf-convolution regularization is replaced by the Laplacian Δ.Another attempt to combine first- and second-order regularization originates from Chan et al. [2], who consider TV regularization together with weighted versions of the Laplacian. More precisely, they proposed the following regularization:
(7)ℛ(u)=∫Ω(α|∇u|+βg(|∇u|)|Δu|2)dxdy,
where g(·) must be a function with certain growth conditions at infinity in order to allow jumps.It is particularly worth mentioning that, in [10], Bredies et al. proposed the following total generalized variation (TGV), which is recognized as the currently best higher-order regularization function:
(8)ℛ(u)=α∥∇u-v∥1+β∥∇·v∥1,
where the variable v=(v1,v2) varies in the space of symmetric tensors of order 2 and
(9)∥∇·v∥1:=∫Ω(∂xv1)2+(∂yv1)2+(∂xv2)2+(∂yv2)2dxdy.The TGV regularization generalizes the inf-convolution regularization in the sense that it coincides with the inf-convolution regularization whenv only varies in the range of ∇. The TGV regularization improves the results for image denoising over TV regularization and inf-convolution regularization. The interested reader can consult [10] for details.More recently, in [9], the authors choose the following combined first- and second-order regularization function:
(10)ℛ(u)=∫Ω(αϕ1(∇u)+βϕ1(∇2u))dxdy,
where ϕ1:ℝ2→ℝ+ and ϕ2:ℝ4→ℝ+ are all convex function with at most linear growth at infinity. When one considers ϕ1(x)=|x| and ϕ2(x)=|x|, the regularization procedure essentially is a combination of the total variation and the total variation of the gradient; it is referred to as TV+TV2 regularization. The idea of TV+TV2 regularization is to regularize with a fairly large weight α in the first-order term—preserving the jumps as good as possible—and use a small weight β for the second-order term such that staircase artifacts created by the first-order regularization are eliminated without introducing any serious blur in the reconstructed image.In [9], the authors show that for image denoising, deblurring, and inpainting TV+TV2 regularization offers solutions whose quality is not far off from the ones produced by TGV regularization. Moreover, the computational effort needed for its numerical solution is not much more than the one needed for solving TV regularization [30]. For comparison the numerical solution for TGV regularization is in general about ten times slower than TV+TV2 regularization.In addition, it is well-known that TV regularization is not optimal for the restoration of textured regions. In recent years directional representation systems were proposed to cope with curved singularities in images. In particular, curvelets and shearlets provide an optimally sparse approximation in the class of piecewise smooth functions withC2 singularity boundaries. In the next subsection, we will briefly review the discrete shearlet transform that is suited as sparsity operator for the restoration of curved structures.


### 2.2. Shearlet Transform

The classical wavelet transform is unable to provide additional information about the geometry of the set of singularities of image. This is due to the fact that this transform is isotropic. As a result, it has a very limited ability to resolve edges and other distributed discontinuities which usually are very important for image. Shearlet transform is a directional representation system that provides more geometrical information for multidimensional functions such as images [31].Letψ∈L2(ℝ2), Aa=[a00a]  (a∈ℝ+) is parabolic scaling matrix, and Ss=[1s01] (s∈ℝ) is shear matrix. The shearlet system {ψast,a∈ℝ+,s∈ℝ,t∈ℝ2} is generated by applying the operations of dilation, shear transformation, and translation on ψ:
(11)ψast(x):=a-3/4ψ(Aa-1Ss-1(x-t)).The continuous shearlet transform of functionf∈L2(ℝ2) is defined as
(12)Sℋψ(f)(a,s,t)=〈f,ψast〉.The shearlet transform is invertible if the functionψ satisfies the admissibility property
(13)∫ℝ2|ψ^(ω1,ω2)|2|ω1|2dω1dω2,
where ψ^ is Fourier transform of ψ.To obtain a discrete shearlet transform, we consider digital image as function sampled on the grid(1/N)ℐ:(1/N){(m1,m2):mi=0,…,N-1,i=1,2}, where we assume quadratic image for simplicity. Further, let j0:=⌊(1/2)log2N⌋ be the number of considered scales. We discretize the scaling, shear, and translation parameters as
(14)aj:=2-2j=14j,j=0,…,j0-1,sj,k:=k2-j,-2j≤k≤2j,tm:=mN,m∈ℐ.
With these notations the shearlet becomes ψj,k,m:=ψaj,sk,tm(x)=ψ(Aaj-1Ssj,k-1(x-tm)).We select shearlet in this paper due to its directional sensitivity, availability of efficient implementation. Shearlet transform is by far the most optimal sparsifying transform for piecewise smooth image containing rich geometric information such as edges, corners, and spikes. It combines the power of multiscale methods with the ability of extracting geometry of image.


## 2.1. TV+TV2 Regularization

In image restoration model (2), the most frequently chosen regularization term is TV seminorm
(4)ℛ(u)=TV(u)=∫Ω|∇u|dxdy,
where |∇u|=:|∂xu|2+|∂yu|2. It is well-known that TV seminorm is tuned towards the preservation of edges and performs very well if the reconstructed image is piecewise constant. However, as an image does not only consist of flat regions and jumps but also possess slanted regions, that is, piecewise linear parts, such regularization procedure will cause undesired staircase effect. One way to reduce this staircase effect is in fact to “round off" TV seminorm by using Huber regularization [25]. However, such procedure can only eliminate the staircase effect in areas with small gradient. Another way of improving TV seminorm is the introduction of higher-order derivatives in the regularization function. Especially in recent years, higher-order versions of nonsmooth image enhancing methods have been considered.To reduce the staircase effect of TV regularization, Chambolle and Lions [4] understood the image u as an addition of two image layers, that is, u=u1+u2 (u1 is the piecewise constant part of u and u2 is the piecewise smooth part), and proposed the inf-convolution regularization function:
(5)ℛ(u)=∫Ω(α|∇u1|+β|∇2u2|)dxdy,
where |∇u1|=:|∂xu1|2+|∂yu1|2 and |∇2u2|=:|∂xxu2|2+|∂yyu2|2+2|∂xyu2|2 is the distributional Hessian of u2. We denote TV(u1)=∫Ω|∇u1|dxdy and TV2(u2)=∫Ω|∇2u2|dxdy. The former is TV regularization and the latter is TV2 regularization, which is total variation of the gradient. In [3, 7, 26–28] the authors also consider using TV2 regularization.Following the idea of the inf-convolution regularization, Chan et al. proposed the following regularization [29]:
(6)ℛ(u)=∫Ω(α|∇u1|+β|Δu2|2)dxdy,p=1,2,
where the Hessian ∇2 used in the inf-convolution regularization is replaced by the Laplacian Δ.Another attempt to combine first- and second-order regularization originates from Chan et al. [2], who consider TV regularization together with weighted versions of the Laplacian. More precisely, they proposed the following regularization:
(7)ℛ(u)=∫Ω(α|∇u|+βg(|∇u|)|Δu|2)dxdy,
where g(·) must be a function with certain growth conditions at infinity in order to allow jumps.It is particularly worth mentioning that, in [10], Bredies et al. proposed the following total generalized variation (TGV), which is recognized as the currently best higher-order regularization function:
(8)ℛ(u)=α∥∇u-v∥1+β∥∇·v∥1,
where the variable v=(v1,v2) varies in the space of symmetric tensors of order 2 and
(9)∥∇·v∥1:=∫Ω(∂xv1)2+(∂yv1)2+(∂xv2)2+(∂yv2)2dxdy.The TGV regularization generalizes the inf-convolution regularization in the sense that it coincides with the inf-convolution regularization whenv only varies in the range of ∇. The TGV regularization improves the results for image denoising over TV regularization and inf-convolution regularization. The interested reader can consult [10] for details.More recently, in [9], the authors choose the following combined first- and second-order regularization function:
(10)ℛ(u)=∫Ω(αϕ1(∇u)+βϕ1(∇2u))dxdy,
where ϕ1:ℝ2→ℝ+ and ϕ2:ℝ4→ℝ+ are all convex function with at most linear growth at infinity. When one considers ϕ1(x)=|x| and ϕ2(x)=|x|, the regularization procedure essentially is a combination of the total variation and the total variation of the gradient; it is referred to as TV+TV2 regularization. The idea of TV+TV2 regularization is to regularize with a fairly large weight α in the first-order term—preserving the jumps as good as possible—and use a small weight β for the second-order term such that staircase artifacts created by the first-order regularization are eliminated without introducing any serious blur in the reconstructed image.In [9], the authors show that for image denoising, deblurring, and inpainting TV+TV2 regularization offers solutions whose quality is not far off from the ones produced by TGV regularization. Moreover, the computational effort needed for its numerical solution is not much more than the one needed for solving TV regularization [30]. For comparison the numerical solution for TGV regularization is in general about ten times slower than TV+TV2 regularization.In addition, it is well-known that TV regularization is not optimal for the restoration of textured regions. In recent years directional representation systems were proposed to cope with curved singularities in images. In particular, curvelets and shearlets provide an optimally sparse approximation in the class of piecewise smooth functions withC2 singularity boundaries. In the next subsection, we will briefly review the discrete shearlet transform that is suited as sparsity operator for the restoration of curved structures.


## 2.2. Shearlet Transform

The classical wavelet transform is unable to provide additional information about the geometry of the set of singularities of image. This is due to the fact that this transform is isotropic. As a result, it has a very limited ability to resolve edges and other distributed discontinuities which usually are very important for image. Shearlet transform is a directional representation system that provides more geometrical information for multidimensional functions such as images [31].Letψ∈L2(ℝ2), Aa=[a00a]  (a∈ℝ+) is parabolic scaling matrix, and Ss=[1s01] (s∈ℝ) is shear matrix. The shearlet system {ψast,a∈ℝ+,s∈ℝ,t∈ℝ2} is generated by applying the operations of dilation, shear transformation, and translation on ψ:
(11)ψast(x):=a-3/4ψ(Aa-1Ss-1(x-t)).The continuous shearlet transform of functionf∈L2(ℝ2) is defined as
(12)Sℋψ(f)(a,s,t)=〈f,ψast〉.The shearlet transform is invertible if the functionψ satisfies the admissibility property
(13)∫ℝ2|ψ^(ω1,ω2)|2|ω1|2dω1dω2,
where ψ^ is Fourier transform of ψ.To obtain a discrete shearlet transform, we consider digital image as function sampled on the grid(1/N)ℐ:(1/N){(m1,m2):mi=0,…,N-1,i=1,2}, where we assume quadratic image for simplicity. Further, let j0:=⌊(1/2)log2N⌋ be the number of considered scales. We discretize the scaling, shear, and translation parameters as
(14)aj:=2-2j=14j,j=0,…,j0-1,sj,k:=k2-j,-2j≤k≤2j,tm:=mN,m∈ℐ.
With these notations the shearlet becomes ψj,k,m:=ψaj,sk,tm(x)=ψ(Aaj-1Ssj,k-1(x-tm)).We select shearlet in this paper due to its directional sensitivity, availability of efficient implementation. Shearlet transform is by far the most optimal sparsifying transform for piecewise smooth image containing rich geometric information such as edges, corners, and spikes. It combines the power of multiscale methods with the ability of extracting geometry of image.


## 3. Our Model and Algorithm

### 3.1. Unconstrained Formulation for Our Model

In this subsection, we present our new model as follows that employs TV+TV2 regularization with a nonconvex sparseness-inducing penalty based on shearlet transform:
(15)min⁡u∈ℝN12∥𝒜u-f∥22+α∫Ω|∇u|dx+β∫Ω|∇2u|dx+γ∥Sℋψ(u)∥pp,
where 0<p<1. This is a nonconvex optimization problem. One approach to handle this kind of nonconvex optimization problem is the iteratively reweighted least-squares (IRLS) approach which was previously taken by Rao and Kreutz-Delgado [32], but with unimpressive results due to getting trapped in local minima. The alternative herein proposed is to use an operator splitting method and the Bregman iteration method, as well as a novel p-shrinkage operator to solve (15).To facilitate the discussion of our algorithm, without the loss of generality, we represent a gray image asn×n matrix, and the Euclidean space ℝn×n is denoted as V, so u∈V, f∈V. The linear operator 𝒜 is a mapping 𝒜:V→V. Using the above notations, the model (15) can be expressed as follows:
(16)min⁡u∈V12∥𝒜u-f∥22+α∥∇u∥1+β∥∇2u∥1+γ∥Sℋd(u)∥pp,
where Sℋd(u)=(Sℋ1(u),Sℋ2(u),…,SℋJ(u)) is the discrete shearlet transform of u and Sℋj(u)(j=1,2,…,J) denotes its jth subband. ∇u=(Dxu,Dyu) and ∇2u=(Dxxu,Dyyu,Dxyu,Dyxu) are the discrete differential operators of first- and second-order, respectively. By assuming periodic boundary conditions for u, we define Dxu and Dyu as a forward difference operator, respectively, as follows:
(17)(Dxu)i,j=(Dx+u)i,j={ui,j+1-ui,j1≤j≤n-1ui,1-ui,nj=n,(Dyu)i,j=(Dy+u)i,j={ui+1,j-ui,j1≤i≤n-1u1,j-un,ji=n.
We also need the discrete divergence operator div⁡:V2→V that has the adjointness property
(18)-div⁡q·u=q·∇u,∀u∈V,q∈V2.
This property forms the discrete analogue of integration by parts. For q=(q1,q2)∈V2 we define
(19)(div⁡q)i,j=(Dx-q1)i,j+(Dy-q2)i,j,
where Dx- and Dy- are the following backward difference operators:
(20)(Dx-u)i,j={ui,j-ui,j-12≤j≤nui,1-ui,nj=1,(Dy-u)i,j={ui,j-ui-1,j2≤i≤nu1,j-un,ji=1.Analogously we define the second-order discrete differential operatorsDxxu, Dyyu, and Dxyu. These are just the corresponding compositions of the first-order operators:
(21)Dxxu=Dxx+u={ui,n-2ui,j+ui,j+1j=1ui,j-1-2ui,j+ui,j+11<j<nui,j-1-2ui,j+ui,1j=n,(22)Dxxu=Dyy+u={un,j-2ui,j+ui+1,ji=1ui-1,j-2ui,j+ui+1,j1<i<nui-1,j-2ui,j+u1,ji=n,(23)Dxyu=Dxy+u={ui,j-ui+1,j-ui,j+1+ui+1,j+11≤i<n,1≤j<1ui,j-u1,j-ui,j-1-2ui,j+u1,j+1i=n,1≤j<nui,j-ui+1,j-ui,1+ui+1,11≤i<n,j=nui,j-u1,j-ui,1+u1,1i=n,j=n.One can easily check thatDxy+=Dx+(Dy+)=Dy+(Dx+)=Dyx+. As before, we also need the discrete second-order divergence operator div⁡2:V4→V with the adjointness property
(24)div⁡2q·u=q·∇2u,∀u∈V,q∈V4.For aq=(q11,q22,q12,q21)∈V4 we define
(25)(div⁡2q)(i,j)=Dxx-q11(i,j)+Dyy-q22(i,j)+Dxy-q12(i,j)+Dyx-q21(i,j),
where Dxx-=Dxx+, Dyy-=Dyy+, and Dxy-=Dyx- with
(26)Dxy-u={ui,j-ui,n-un,j+un,ni=1,j=1ui,j-ui,j-1-un,j+un,j-1i=1,1<j≤nui,j-ui-1,j-ui,n+ui-1,n1<i≤n,j=1ui,j-ui,j-1-ui-1,j+ui-1,j-11<i≤n,1<j≤n.Remark 1.
By choosing periodic boundary conditions, the action of each of the discrete differential operators can be regarded as a circular convolution of u and allows the use of fast Fourier transform, which our algorithm is relying on.


### 3.2. Constrained Formulation for Our Model

In the following we use what is frequently called the split Bregman iteration to solve (16); see [33, 34]. To do so, we firstly formulate the unconstraint minimization problem (16) as a constraint problem and then apply the Bregman iteration to solve it. We introduce auxiliary variables v=(v1,v2)∈V2, w=(w1,w2,w3,w4)∈V4, and c=(c1,c2,…,cJ)∈VJ such that (16) is equivalent to the following constrained minimization problem:
(27)min⁡u∈V,v∈V2,w∈V4,c∈VJ12∥𝒜u-f∥22+α∥v∥1+β∥w∥1+γ∥c∥pps.t.v=∇u,w=∇2u,c=Sℋd(u).Here, we define(28)∥𝒜u-f∥2=(∑i=1n∑j=1n(𝒜u(i,j)-f(i,j))2)1/2∥v∥1=∑i=1n∑j=1n(v1(i,j)2+v2(i,j)2)1/2∥w∥1=∑i=1n∑j=1n(w11(i,j)2+w22(i,j)2+w12(i,j)2+w21(i,j)2)1/2∥c∥p=(∑i=1n∑j=1n∑l=1Jcl(i,j)p)1/p.Note that the problems (16) and (27) are equivalent. The Bregman iteration that corresponds to (27) is the following:
(29)[uk+1,vk+1,wk+1,ck+1]=arg minu,v,w,c12∥𝒜u-f∥22+α∥v∥1+β∥w∥1+γ∥c∥pp+λ12∥b1k+∇u-v∥22+λ22∥b2k+∇2u-w∥22+λ32∥b3k+Sℋd(u)-c∥22,(30)b1k+1=b1k+∇uk+1-vk+1,(31)b2k+1=b2k+∇2uk+1-wk+1,(32)b3k+1=b3k+Sℋd(uk+1)-ck+1,
where b1k=(b1,1k,b1,2k)∈V2, b2k=(b2,1,1k,b2,2,2k,b2,1,2k,b2,2,1k)∈V4, b3k=(b3,1k,b3,2k,…,b3,Jk)∈VJ, and λ1, λ2, λ3 are positive constants.It is difficult to solve the minimization problem (29) directly; thus we split it into four separate subproblems for u, v, w, and c, each of which can be solved fast as follows.Subproblem  1:
(33)
u
k
+
1
=
arg min
u
1
2
∥
𝒜
u
-
f
∥
2
2
+
λ
1
2
∥
b
1
k
+
∇
u
-
v
k
∥
2
2
+
λ
2
2
∥
b
2
k
+
∇
2
u
-
w
k
∥
2
2
+
λ
3
2
∥
b
3
k
+
S
ℋ
d
(
u
)
-
c
k
∥
2
2
.Subproblem  2:
(34)
v
k
+
1
=
arg min
v
α
∥
v
∥
1
+
λ
1
2
∥
b
1
k
+
∇
u
k
+
1
-
v
∥
2
2
.Subproblem  3:
(35)
w
k
+
1
=
arg min
w
β
∥
w
∥
1
+
λ
2
2
∥
b
2
k
+
∇
2
u
k
+
1
-
w
∥
2
2
.Subproblem  4:
(36)
c
k
+
1
=
arg min
c
γ
∥
c
∥
p
p
+
λ
3
2
∥
b
3
k
+
𝒮
ℋ
d
(
u
k
+
1
)
-
c
∥
2
2
.We now describe how we solve each of the four subproblems (33)–(36).


### 3.3. To Solve Subproblem 1

Subproblem 1 can be solved through its optimality condition. The optimality condition reads as follows:(37)∂u(12∥𝒜u-f∥22+λ12∥b1k+∇u-vk∥22+λ22∥b2k+∇2u-wk∥22+λ32∥b3k+Sℋd(u)-ck∥22)=0.Thejth subband of shearlet transform of u is able to be efficiently implemented as component-wise multiplication with a mask matrix denoted by Hj in frequency domain:
(38)Sℋj(u)=ℱ-1Hjℱu=MHju
with MHl=ℱ-1Hjℱ. Here ℱ and ℱ-1 are the 2D discrete Fourier transform and inverse Fourier transform, respectively.Now by solving the above optimality condition, we get(39)𝒜*𝒜u-λ1(Dx-Dx++Dy-Dy+)u+λ2(Dxx-Dxx++Dyy-Dyy++2Dxy-Dxy+)u+λ3∑l=1JMHl*MHlu=𝒜*f+λ1(Dx-Dx+(b1,1k-v1k)+Dy-Dy+(b1,2k-v2k))-λ2(Dxx-Dxx+(b2,1,1k-w1,1k)+Dyy-Dyy+(b2,2,2k-w2,2k)+2Dxy-Dxy+(b2,1,2k-w1,2k))-λ3∑j=1JMHj*(b3,jk-cjk).
Taking into account that the discrete differential operators act like circular convolutions, we can use the 2D discrete Fourier transform to solve system (39) which has a closed-form solution
(40)FD·ℱ(uk+1)=ℱ(Rk),
where “·” means pointwise multiplication of matrices. Consider
(41)FD≔ℱ(𝒜*𝒜)-λ1(ℱ(Dx-(Dx+))+ℱ(Dy-(Dy+)))+λ2(ℱ(Dxx-(Dxx+))+ℱ(Dyy-(Dyy+))+ℱ(Dxy-(Dxy+)))+λ3(ℱ(∑j=1JMHj*MHj))=ℱ(𝒜*𝒜)-λ1(ℱ(Dx-(Dx+))+ℱ(Dy-(Dy+)))+λ2(ℱ(Dxx-(Dxx+))+ℱ(Dyy-(Dyy+))+ℱ(Dxy-(Dxy+)))+λ3∑j=1J|Hj|2
and Rk is the right hand side of (39). The term FD needs to be calculated only once at the beginning of the algorithm. Thus, we have
(42)uk+1=ℱ-1(ℱ(Rk)FD).


### 3.4. To Solve Subproblems 2 and 3

Subproblems 2 and 3 are similar and the solutions can be explicitly expressed in the form of shrinkage(43)vk+1=Shink(b1k+∇uk+1,αλ1),(44)wk+1=Shink(b2k+∇2uk+1,βλ2),
where Shink(t,σ)=sgn⁡(t)max⁡(|t|-σ,0).


### 3.5. To Solve Subproblem 4

Now we present thep-shrinkage operator to solve subproblem 4. To avoid nondifferentiability caused by the ℓp penalty term in subproblem 4, we consider its smooth approximation problem
(45)ck+1=arg mincγφτ(c)+λ32∥b3k+𝒮ℋd(uk+1)-c∥22,
where φτ(·) is defined by φτ(t)={κ|t|2if|t|≤τ|t|p/p-δif|t|>τ, for τ>0 and any real p. Here the parameters κ and δ are chosen to make φτ(·) a C1 function, and τ will be chosen shortly. For convenience we adopt the unusual convention that, when p=0, |t|p/p will mean log⁡(|t|).Next, we show that there is a functionψτ satisfying
(46)φτ(t)=min⁡s{ψτ(s)+η2∥s-t∥22}.
Let
(47)θτ(t)=|t|22-1ηφτ(t),ϑτ(s)=|s|22+1ηψτ(s);φτ(t) is a smooth approximation of ∥t∥pp; it still is nonconvex, but by taking τ≥η1/(p-2), above the auxiliary function θτ(t) is convex. Simple manipulation shows that, for φτ and ψτ to satisfy (46), it is necessary and sufficient to have θτ=ϑτ*, where ϑτ* is the conjugate (or Legendre-Fenchel transform) of ϑτ defined as
(48)ϑτ*(s)=sup⁡t{〈s,t〉-ϑτ(t)}
which shows how to construct ψτ from φτ through computing ϑτ from θτ. Simple manipulations then show that (46) is satisfied by
(49)ψτ(s)=η(ϑτ(s)-|s|22).
Moreover, the minimizer of (46) will be given by s*=∇θτ(t), which can be shown to be given by what we call the p-shrinkage operator as follows:
(50)Sτp(t)=sgn⁡(t)max⁡{|t|-τ|t|p-1,0}.
This shrinkage operator will not be necessary to compute ψτ(s) explicitly. We now incorporate the above into our algorithm. By using the half-quadratic technique originally introduced by Geman and Yang in [35], the approximation problem (45) becomes an equivalent augmented problem as follows:
(51)ck+1=arg minc,sψτ(s)+η2∥s-c∥22+λ32γ∥b3k+𝒮ℋd(uk+1)-c∥22,
where s=(s1,s2,…,sJ)∈VJ.By lettingη→∞, s is forced to approach c, making the solution to (51) approach that of (45). There is no interaction between c and s in (51), so (51) can become the following two subproblems which are alternately computed:
(52)sk+1=arg minsψτ(s)+η2∥s-ck∥22,(53)ck+1=arg mincη2∥sk+1-c∥22+λ32γ∥b3k+𝒮ℋd(uk+1)-c∥22.Ma et al. [36] use a continuation approach to solve (52), starting with small η values and then increasing them geometrically, using the solution at each stage to initialize the next stage. Instead, we adopt the Bregman iteration approach of Goldstein and Osher [33].For fixedck, the resulting algorithm for (52) is as follows:
(54)sk+1=S1/ηp(ck+b4k),b4k+1=b4k+ck-sk+1,
where b4k=(b4,1k,b4,2k,…,b4,Jk)∈VJ.For fixedsk=1, (53) has a closed-form solution:
(55)ck+1=ηγηγ+λ3sk+1+λ3ηγ+λ3(b3k+Sℋj(uk+1)).The resulting algorithm is summarized in Algorithm1.Algorithm 1:TV+TV2+ℓp sparseness-inducing penalty.
Initialisation:
(1)  uk+1 is given by (42) for solving Subproblem 1;
(2) vk+1 is given by (43) for solving Subproblem 2;
(3) wk+1 is given by (44) for solving Subproblem 3;
(4) ck+1 is given by (54) and (55) for solving Subproblem 3;
(5) Update on b1: b1k+1=b1k+∇uk+1-vk+1;
(6) Update on b2: b2k+1=b2k+∇2uk+1-wk+1;
(7) Update on b3: b3k+1=b3k+Sℋd(uk+1)-ck+1.


## 3.1. Unconstrained Formulation for Our Model

In this subsection, we present our new model as follows that employs TV+TV2 regularization with a nonconvex sparseness-inducing penalty based on shearlet transform:
(15)min⁡u∈ℝN12∥𝒜u-f∥22+α∫Ω|∇u|dx+β∫Ω|∇2u|dx+γ∥Sℋψ(u)∥pp,
where 0<p<1. This is a nonconvex optimization problem. One approach to handle this kind of nonconvex optimization problem is the iteratively reweighted least-squares (IRLS) approach which was previously taken by Rao and Kreutz-Delgado [32], but with unimpressive results due to getting trapped in local minima. The alternative herein proposed is to use an operator splitting method and the Bregman iteration method, as well as a novel p-shrinkage operator to solve (15).To facilitate the discussion of our algorithm, without the loss of generality, we represent a gray image asn×n matrix, and the Euclidean space ℝn×n is denoted as V, so u∈V, f∈V. The linear operator 𝒜 is a mapping 𝒜:V→V. Using the above notations, the model (15) can be expressed as follows:
(16)min⁡u∈V12∥𝒜u-f∥22+α∥∇u∥1+β∥∇2u∥1+γ∥Sℋd(u)∥pp,
where Sℋd(u)=(Sℋ1(u),Sℋ2(u),…,SℋJ(u)) is the discrete shearlet transform of u and Sℋj(u)(j=1,2,…,J) denotes its jth subband. ∇u=(Dxu,Dyu) and ∇2u=(Dxxu,Dyyu,Dxyu,Dyxu) are the discrete differential operators of first- and second-order, respectively. By assuming periodic boundary conditions for u, we define Dxu and Dyu as a forward difference operator, respectively, as follows:
(17)(Dxu)i,j=(Dx+u)i,j={ui,j+1-ui,j1≤j≤n-1ui,1-ui,nj=n,(Dyu)i,j=(Dy+u)i,j={ui+1,j-ui,j1≤i≤n-1u1,j-un,ji=n.
We also need the discrete divergence operator div⁡:V2→V that has the adjointness property
(18)-div⁡q·u=q·∇u,∀u∈V,q∈V2.
This property forms the discrete analogue of integration by parts. For q=(q1,q2)∈V2 we define
(19)(div⁡q)i,j=(Dx-q1)i,j+(Dy-q2)i,j,
where Dx- and Dy- are the following backward difference operators:
(20)(Dx-u)i,j={ui,j-ui,j-12≤j≤nui,1-ui,nj=1,(Dy-u)i,j={ui,j-ui-1,j2≤i≤nu1,j-un,ji=1.Analogously we define the second-order discrete differential operatorsDxxu, Dyyu, and Dxyu. These are just the corresponding compositions of the first-order operators:
(21)Dxxu=Dxx+u={ui,n-2ui,j+ui,j+1j=1ui,j-1-2ui,j+ui,j+11<j<nui,j-1-2ui,j+ui,1j=n,(22)Dxxu=Dyy+u={un,j-2ui,j+ui+1,ji=1ui-1,j-2ui,j+ui+1,j1<i<nui-1,j-2ui,j+u1,ji=n,(23)Dxyu=Dxy+u={ui,j-ui+1,j-ui,j+1+ui+1,j+11≤i<n,1≤j<1ui,j-u1,j-ui,j-1-2ui,j+u1,j+1i=n,1≤j<nui,j-ui+1,j-ui,1+ui+1,11≤i<n,j=nui,j-u1,j-ui,1+u1,1i=n,j=n.One can easily check thatDxy+=Dx+(Dy+)=Dy+(Dx+)=Dyx+. As before, we also need the discrete second-order divergence operator div⁡2:V4→V with the adjointness property
(24)div⁡2q·u=q·∇2u,∀u∈V,q∈V4.For aq=(q11,q22,q12,q21)∈V4 we define
(25)(div⁡2q)(i,j)=Dxx-q11(i,j)+Dyy-q22(i,j)+Dxy-q12(i,j)+Dyx-q21(i,j),
where Dxx-=Dxx+, Dyy-=Dyy+, and Dxy-=Dyx- with
(26)Dxy-u={ui,j-ui,n-un,j+un,ni=1,j=1ui,j-ui,j-1-un,j+un,j-1i=1,1<j≤nui,j-ui-1,j-ui,n+ui-1,n1<i≤n,j=1ui,j-ui,j-1-ui-1,j+ui-1,j-11<i≤n,1<j≤n.Remark 1.
By choosing periodic boundary conditions, the action of each of the discrete differential operators can be regarded as a circular convolution of u and allows the use of fast Fourier transform, which our algorithm is relying on.


## 3.2. Constrained Formulation for Our Model

In the following we use what is frequently called the split Bregman iteration to solve (16); see [33, 34]. To do so, we firstly formulate the unconstraint minimization problem (16) as a constraint problem and then apply the Bregman iteration to solve it. We introduce auxiliary variables v=(v1,v2)∈V2, w=(w1,w2,w3,w4)∈V4, and c=(c1,c2,…,cJ)∈VJ such that (16) is equivalent to the following constrained minimization problem:
(27)min⁡u∈V,v∈V2,w∈V4,c∈VJ12∥𝒜u-f∥22+α∥v∥1+β∥w∥1+γ∥c∥pps.t.v=∇u,w=∇2u,c=Sℋd(u).Here, we define(28)∥𝒜u-f∥2=(∑i=1n∑j=1n(𝒜u(i,j)-f(i,j))2)1/2∥v∥1=∑i=1n∑j=1n(v1(i,j)2+v2(i,j)2)1/2∥w∥1=∑i=1n∑j=1n(w11(i,j)2+w22(i,j)2+w12(i,j)2+w21(i,j)2)1/2∥c∥p=(∑i=1n∑j=1n∑l=1Jcl(i,j)p)1/p.Note that the problems (16) and (27) are equivalent. The Bregman iteration that corresponds to (27) is the following:
(29)[uk+1,vk+1,wk+1,ck+1]=arg minu,v,w,c12∥𝒜u-f∥22+α∥v∥1+β∥w∥1+γ∥c∥pp+λ12∥b1k+∇u-v∥22+λ22∥b2k+∇2u-w∥22+λ32∥b3k+Sℋd(u)-c∥22,(30)b1k+1=b1k+∇uk+1-vk+1,(31)b2k+1=b2k+∇2uk+1-wk+1,(32)b3k+1=b3k+Sℋd(uk+1)-ck+1,
where b1k=(b1,1k,b1,2k)∈V2, b2k=(b2,1,1k,b2,2,2k,b2,1,2k,b2,2,1k)∈V4, b3k=(b3,1k,b3,2k,…,b3,Jk)∈VJ, and λ1, λ2, λ3 are positive constants.It is difficult to solve the minimization problem (29) directly; thus we split it into four separate subproblems for u, v, w, and c, each of which can be solved fast as follows.Subproblem  1:
(33)
u
k
+
1
=
arg min
u
1
2
∥
𝒜
u
-
f
∥
2
2
+
λ
1
2
∥
b
1
k
+
∇
u
-
v
k
∥
2
2
+
λ
2
2
∥
b
2
k
+
∇
2
u
-
w
k
∥
2
2
+
λ
3
2
∥
b
3
k
+
S
ℋ
d
(
u
)
-
c
k
∥
2
2
.Subproblem  2:
(34)
v
k
+
1
=
arg min
v
α
∥
v
∥
1
+
λ
1
2
∥
b
1
k
+
∇
u
k
+
1
-
v
∥
2
2
.Subproblem  3:
(35)
w
k
+
1
=
arg min
w
β
∥
w
∥
1
+
λ
2
2
∥
b
2
k
+
∇
2
u
k
+
1
-
w
∥
2
2
.Subproblem  4:
(36)
c
k
+
1
=
arg min
c
γ
∥
c
∥
p
p
+
λ
3
2
∥
b
3
k
+
𝒮
ℋ
d
(
u
k
+
1
)
-
c
∥
2
2
.We now describe how we solve each of the four subproblems (33)–(36).


## 3.3. To Solve Subproblem 1

Subproblem 1 can be solved through its optimality condition. The optimality condition reads as follows:(37)∂u(12∥𝒜u-f∥22+λ12∥b1k+∇u-vk∥22+λ22∥b2k+∇2u-wk∥22+λ32∥b3k+Sℋd(u)-ck∥22)=0.Thejth subband of shearlet transform of u is able to be efficiently implemented as component-wise multiplication with a mask matrix denoted by Hj in frequency domain:
(38)Sℋj(u)=ℱ-1Hjℱu=MHju
with MHl=ℱ-1Hjℱ. Here ℱ and ℱ-1 are the 2D discrete Fourier transform and inverse Fourier transform, respectively.Now by solving the above optimality condition, we get(39)𝒜*𝒜u-λ1(Dx-Dx++Dy-Dy+)u+λ2(Dxx-Dxx++Dyy-Dyy++2Dxy-Dxy+)u+λ3∑l=1JMHl*MHlu=𝒜*f+λ1(Dx-Dx+(b1,1k-v1k)+Dy-Dy+(b1,2k-v2k))-λ2(Dxx-Dxx+(b2,1,1k-w1,1k)+Dyy-Dyy+(b2,2,2k-w2,2k)+2Dxy-Dxy+(b2,1,2k-w1,2k))-λ3∑j=1JMHj*(b3,jk-cjk).
Taking into account that the discrete differential operators act like circular convolutions, we can use the 2D discrete Fourier transform to solve system (39) which has a closed-form solution
(40)FD·ℱ(uk+1)=ℱ(Rk),
where “·” means pointwise multiplication of matrices. Consider
(41)FD≔ℱ(𝒜*𝒜)-λ1(ℱ(Dx-(Dx+))+ℱ(Dy-(Dy+)))+λ2(ℱ(Dxx-(Dxx+))+ℱ(Dyy-(Dyy+))+ℱ(Dxy-(Dxy+)))+λ3(ℱ(∑j=1JMHj*MHj))=ℱ(𝒜*𝒜)-λ1(ℱ(Dx-(Dx+))+ℱ(Dy-(Dy+)))+λ2(ℱ(Dxx-(Dxx+))+ℱ(Dyy-(Dyy+))+ℱ(Dxy-(Dxy+)))+λ3∑j=1J|Hj|2
and Rk is the right hand side of (39). The term FD needs to be calculated only once at the beginning of the algorithm. Thus, we have
(42)uk+1=ℱ-1(ℱ(Rk)FD).


## 3.4. To Solve Subproblems 2 and 3

Subproblems 2 and 3 are similar and the solutions can be explicitly expressed in the form of shrinkage(43)vk+1=Shink(b1k+∇uk+1,αλ1),(44)wk+1=Shink(b2k+∇2uk+1,βλ2),
where Shink(t,σ)=sgn⁡(t)max⁡(|t|-σ,0).


## 3.5. To Solve Subproblem 4

Now we present thep-shrinkage operator to solve subproblem 4. To avoid nondifferentiability caused by the ℓp penalty term in subproblem 4, we consider its smooth approximation problem
(45)ck+1=arg mincγφτ(c)+λ32∥b3k+𝒮ℋd(uk+1)-c∥22,
where φτ(·) is defined by φτ(t)={κ|t|2if|t|≤τ|t|p/p-δif|t|>τ, for τ>0 and any real p. Here the parameters κ and δ are chosen to make φτ(·) a C1 function, and τ will be chosen shortly. For convenience we adopt the unusual convention that, when p=0, |t|p/p will mean log⁡(|t|).Next, we show that there is a functionψτ satisfying
(46)φτ(t)=min⁡s{ψτ(s)+η2∥s-t∥22}.
Let
(47)θτ(t)=|t|22-1ηφτ(t),ϑτ(s)=|s|22+1ηψτ(s);φτ(t) is a smooth approximation of ∥t∥pp; it still is nonconvex, but by taking τ≥η1/(p-2), above the auxiliary function θτ(t) is convex. Simple manipulation shows that, for φτ and ψτ to satisfy (46), it is necessary and sufficient to have θτ=ϑτ*, where ϑτ* is the conjugate (or Legendre-Fenchel transform) of ϑτ defined as
(48)ϑτ*(s)=sup⁡t{〈s,t〉-ϑτ(t)}
which shows how to construct ψτ from φτ through computing ϑτ from θτ. Simple manipulations then show that (46) is satisfied by
(49)ψτ(s)=η(ϑτ(s)-|s|22).
Moreover, the minimizer of (46) will be given by s*=∇θτ(t), which can be shown to be given by what we call the p-shrinkage operator as follows:
(50)Sτp(t)=sgn⁡(t)max⁡{|t|-τ|t|p-1,0}.
This shrinkage operator will not be necessary to compute ψτ(s) explicitly. We now incorporate the above into our algorithm. By using the half-quadratic technique originally introduced by Geman and Yang in [35], the approximation problem (45) becomes an equivalent augmented problem as follows:
(51)ck+1=arg minc,sψτ(s)+η2∥s-c∥22+λ32γ∥b3k+𝒮ℋd(uk+1)-c∥22,
where s=(s1,s2,…,sJ)∈VJ.By lettingη→∞, s is forced to approach c, making the solution to (51) approach that of (45). There is no interaction between c and s in (51), so (51) can become the following two subproblems which are alternately computed:
(52)sk+1=arg minsψτ(s)+η2∥s-ck∥22,(53)ck+1=arg mincη2∥sk+1-c∥22+λ32γ∥b3k+𝒮ℋd(uk+1)-c∥22.Ma et al. [36] use a continuation approach to solve (52), starting with small η values and then increasing them geometrically, using the solution at each stage to initialize the next stage. Instead, we adopt the Bregman iteration approach of Goldstein and Osher [33].For fixedck, the resulting algorithm for (52) is as follows:
(54)sk+1=S1/ηp(ck+b4k),b4k+1=b4k+ck-sk+1,
where b4k=(b4,1k,b4,2k,…,b4,Jk)∈VJ.For fixedsk=1, (53) has a closed-form solution:
(55)ck+1=ηγηγ+λ3sk+1+λ3ηγ+λ3(b3k+Sℋj(uk+1)).The resulting algorithm is summarized in Algorithm1.Algorithm 1:TV+TV2+ℓp sparseness-inducing penalty.
Initialisation:
(1)  uk+1 is given by (42) for solving Subproblem 1;
(2) vk+1 is given by (43) for solving Subproblem 2;
(3) wk+1 is given by (44) for solving Subproblem 3;
(4) ck+1 is given by (54) and (55) for solving Subproblem 3;
(5) Update on b1: b1k+1=b1k+∇uk+1-vk+1;
(6) Update on b2: b2k+1=b2k+∇2uk+1-wk+1;
(7) Update on b3: b3k+1=b3k+Sℋd(uk+1)-ck+1.


## 4. Numerical Experiments

In this section, we will perform some image restoration experiments using our proposed TV+TV2+ℓp sparseness-inducing penalty model (16). The image restoration scenarios will be denoising, deblurring, and reconstruction from incomplete data. We will use Algorithm 1 to solve (16). We compare our method with the one proposed in [30] (TV method) and the one proposed in [9] (TV+TV2 method). All computations of the present paper are carried out in Matlab 7.0. The results are obtained by running the Matlab codes on an Intel(R) Core (TM) i3 CPU (3.20 GHz, 3.19 GHz) computer with RAM of 1.92 GB.The quality of the restoration results with different methods is compared quantitatively by using the peak signal-to-noise ratio (PSNR), the mean-squared error (MSE), and the structural similarity index (SSIM). They are defined as follows:(56)PSNR=10·log10(255×255MSE),MSE=1mn∥u^-u∥22,
where u and u^ are the expected image of size m×n and the restored one, respectively, and
(57)SSIM=(2μ1μ2+C1)(2σ12+C2)(2μ12μ22+C1)(σ12+σ22+C2),
where μ1 and μ2 are averages of u and u^, respectively, σ1 and σ2 are the variance of u and u^, respectively, and σ12 is the covariance of u and u^. The positive constants C1 and C2 can be thought of as stabilizing constants for near-zero denominator values. The SSIM is a well-known quality metric used to measure the similarity between two images. The SSIM metric was developed by Wang et al. [37] and is based on three specific statistical measures that are much closer to how the human eye perceives differences between images. The SSIM index also assesses the conservation of the structural information of the reconstructed image. Note that the perfect reconstruction would have SSIM value equal to 1. In the following experiments, we will use SSIM map to reveal areas of high or low similarity between two images, the whiter SSIM map, and the closer between the two images.


### 4.1. Denoising Experiments

In our denoising implementation, the linear operator𝒜 equals the identity. We perform experiments to images that are corrupted a white Gaussian noise with standard deviation σ. Note that, for α≠0, β=0, and γ=0, our model (16) corresponds to the classical Rudin-Osher-Fatemi denoising model [30], and, for α≠0, β≠0, and γ=0, it corresponds to the TV+TV2 restoration model [9]. In the following we will examine whether the introduction of the higher-order term (β≠0) and the nonconvex sparseness-inducing penalty term (γ≠0) in the denoising procedure produces results of higher quality.We test our method on the cropped Lena image of size256×256. The image is contaminated with Gaussian noise of variance 0.05. For comparison, we report the results of our proposed method (taking parameters α=0.02, β=0.01, γ=0.05, λ1=0.2, λ2=0.1, λ3=0.0001, and μ=1000), the TV regularization method (taking parameters α=0.02, β=0, γ=0, λ1=0.2, λ2=0, and λ3=0), and the TV+TV2 regularization method (taking parameters α=0.02, β=0.01, γ=0, λ1=0.2, λ2=0.1, and λ3=0). In general, we have found that λ1=10α, λ2=10β, and p≥0.5 usually result in good convergence. Iterations are terminated when the condition ∥uk+1-uk∥2/∥uk+1∥2≤10-3 is met. The denoising results by all three methods are shown in Figures 1 and 2. For better visualization, we include the middle row slices and the patch 3D profile. Note that the nonconvex parameter is p=0.5 for our proposed method. We also have found that the convergence speed will become slow if p<0.5, and the results are roughly the same if p≥0.5. The following deblurring and reconstructing experiments are also this kind of circumstance. Compared with the other methods, the proposed method yields better visual results since it is able to preserve edges and fine features and provide more “natural-looking” images.Figure 1
Comparison of Gaussian noise removal. First column from up to down: original image, noisy image (PSNR = 58.1704, MSE = 0.0998, and SSIM = 0.4210), results of the TV method (PSNR = 63.3717, MSE = 0.0521, and SSIM = 0.8224), results of the TV+TV2 method (PSNR = 63.9644, MSE = 0.0298, and SSIM = 0.8512), and results of our proposed TV+TV2+ℓp method (PSNR = 64.1546, MSE = 0.0161, and SSIM = 0.8907). Second column: middle row slices on the blue line. Third column: 3D close-up intensity profile in the red box. Fourth column: SSIM map.Denoising performance comparisons on the Lena image with different method.
(a)
Iterations versus PSNR (dB)
(b)
Iterations versus SSIM


### 4.2. Deblurring Experiments

In our deblurring implementation,𝒜 denotes a circular convolution with a discrete approximation of a Gaussian kernel of size 11×11. As an additional degradation factor, the initial image is corrupted by additive Gaussian noise of variance 0.001. The quality of the resulting images is measured in terms of an increase in the signal-to-noise ratio (ISNR) measured in dB. The ISNR is defined as
(58)ISNR=10log10(MSE1MSE2),
where MSE1 and MSE2 are the mean-squared error between the degraded and the original image and the restored and the original image, respectively. The test image takes the Goldhill image of size 256×256. In this experiment, we take α=0.00001, β=0.00002, γ=0.001, λ1=0.0001, λ2=0.0002, λ3=0.0001, and μ=5000 for our method, α=0.00001, β=0, γ=0, λ1=0.0001, λ2=0, and λ3=0 for the TV method, and α=0.00001, β=0.00002, γ=0, λ1=0.0001, λ2=0.0002, and λ3=0 for the TV+TV2 model. As the previous experiments, in these deblurring experiments, we still take p=0.5 for our proposed method.For visual appreciation, the corresponding restoration results are reported in Figures3 and 4. The results show the consistent advantage of the proposed method over all the others. The TV solution is characterized by sharper edges, as expected, but also results in staircase effects appearing in the smoother regions of the image. The TV+TV2 solution produces the smoother results and loses some details of the image. The result of our method achieves better visual effect and quantitative index than other methods.Figure 3
Comparison of the blurry and noisy image restoration. First column from up to down: original image, degraded image (PSNR = 59.0771, ISNR = 0.0, MSE = 0.0804, and SSIM = 0.4112), results of the TV method (PSNR = 62.2991, MSE = 0.0401, ISNR = 3.2220, and SSIM = 0.6907), results of the TV+TV2 method (PSNR = 62.8713, MSE = 0.0344, ISNR = 3.7942, and SSIM = 0.8149), and results of our proposed TV+TV2+ℓp method (PSNR = 63.1684, MSE = 0.0286, ISNR = 4.5042, and SSIM = 0.8547). Second column: close-up image in the red box. Third column: 3D close-up intensity profile in the red box. Fourth column: SSIM map.Deblurring performance comparisons on the Goldhill image with different method.
(a)
Iterations versus ISNR (dB)
(b)
Iterations versus SSIM


### 4.3. Reconstructing Experiments from Incomplete Data

Now we show the performance of the proposed method on the incomplete spectral data. Here, the linear operator𝒜 denotes the incomplete measurements on the Fourier transform domain of image. To simulate the incomplete measurements, we start with a clean image, normalize its intensity values to [0,1], and then apply discrete Fourier transform to the clean image. The incomplete measurements are obtained by keeping the samples in some radial lines passing through the center in Fourier domain. In our reconstructing experiments, we use Shepp-Logan phantom image and Brain image as the test images. We consider 22 radial lines k-space sampling (9.36% DCT data). The reconstructing parameters are taken as α=0.3, β=0.01, γ=0.01, and λ3=0.005, for our method, α=0.3, β=0, γ=0, and λ3=0 for the TV method, and α=0.3, β=0.01, γ=0, λ1=0.03, λ2=0.001, and λ3=0 for the TV+TV2 method.The corresponding reconstruction results are reported in Figures5, 6, and 7. Because the phantom has a very sparse gradient, all methods obtain good visual effect in the reconstruction of Shepp-Logan phantom image from 9.36% DCT data. Although our proposed method achieves a better quantitative index than the TV and the TV+TV2, the results look very similar and it is difficult to observe any significant visual differences. However, the proposed method demonstrates the distinct advantage in the reconstruction of the Brain image.Figure 5
Comparison of Shepp-Logan phantom image reconstruction from 9.36% DCT data. First column from up to down: original image, back projection image (PSNR = 56.9194, MSE = 0.1322, and SSIM = 0.2790), results of the TV method (PSNR = 65.1352, MSE = 0.0274, and SSIM = 0.9768), results of the TV+TV2 method (PSNR = 65.3658, MSE = 0.0189, and SSIM = 0.9820), and results of our proposed method (PSNR = 66.6062, MSE = 0.0142, and SSIM = 0.9896). Second column: close-up image in the red box. Third column: 3D close-up intensity profile in the red box. Fourth column: SSIM map.Figure 6
Comparison of reconstruction of Brain image from 9.36% DCT data. First column from up to down: original image, back projection image (PSNR = 58.0251, MSE = 0.1025, and SSIM = 0.3919), results of the TV method (PSNR = 59.8861, MSE = 0.0668, and SSIM = 0.5376), results of the TV+TV2 method (PSNR = 62.5794, MSE = 0.0395, and SSIM = 0.7519), and results of our proposed TV+TV2+ℓp method (PSNR = 64.4744, MSE = 0.0232, and SSIM = 0.8713). Second column: close-up image in the red box. Third column: 3D close-up intensity profile in the red box. Fourth column: SSIM map.Reconstruction performance (Iterations versus MSE) comparisons on two test images.
(a)
Shepp-Logan phantom image
(b)
Brain image


## 4.1. Denoising Experiments

In our denoising implementation, the linear operator𝒜 equals the identity. We perform experiments to images that are corrupted a white Gaussian noise with standard deviation σ. Note that, for α≠0, β=0, and γ=0, our model (16) corresponds to the classical Rudin-Osher-Fatemi denoising model [30], and, for α≠0, β≠0, and γ=0, it corresponds to the TV+TV2 restoration model [9]. In the following we will examine whether the introduction of the higher-order term (β≠0) and the nonconvex sparseness-inducing penalty term (γ≠0) in the denoising procedure produces results of higher quality.We test our method on the cropped Lena image of size256×256. The image is contaminated with Gaussian noise of variance 0.05. For comparison, we report the results of our proposed method (taking parameters α=0.02, β=0.01, γ=0.05, λ1=0.2, λ2=0.1, λ3=0.0001, and μ=1000), the TV regularization method (taking parameters α=0.02, β=0, γ=0, λ1=0.2, λ2=0, and λ3=0), and the TV+TV2 regularization method (taking parameters α=0.02, β=0.01, γ=0, λ1=0.2, λ2=0.1, and λ3=0). In general, we have found that λ1=10α, λ2=10β, and p≥0.5 usually result in good convergence. Iterations are terminated when the condition ∥uk+1-uk∥2/∥uk+1∥2≤10-3 is met. The denoising results by all three methods are shown in Figures 1 and 2. For better visualization, we include the middle row slices and the patch 3D profile. Note that the nonconvex parameter is p=0.5 for our proposed method. We also have found that the convergence speed will become slow if p<0.5, and the results are roughly the same if p≥0.5. The following deblurring and reconstructing experiments are also this kind of circumstance. Compared with the other methods, the proposed method yields better visual results since it is able to preserve edges and fine features and provide more “natural-looking” images.Figure 1
Comparison of Gaussian noise removal. First column from up to down: original image, noisy image (PSNR = 58.1704, MSE = 0.0998, and SSIM = 0.4210), results of the TV method (PSNR = 63.3717, MSE = 0.0521, and SSIM = 0.8224), results of the TV+TV2 method (PSNR = 63.9644, MSE = 0.0298, and SSIM = 0.8512), and results of our proposed TV+TV2+ℓp method (PSNR = 64.1546, MSE = 0.0161, and SSIM = 0.8907). Second column: middle row slices on the blue line. Third column: 3D close-up intensity profile in the red box. Fourth column: SSIM map.Denoising performance comparisons on the Lena image with different method.
(a)
Iterations versus PSNR (dB)
(b)
Iterations versus SSIM


## 4.2. Deblurring Experiments

In our deblurring implementation,𝒜 denotes a circular convolution with a discrete approximation of a Gaussian kernel of size 11×11. As an additional degradation factor, the initial image is corrupted by additive Gaussian noise of variance 0.001. The quality of the resulting images is measured in terms of an increase in the signal-to-noise ratio (ISNR) measured in dB. The ISNR is defined as
(58)ISNR=10log10(MSE1MSE2),
where MSE1 and MSE2 are the mean-squared error between the degraded and the original image and the restored and the original image, respectively. The test image takes the Goldhill image of size 256×256. In this experiment, we take α=0.00001, β=0.00002, γ=0.001, λ1=0.0001, λ2=0.0002, λ3=0.0001, and μ=5000 for our method, α=0.00001, β=0, γ=0, λ1=0.0001, λ2=0, and λ3=0 for the TV method, and α=0.00001, β=0.00002, γ=0, λ1=0.0001, λ2=0.0002, and λ3=0 for the TV+TV2 model. As the previous experiments, in these deblurring experiments, we still take p=0.5 for our proposed method.For visual appreciation, the corresponding restoration results are reported in Figures3 and 4. The results show the consistent advantage of the proposed method over all the others. The TV solution is characterized by sharper edges, as expected, but also results in staircase effects appearing in the smoother regions of the image. The TV+TV2 solution produces the smoother results and loses some details of the image. The result of our method achieves better visual effect and quantitative index than other methods.Figure 3
Comparison of the blurry and noisy image restoration. First column from up to down: original image, degraded image (PSNR = 59.0771, ISNR = 0.0, MSE = 0.0804, and SSIM = 0.4112), results of the TV method (PSNR = 62.2991, MSE = 0.0401, ISNR = 3.2220, and SSIM = 0.6907), results of the TV+TV2 method (PSNR = 62.8713, MSE = 0.0344, ISNR = 3.7942, and SSIM = 0.8149), and results of our proposed TV+TV2+ℓp method (PSNR = 63.1684, MSE = 0.0286, ISNR = 4.5042, and SSIM = 0.8547). Second column: close-up image in the red box. Third column: 3D close-up intensity profile in the red box. Fourth column: SSIM map.Deblurring performance comparisons on the Goldhill image with different method.
(a)
Iterations versus ISNR (dB)
(b)
Iterations versus SSIM


## 4.3. Reconstructing Experiments from Incomplete Data

Now we show the performance of the proposed method on the incomplete spectral data. Here, the linear operator𝒜 denotes the incomplete measurements on the Fourier transform domain of image. To simulate the incomplete measurements, we start with a clean image, normalize its intensity values to [0,1], and then apply discrete Fourier transform to the clean image. The incomplete measurements are obtained by keeping the samples in some radial lines passing through the center in Fourier domain. In our reconstructing experiments, we use Shepp-Logan phantom image and Brain image as the test images. We consider 22 radial lines k-space sampling (9.36% DCT data). The reconstructing parameters are taken as α=0.3, β=0.01, γ=0.01, and λ3=0.005, for our method, α=0.3, β=0, γ=0, and λ3=0 for the TV method, and α=0.3, β=0.01, γ=0, λ1=0.03, λ2=0.001, and λ3=0 for the TV+TV2 method.The corresponding reconstruction results are reported in Figures5, 6, and 7. Because the phantom has a very sparse gradient, all methods obtain good visual effect in the reconstruction of Shepp-Logan phantom image from 9.36% DCT data. Although our proposed method achieves a better quantitative index than the TV and the TV+TV2, the results look very similar and it is difficult to observe any significant visual differences. However, the proposed method demonstrates the distinct advantage in the reconstruction of the Brain image.Figure 5
Comparison of Shepp-Logan phantom image reconstruction from 9.36% DCT data. First column from up to down: original image, back projection image (PSNR = 56.9194, MSE = 0.1322, and SSIM = 0.2790), results of the TV method (PSNR = 65.1352, MSE = 0.0274, and SSIM = 0.9768), results of the TV+TV2 method (PSNR = 65.3658, MSE = 0.0189, and SSIM = 0.9820), and results of our proposed method (PSNR = 66.6062, MSE = 0.0142, and SSIM = 0.9896). Second column: close-up image in the red box. Third column: 3D close-up intensity profile in the red box. Fourth column: SSIM map.Figure 6
Comparison of reconstruction of Brain image from 9.36% DCT data. First column from up to down: original image, back projection image (PSNR = 58.0251, MSE = 0.1025, and SSIM = 0.3919), results of the TV method (PSNR = 59.8861, MSE = 0.0668, and SSIM = 0.5376), results of the TV+TV2 method (PSNR = 62.5794, MSE = 0.0395, and SSIM = 0.7519), and results of our proposed TV+TV2+ℓp method (PSNR = 64.4744, MSE = 0.0232, and SSIM = 0.8713). Second column: close-up image in the red box. Third column: 3D close-up intensity profile in the red box. Fourth column: SSIM map.Reconstruction performance (Iterations versus MSE) comparisons on two test images.
(a)
Shepp-Logan phantom image
(b)
Brain image


## 5. Conclusion

We present a TV+TV2+ℓp sparseness-inducing penalty image restoration model. The TV+TV2 regularization combines convex functions of the total variation and the total variation of the first derivatives. It leads to a significant reduction of the staircase effects commonly seen in traditional TV based image processing algorithms. The introduction of the nonconvex ℓp sparseness-inducing penalty based on the shearlet transform increases robustness to image nonsparsity. Our algorithm makes use of previous operator splitting methods and the Bregman iteration method, as well as a novel p-shrinkage operator. Numerical results show the consistent excellence of our proposed method in preserving diverse details during the restoration.Future researches aim at investigating the theoretical rules on selecting the regularization parameters and the nonconvex parameterp(0<p<1) provided by our model.



---

*Source: 790547-2014-03-04.xml*