# Fractional-Order Total Variation Image Restoration Based on Primal-Dual Algorithm

Authors: Dali Chen; YangQuan Chen; Dingyu Xue
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

This paper proposes a fractional-order total variation image denoising algorithm based on the primal-dual method, which provides a much more elegant and effective way of treating problems of the algorithm implementation, ill-posed inverse, convergence rate, and blocky effect. The fractional-order total variation model is introduced by generalizing the first-order model, and the corresponding saddle-point and dual formulation are constructed in theory. In order to guaranteeO1/N2 convergence rate, the primal-dual algorithm was used to solve the constructed saddle-point problem, and the final numerical procedure is given for image denoising. Finally, the experimental results demonstrate that the proposed
methodology avoids the blocky effect, achieves state-of-the-art performance, and guarantees O1/N2 convergence rate.


---

## Body

## 1. Introduction

Since the work of Rudin et al. [1], total variation (TV) minimization problems arise in many image processing applications for regularizing inverse problems where one expects the recovered image or signal to be piecewise constant [2–5]. The typical total variation model, for example, ROF model, has been proved to be able to achieve a good tradeoff between edge preservation and noise removal [1]. However, it tends to produce the so-called blocky (staircase) effects on the images because it favors a piecewise constant solution in bounded variation (BV) space [6]. In order to deal with blocky effects, the modification of TV model, which generalizes the differential order in regularization term, has aroused the more and more attentions of numerous scholars.The improved methods of TV model are divided into two kinds: the high-order derivative and the fractional-order derivative. The first one replaces the first-order derivative in regularization term by the high-order derivative. For example, a fourth-order partial differential equation-(PDE-based) denoising mode was proposed by [7], in which the regularized solution is obtained by solving the minimization of potential function of second-order derivative of the image. It has been proved that this model is able to deal with the blocky problem. However, it tends to cause the sign of uplifting effect and formation of artifacts around edges. For this problem, an improved fourth-order PDE model, replacing the Laplacian operator of diffusivity function by the gradient operator, was proposed in [8].In this paper, our interest focuses on the second generalization which deals with fractional-order differentiation. Fractional calculus is a rapidly growing mathematical discipline, which provides an important tool for nonlocal field theories [9]. Recently, it has been greatly studied in computer vision [10–16], and the main reason for this development is the expectation that the use of this theory will lead to a much more elegant and effective way of treating problems of blocky effect and detailed information protection. Specially, the fractional-order total variation (TV) models play an important role for image denoising, inpainting, and motion estimation [17–21]. So far, the methods adopted to deal with the fractional-order total variation problem are divided into three kinds: (1) solving the associated Euler-Lagrange equation, which is a nonlinear partial differential equation (PDE) [17, 18]; (2) using the methods based on duality [6]; (3) using the method based on majorization-minimization (MM) algorithm [22]. Although the convergence of these methods has been studied in the corresponding literatures, the rate of convergence for them has not been discussed yet. For this problem, in this paper, a primal-dual algorithm is used to solve the fractional-order total variation problem, which is able to guarantee O(1/N2) convergence rate.The primal-dual algorithm was first presented by [23] and named as the primal-dual hybrid gradient (PDHG) method in [24]. In this method, each iteration updates both a primal and a dual variable. It is thus able to avoid some of the difficulties that arise when working only on the primal or dual side [25, 26]. The convergence of the primal-dual algorithm has been studied in [27, 28]. Recently a unified form of primal-dual algorithm was presented by [29, 30], which demonstrated that, in some case, these algorithms can also achieve the O(1/N2) rate of convergence. In our paper, a new image denoising method is proposed, in which the primal-dual algorithm is used to solve the fractional-order total variation denoising model. The proposed denoising method is able to avoid the blocky effect, achieves state-of-the-art performance, and guarantees the O(1/N2) rate of convergence.This paper is organized as follows. Section1 introduces prior work, focusing on the main problems with existing methods that are addressed by our model. In Section 2, the fractional-order total variation denoising model is described, and the corresponding saddle-point and dual formulation are constructed in theory. Based on this, the primal-dual algorithm was used to solve the constructed saddle-point problem, and the final numerical implementation is presented for image denoising. Experimental evaluation is presented in Section 3 and the paper is concluded in Section 4.


## 2. Fractional-Order Total Variation Model and Primal-Dual Algorithm

### 2.1. Model Description

Letfi,j=ui,j+vi,j denote the observed noisy image, where (i,j)T denotes the location with a rectangular image domain Ω∈Rn×n and v is white Gaussian noise. The typical total variation (TV) denoising model estimates the desired clean image ui,j by solving the following finite-dimensional optimization problem:
(1)u^=argminu∈Ω{E(u):=TV(u)+λ2∥f-u∥22},
where TV(u) is the regularization term, ∥f-u∥22 is the data fidelity term, ∥u∥ν is ν-norm of u, and λ is regularization parameter which controls the degree of smoothing. In this paper, we consider a fractional-order total variation model, defined as
(2)u^=argminu∈Ω{E(u):=TVα(u)+λ2∥f-u∥22},
where TVα(u) is obtained by the following formula:
(3)TVα(u)=∥(D1αu)2+(D2αu)2∥1,
where (D1αu)2+(D2αu)2 is the discretization of |Dαu|, Dα is the fractional-order derivative operator, and Dαu:=(D1αu,D2αu)T. D1α and D2α are linear operators corresponding to horizontal and vertical fractional-order derivative.From Grünwald-Letnikov fractional derivative definition [31], the finite fractional-order forward difference can be obtained by
(4)D1αui,j=∑k=0K-1‍wk(α)ui+k,j,D2αui,j=∑k=0K-1‍wk(α)ui,j+k,
where wk(α)=(-1)k+1𝒞kα, 𝒞kα=Γ(α+1)/(Γ(k+1)Γ(α-k+1)) denotes the generalized binomial coefficient, and Γ(x) is the Gamma function. In addition, the coefficients w(α) can also be obtained recursively from
(5)w0(α)=-1,wk(α)=(1-α+1k)wk-1(α),k=1,2,….
When α=1, wk1=0 for k>1 and (4) is the first-order forward derivative as usual.To describe (2) in matrix algebra language, we reord the image matrix u and f row wisely into the vector x and z, associating the (i,j) element of the two-dimensional structure with the element (j-1)n+i of the vector structure, x(j-1)n+i=ui,j, and z(j-1)n+i=fi,j. We have x∈RN and z∈RN, where N=n2. The (i,j) component of the fractional-order derivative Dαu can thus be represented as a multiplication of the vector x∈RN by a matrix Am∈R2×N, for m=1,2,…,N:
(6)Amx={(∑k=0K-1‍wk(α)xm+k,∑k=0K-1‍wk(α)xm+nk,)Tif(mmodn)≠0andm≤N-n(0,∑k=0K-1‍wk(α)xm+nk)T,if(mmodn)=0andm≤N-n(∑k=0K-1‍wk(α)xm+k,0)T,if(mmodn)≠0andm>N-n(0,0)T,if(mmodn)=0andm>N-n.
From this definition, the discrete version of the prime fractional-order total variation (FOTV) model (2) can be written as
(7)minx∈RN{E(x):=∑m=1N‍∥Amx∥2+λ2∥x-z∥22}.
Define a norm ∥v∥E=∑m=1N‍v12+v22 where v=(v1,v2)T, and (7) can be rewritten in term of ∥·∥E as
(8)minx∈RN{E(x):=∥Ax∥E+λ2∥x-z∥22},
Here A=[A1,A2,…,AN]T∈R2N×N.


### 2.2. Dual Problem and Saddle-Point Problem

For any vector,pm=(pm1,pm2)T∈R2 and p=(p1,p2,…,pN)T∈R2N. We have ∥Ax∥E=max∥p∥E*≤1〈Ax,p〉, where ∥·∥E* is the dual norm of norm ∥·∥E, and the ∥·∥E and ∥·∥E* are analogous to ∥·∥1 and ∥·∥∞, respectively.Let the functionδ(y):R2N→R be the conjugate of ∥Ax∥E; we have
(9)δ(y)=supx〈y,Ax〉-∥Ax∥E=supx〈y,Ax〉-max∥p∥E*≤1〈Ax,p〉={0,if∥v∥E*≤1∞,ifotherwise,
in which case we can interpret δ(·) as the indicator function for the unit ball in the dual norm, and y=(y1,y2,…,yN)T∈R2N with ym=(ym1,ym2)T∈R2. Since F=F** if F is closed and convex, we have
(10)∥Ax∥E=supy∈R2N〈y,Ax〉-δ(y).
Substituting (10) into (8), we can obtain the saddle-point formulation of FOTV model defined by
(11)minx∈RNsupy∈R2N〈y,Ax〉-δ(y)+λ2∥x-z∥22.
The minimization problem of (11) can be solved exactly as
(12)x=z-1λATy.
Substituting (12) into (11) yields the following dual problem:
(13)maxy∈R2N-(∥ATy-λz∥22+δ(y)),
where -AT is the corresponding discretization of the fractional-order divergence, and when α=1, it is the discrete divergence as usual.For primal-dual feasible pair(x,y), the partial primal-dual gap G(x,y) is defined by
(14)G(x,y)={maxy′∈R2N〈y′,Ax〉-δ(y′)+λ2∥x-z∥22}-{minx′∈RN〈y,Ax′〉-δ(y)+λ2∥x′-z∥22}.
If (x^,y^) is the saddle point of problem (11), then x^ is primal optimal, y^ is dual optimal [32], and we have
(15)G(x,y)≥{〈y^,Ax〉-δ(y^)+λ2∥x-z∥22}-{〈y,Ax^〉-δ(y)+λ2∥x^-z∥22}≥0,
which vanishes only if (x,y) is itself the saddle point [29]. The primal-dual gap G(x,y) is a measure of closeness of the primal-dual (x,y) to the primal-dual solution, and we use it to design the stopping criterion for our numerical algorithm in this paper.


### 2.3. Primal-Dual Method for Solving Saddle-Point Problem

The first-order primal-dual method summarized in [29] for convex problem was applied in this paper to solve the saddle-point problem described by (11), which is described by Algorithm 1.Algorithm 1
Step 1. Choose τ0,σ0>0 with τ0σ0L2≤1, x0∈RN, y0∈R2N, and x¯0=x0.
Step 2. Update xt, yt, x¯t, θt, τt, σt as follows:
y
t
+
1
=
arg
min
y
∈
R
2
N
〈
y
,
-
A
x
¯
t
〉
+
δ
(
y
)
+
1
2
σ
t
∥
y
-
y
t
∥
2
2
,
(
A
)
x
t
+
1
=
arg
min
x
∈
R
N
λ
2
∥
x
-
z
∥
2
2
+
〈
A
x
,
y
t
+
1
〉
+
1
2
τ
t
∥
x
-
x
t
∥
2
2
,
(
B
)
θ
t
=
1
1
+
2
γ
τ
t
,
τ
t
+
1
=
θ
t
τ
t
,
σ
t
+
1
=
σ
t
θ
t
,
x
¯
t
+
1
=
x
t
+
1
+
θ
t
(
x
t
+
1
-
x
t
)
,LetL=∥A∥ and τtσtL2≤1; then there exists a saddle-point (x^,y^) such that xt→x^ and yt→y^. Since (λ/2)∥x-z∥22 is uniformly convex such that it has a Lipschitz continuous gradient, the output (xt,yt) possesses a rate of convergence of O(1/N2) [29]. More specifically, we assume that the constant stepsizes are used; that is, τt=τ>0, σt=σ>0, and θt=θ for all t≥0. If τσL2<1; then the O(1/N) convergence rate can be guaranteed for θ=1, and O(1/N) convergence rate can be guaranteed for θ=0 [30].To summarize, our entire noise removal algorithm in a form of a pseudocode is done in following:(1)
initialization:
(i)
for a given noisy imagef of size n×n, we reord the image f row wisely into the vector z∈RN with N=n×n, z(j-1)n-i=fi,j;(ii)
compute the fractional-order discrete operatorA according to (6);(iii)
initialize the regularization parameterλ, iteration number T, and maximum permissible error Empe;(iv)
sett=0, x0=z, L=∥A∥, τ0=1/L2, σ0=1/L2, and γ=0.7λ;(2)
iteration: computex by the following steps.Step 1.
Solve the problem (A) by(16)ymt+1=ymt+σtAmx¯tmax{1,∥ymt+σtAmx¯t∥2},m=1,2,…,N.Step 2.
Substituteyt+1 into the (B) and solve the problem (B) by
(17)xt+1=xt-τtATyt+1+τtλz1+τtλ.Step 3.
Computeθt=1/1+2γτt, τt+1=τtθt, and σt+1=σt/θt.Step 4.
Updatex¯ by x¯t+1=(1+θt)xt+1-θtxt.Step 5.
Compute the primal-dual gapgt+1 according to (15) as follows:
(18)gt+1=∑m=1N‍∥Amxt+1∥2-(yt+1)TAxt+1+λ2∥xt+1-z+1λATyt+1∥22.Step 6.
If (gt+1<Empe) or (t>T), then we terminate the iteration and output xt+1; otherwise, go back to Step 1.Due to the simplicity of our model, it is easy to be implemented and can be effectively accelerated on parallel hardware such as field programmable gata array (FPGA) and graphics processing unit (GPU). In order to employ the huge computational power and the parallel processing capabilities of FPGA to obtain a fully accelerated implementation of our denoising method, the numerical method should be working on the regular grids. For this purpose, the implementation of our algorithm on FPGA is described as the following strategy. First, an image frame is captured and read in the processor from a live camera. For the captured framef∈Rn×n, the kth iteration of our method performs the following updates on prime variable u∈Rn×n and dual variable p∈Rn×n×2:
(19)pk+1=pk+σkDαu¯kmax{1,|pk+σkDαu¯k|},(20)uk+1=uk-τkdivαpk+1+τkλf1+τkλ,
where the dual variable p=(p1,p2)T corresponds with y and the discrete fractional divergence operator divαis adjoint to the discrete fractional gradient operator Dα. Note that (19) and (20), respectively, correspond with the solutions of problem (A) and (B). The updates of the other parameters are the same as the description in the foregoing pseudocode, and the number of iteration is fixed to ensure the correct timing sequence. After the iteration, the denoised image u is outputted, the next image frame is captured, and the same process is repeated until the assignment is over.


## 2.1. Model Description

Letfi,j=ui,j+vi,j denote the observed noisy image, where (i,j)T denotes the location with a rectangular image domain Ω∈Rn×n and v is white Gaussian noise. The typical total variation (TV) denoising model estimates the desired clean image ui,j by solving the following finite-dimensional optimization problem:
(1)u^=argminu∈Ω{E(u):=TV(u)+λ2∥f-u∥22},
where TV(u) is the regularization term, ∥f-u∥22 is the data fidelity term, ∥u∥ν is ν-norm of u, and λ is regularization parameter which controls the degree of smoothing. In this paper, we consider a fractional-order total variation model, defined as
(2)u^=argminu∈Ω{E(u):=TVα(u)+λ2∥f-u∥22},
where TVα(u) is obtained by the following formula:
(3)TVα(u)=∥(D1αu)2+(D2αu)2∥1,
where (D1αu)2+(D2αu)2 is the discretization of |Dαu|, Dα is the fractional-order derivative operator, and Dαu:=(D1αu,D2αu)T. D1α and D2α are linear operators corresponding to horizontal and vertical fractional-order derivative.From Grünwald-Letnikov fractional derivative definition [31], the finite fractional-order forward difference can be obtained by
(4)D1αui,j=∑k=0K-1‍wk(α)ui+k,j,D2αui,j=∑k=0K-1‍wk(α)ui,j+k,
where wk(α)=(-1)k+1𝒞kα, 𝒞kα=Γ(α+1)/(Γ(k+1)Γ(α-k+1)) denotes the generalized binomial coefficient, and Γ(x) is the Gamma function. In addition, the coefficients w(α) can also be obtained recursively from
(5)w0(α)=-1,wk(α)=(1-α+1k)wk-1(α),k=1,2,….
When α=1, wk1=0 for k>1 and (4) is the first-order forward derivative as usual.To describe (2) in matrix algebra language, we reord the image matrix u and f row wisely into the vector x and z, associating the (i,j) element of the two-dimensional structure with the element (j-1)n+i of the vector structure, x(j-1)n+i=ui,j, and z(j-1)n+i=fi,j. We have x∈RN and z∈RN, where N=n2. The (i,j) component of the fractional-order derivative Dαu can thus be represented as a multiplication of the vector x∈RN by a matrix Am∈R2×N, for m=1,2,…,N:
(6)Amx={(∑k=0K-1‍wk(α)xm+k,∑k=0K-1‍wk(α)xm+nk,)Tif(mmodn)≠0andm≤N-n(0,∑k=0K-1‍wk(α)xm+nk)T,if(mmodn)=0andm≤N-n(∑k=0K-1‍wk(α)xm+k,0)T,if(mmodn)≠0andm>N-n(0,0)T,if(mmodn)=0andm>N-n.
From this definition, the discrete version of the prime fractional-order total variation (FOTV) model (2) can be written as
(7)minx∈RN{E(x):=∑m=1N‍∥Amx∥2+λ2∥x-z∥22}.
Define a norm ∥v∥E=∑m=1N‍v12+v22 where v=(v1,v2)T, and (7) can be rewritten in term of ∥·∥E as
(8)minx∈RN{E(x):=∥Ax∥E+λ2∥x-z∥22},
Here A=[A1,A2,…,AN]T∈R2N×N.


## 2.2. Dual Problem and Saddle-Point Problem

For any vector,pm=(pm1,pm2)T∈R2 and p=(p1,p2,…,pN)T∈R2N. We have ∥Ax∥E=max∥p∥E*≤1〈Ax,p〉, where ∥·∥E* is the dual norm of norm ∥·∥E, and the ∥·∥E and ∥·∥E* are analogous to ∥·∥1 and ∥·∥∞, respectively.Let the functionδ(y):R2N→R be the conjugate of ∥Ax∥E; we have
(9)δ(y)=supx〈y,Ax〉-∥Ax∥E=supx〈y,Ax〉-max∥p∥E*≤1〈Ax,p〉={0,if∥v∥E*≤1∞,ifotherwise,
in which case we can interpret δ(·) as the indicator function for the unit ball in the dual norm, and y=(y1,y2,…,yN)T∈R2N with ym=(ym1,ym2)T∈R2. Since F=F** if F is closed and convex, we have
(10)∥Ax∥E=supy∈R2N〈y,Ax〉-δ(y).
Substituting (10) into (8), we can obtain the saddle-point formulation of FOTV model defined by
(11)minx∈RNsupy∈R2N〈y,Ax〉-δ(y)+λ2∥x-z∥22.
The minimization problem of (11) can be solved exactly as
(12)x=z-1λATy.
Substituting (12) into (11) yields the following dual problem:
(13)maxy∈R2N-(∥ATy-λz∥22+δ(y)),
where -AT is the corresponding discretization of the fractional-order divergence, and when α=1, it is the discrete divergence as usual.For primal-dual feasible pair(x,y), the partial primal-dual gap G(x,y) is defined by
(14)G(x,y)={maxy′∈R2N〈y′,Ax〉-δ(y′)+λ2∥x-z∥22}-{minx′∈RN〈y,Ax′〉-δ(y)+λ2∥x′-z∥22}.
If (x^,y^) is the saddle point of problem (11), then x^ is primal optimal, y^ is dual optimal [32], and we have
(15)G(x,y)≥{〈y^,Ax〉-δ(y^)+λ2∥x-z∥22}-{〈y,Ax^〉-δ(y)+λ2∥x^-z∥22}≥0,
which vanishes only if (x,y) is itself the saddle point [29]. The primal-dual gap G(x,y) is a measure of closeness of the primal-dual (x,y) to the primal-dual solution, and we use it to design the stopping criterion for our numerical algorithm in this paper.


## 2.3. Primal-Dual Method for Solving Saddle-Point Problem

The first-order primal-dual method summarized in [29] for convex problem was applied in this paper to solve the saddle-point problem described by (11), which is described by Algorithm 1.Algorithm 1
Step 1. Choose τ0,σ0>0 with τ0σ0L2≤1, x0∈RN, y0∈R2N, and x¯0=x0.
Step 2. Update xt, yt, x¯t, θt, τt, σt as follows:
y
t
+
1
=
arg
min
y
∈
R
2
N
〈
y
,
-
A
x
¯
t
〉
+
δ
(
y
)
+
1
2
σ
t
∥
y
-
y
t
∥
2
2
,
(
A
)
x
t
+
1
=
arg
min
x
∈
R
N
λ
2
∥
x
-
z
∥
2
2
+
〈
A
x
,
y
t
+
1
〉
+
1
2
τ
t
∥
x
-
x
t
∥
2
2
,
(
B
)
θ
t
=
1
1
+
2
γ
τ
t
,
τ
t
+
1
=
θ
t
τ
t
,
σ
t
+
1
=
σ
t
θ
t
,
x
¯
t
+
1
=
x
t
+
1
+
θ
t
(
x
t
+
1
-
x
t
)
,LetL=∥A∥ and τtσtL2≤1; then there exists a saddle-point (x^,y^) such that xt→x^ and yt→y^. Since (λ/2)∥x-z∥22 is uniformly convex such that it has a Lipschitz continuous gradient, the output (xt,yt) possesses a rate of convergence of O(1/N2) [29]. More specifically, we assume that the constant stepsizes are used; that is, τt=τ>0, σt=σ>0, and θt=θ for all t≥0. If τσL2<1; then the O(1/N) convergence rate can be guaranteed for θ=1, and O(1/N) convergence rate can be guaranteed for θ=0 [30].To summarize, our entire noise removal algorithm in a form of a pseudocode is done in following:(1)
initialization:
(i)
for a given noisy imagef of size n×n, we reord the image f row wisely into the vector z∈RN with N=n×n, z(j-1)n-i=fi,j;(ii)
compute the fractional-order discrete operatorA according to (6);(iii)
initialize the regularization parameterλ, iteration number T, and maximum permissible error Empe;(iv)
sett=0, x0=z, L=∥A∥, τ0=1/L2, σ0=1/L2, and γ=0.7λ;(2)
iteration: computex by the following steps.Step 1.
Solve the problem (A) by(16)ymt+1=ymt+σtAmx¯tmax{1,∥ymt+σtAmx¯t∥2},m=1,2,…,N.Step 2.
Substituteyt+1 into the (B) and solve the problem (B) by
(17)xt+1=xt-τtATyt+1+τtλz1+τtλ.Step 3.
Computeθt=1/1+2γτt, τt+1=τtθt, and σt+1=σt/θt.Step 4.
Updatex¯ by x¯t+1=(1+θt)xt+1-θtxt.Step 5.
Compute the primal-dual gapgt+1 according to (15) as follows:
(18)gt+1=∑m=1N‍∥Amxt+1∥2-(yt+1)TAxt+1+λ2∥xt+1-z+1λATyt+1∥22.Step 6.
If (gt+1<Empe) or (t>T), then we terminate the iteration and output xt+1; otherwise, go back to Step 1.Due to the simplicity of our model, it is easy to be implemented and can be effectively accelerated on parallel hardware such as field programmable gata array (FPGA) and graphics processing unit (GPU). In order to employ the huge computational power and the parallel processing capabilities of FPGA to obtain a fully accelerated implementation of our denoising method, the numerical method should be working on the regular grids. For this purpose, the implementation of our algorithm on FPGA is described as the following strategy. First, an image frame is captured and read in the processor from a live camera. For the captured framef∈Rn×n, the kth iteration of our method performs the following updates on prime variable u∈Rn×n and dual variable p∈Rn×n×2:
(19)pk+1=pk+σkDαu¯kmax{1,|pk+σkDαu¯k|},(20)uk+1=uk-τkdivαpk+1+τkλf1+τkλ,
where the dual variable p=(p1,p2)T corresponds with y and the discrete fractional divergence operator divαis adjoint to the discrete fractional gradient operator Dα. Note that (19) and (20), respectively, correspond with the solutions of problem (A) and (B). The updates of the other parameters are the same as the description in the foregoing pseudocode, and the number of iteration is fixed to ensure the correct timing sequence. After the iteration, the denoised image u is outputted, the next image frame is captured, and the same process is repeated until the assignment is over.


## 3. Experiments and Analysis

### 3.1. Restraint of Block Effect

The blocky effect is the main drawback produced by the typical first-order total variation denoising algorithm. In this section, some experiments are given to assess the capability of reducing block effect of our proposed fractional-order total variation (FOTV) denoising algorithm. Firstly, the one-dimensional signal is used as the test signal,y(t)=4sin(2πt)+8sin(3πt). We contaminate the given signal using the additive white Gaussian noise (AWGN) with standard deviation (SD) of 0.5. The corrupted signal is shown in Figure 1(a). The proposed FOTV denoising algorithm with α=1 and α=1.5 is used to process the contaminated signal, respectively. When α=1, our FOTV algorithm is the typical ROF denoising algorithm based on primal-dual method. In this experiment, we set λ=0.01, and the results at the 1000th iteration are, respectively, shown in Figures 1(b) and 1(c). The difference is obvious: while the first-order TV denoising algorithm approximates the observed signal with a step signal, the fractional-order algorithm with a piecewise planar signal which looks more natural and does not produce false edges.Signal processed by FOTV withα=1 and α=1.5.
(a)
Noisy signal
(b)
Processed by FOTV withα=1
(c)
Processed by FOTV withα=1.5We now consider a two-dimensional image “Lena” of size512×512. The image is degraded by using AWGN with standard deviation of 10 and the result is shown in Figure 2(a). This degraded image is, respectively, fed into the proposed FOTV denoising algorithm with α=1 and α=1.5 as initial condition, and the time evolution of the algorithms begins. The results at the 1000th iteration are, respectively, shown in Figures 2(b) and 2(c). The blocky effects are obvious in Figure 2(b), while the Figure 2(c) looks more natural and does not produce blocky effect.Image processed by FOTV withα=1 and α=1.5.
(a)
Noisy image
(b)
Processed by FOTV withα=1
(c)
Processed by FOTV withα=1.5In conclusion, the fractional-order TV algorithm can reduce blocky effect effectively comparing with the traditional first-order TV algorithm.


### 3.2. Analysis of Denoising Performance

The aim of this section is to analyze the denoising performance of the FOTV algorithm. For this purpose, three famous512×512 test images, called Barbara, Lena, and Peppers, are used in the experiments, as shown in Figure 3. In order to quantify the denoised image, we consider the peak-signal-to-noise ratio (PSNR), which has been largely used in the literature and commonly applied to determine the quality of a processed image. It can be calculated by the following formula:
(21)PSNR=10log10VL2MN∑i=1N∑j=1M(u(i,j)-u^(i,j))2,
where u is the original image, u^ is the denoised image, and VL is the maximal gray level of the image.Original image: Barbara, Lena, and Peppers.
(a)
Barbara
(b)
Lena
(c)
PeppersFirstly, in order to decide the value of the fractional-orderα, we study the relation between the PSNR and α. Figure 4 shows the relation between the PSNR and α on Lena and Peppers image corrupted by the additive white Gaussian noise with standard deviation (SD) of 20. The left image is the result of Lena image and the right one is the result of Peppers image. From the figure, we can obtain the following conclusions. Firstly, the PSNR reaches a maximum between α=1 and α=2. Secondly, the PSNR at α=1 is lower than PSNR at α>1, which is owing to the blocky effect. Finally, the PSNR decreases rapidly as α tends to zero, and when α=0, the PSNR reaches a minimum. According to these facts, we can select the fractional-order α between α=1.2 and α=2.Relation between the fractional-orderα and PSNR on Lena (a) and Peppers (b) image corrupted by the additive white Gaussian noise with standard deviation (SD) of 20.
(a)
Lena
(b)
PeppersIn order to analyze the denoising performance, the additive white Gaussian noise (AWGN) with standard deviation (SD) of10, 20, and 30 is, respectively, added in the three test images. Five denoising algorithms are used to process these given noisy images, which are the improved Perona and Malik (IP-M) model [33], fourth order (F-O) PDE model [7], improved fourth order (IF-O) PDE model [8], ROF model [1], and the proposed FOTV denoising algorithm. In this experiment, we set σ=5 and Δt=0.25 for IP-M, F-O-PDE, and IF-O-PDE models, and we set λ=0.07, τ=0.1, and K=20 for ROF model. These parameters are able to ensure the best denoising performance of the corresponding denoising model. Based on the conclusion of the foregoing experiment, we set α=1.8 for FOTV denoising algorithm, and the denoising results are shown in Table 1. In the table, the first column lists the given image and the second column lists the standard deviation of noisy image. The PSNRs of the denoised images processed by the different denoising algorithms are listed under the corresponding denoising model. The bigger the PSNR is, the better the denoising performance is. For easy observation, the biggest PSNR values are shown in boldface. From the table, it is obvious that the PSNR of our proposed algorithm is bigger than that of the other four algorithms, so we can conclude that our FOTV denoising algorithm outperforms the other algorithms.Table 1
PSNR quantitative comparison among five denoising models.
Image
SD
IP-M
F-O-PDE
IF-O-PDE
ROF
FOTV
Barbara
10
31.2427
29.3776
29.3777
31.0871
31.2850
20
26.6231
24.8155
24.8156
26.8212
27.0976
30
24.4496
22.9267
22.9273
24.7415
24.9479
Lena
10
33.6393
31.5440
31.5442
33.8378
34.4747
20
29.7665
27.8080
27.8098
30.4101
31.2153
30
27.4437
26.0883
26.1019
28.5889
29.3855
Peppers
10
33.6967
31.8199
31.8205
33.8715
34.0619
20
30.0275
28.0437
28.0457
30.1768
31.0285
30
27.4982
26.1613
26.1694
28.2689
29.0046In order to further verify the denoising performance of our proposed algorithm, Figure5 shows the denoised images processed by four denoising algorithms. The first figure is the famous image “Barbara” corrupted by the additive white Gaussian noise with standard deviation (SD) of 20. In order to show the comparison clearly, the partial enlarged view of the noisy image is shown in the second figure. The denoised result of IP-M algorithm is shown in the right-top figure. The second row shows the denoised images processed by IF-O PDE algorithm, ROF algorithm, and our FOTV algorithm. In Figures 5(c) and 5(d), although the detailed information is preserved, there are a lot of noises unremoved. The result in Figure 5(e) looks blocky. Only the result of our model, as shown in Figure 5(f), looks natural and does not produce false edges. So we conclude that our proposed FOTV algorithm is able to achieve a better tradeoff between edge preservation and noise removal.Comparison of four algorithms for the noisy image with additive Gaussian white noise ofσ=20.
(a)
Barbara with 20 SD Gaussian noise
(b)
Partial view of (a)
(c)
IP-M algorithm
(d)
IF-O PDE algorithm
(e)
ROF algorithm
(f)
FOTV algorithm


### 3.3. The Convergence and the Rate of Convergence

The convergence and the rate of convergence are two important factors for evaluating the performance of the denoising method. In this section, some experiments are given to show the convergence and convergence rate of our proposed fractional-order total variation (FOTV) denoising algorithm.For this purpose, the proposed FOTV algorithm is, respectively, used to process the image “Lena” corrupted by the additive white Gaussian noise with standard deviation (SD) of 10, 20, and 30, and the PSNR is recorded at each iteration. Figure6(a) shows the relation between the iteration and PSNR. The horizontal axis of this figure is iteration number and the vertical axis is PSNR. The different colorful curve denotes the result obtained from the different noisy image. The red dotted line is the result of the noisy image with σ=30, the green dashed line is the result of the noisy image with σ=20, and the red real line is the result of the noisy image with σ=10. From the figure, we can see that the PSNR values are stable after 50 iterations, which indicate that our proposed FOTV denoising algorithm is convergent. In addition, it can be seen that the bigger the variance of white noise is, the longer the stopping times are and the lower the PSNRs are.Relation among the iteration, PSNR, and primal-dual gap on the noisy Lena image withσ=10, σ=20, and σ=30.
(a)
Relation between the iteration and PSNR
(b)
Relation between the iteration and primal-dual gapFigure6 plots the convergence of the FOTV algorithm for the above experiment together with the theoretical O(1/N2) and O(1/N4) rate. For easy observation, the horizontal axis is the log curve of iteration number and the vertical axis is the log curve of primal-dual gap calculated in Step 5. The dashed line shows the theoretical O(1/N2) rate and the dotted line shows the theoretical O(1/N4) rate. The red line shows the convergence of FOTV algorithm in the noisy image with σ=30, the green line shows the convergence of FOTV algorithm in the noisy image with σ=20, and the blue line shows the convergence of FOTV algorithm in the noisy image with σ=10. From this figure, it can be seen that the smaller the variance of white noise is, the faster convergence is. In addition, we can see that the FOTV algorithm converges at O(1/N2) rate firstly and O(1/N4) rate in the end.


## 3.1. Restraint of Block Effect

The blocky effect is the main drawback produced by the typical first-order total variation denoising algorithm. In this section, some experiments are given to assess the capability of reducing block effect of our proposed fractional-order total variation (FOTV) denoising algorithm. Firstly, the one-dimensional signal is used as the test signal,y(t)=4sin(2πt)+8sin(3πt). We contaminate the given signal using the additive white Gaussian noise (AWGN) with standard deviation (SD) of 0.5. The corrupted signal is shown in Figure 1(a). The proposed FOTV denoising algorithm with α=1 and α=1.5 is used to process the contaminated signal, respectively. When α=1, our FOTV algorithm is the typical ROF denoising algorithm based on primal-dual method. In this experiment, we set λ=0.01, and the results at the 1000th iteration are, respectively, shown in Figures 1(b) and 1(c). The difference is obvious: while the first-order TV denoising algorithm approximates the observed signal with a step signal, the fractional-order algorithm with a piecewise planar signal which looks more natural and does not produce false edges.Signal processed by FOTV withα=1 and α=1.5.
(a)
Noisy signal
(b)
Processed by FOTV withα=1
(c)
Processed by FOTV withα=1.5We now consider a two-dimensional image “Lena” of size512×512. The image is degraded by using AWGN with standard deviation of 10 and the result is shown in Figure 2(a). This degraded image is, respectively, fed into the proposed FOTV denoising algorithm with α=1 and α=1.5 as initial condition, and the time evolution of the algorithms begins. The results at the 1000th iteration are, respectively, shown in Figures 2(b) and 2(c). The blocky effects are obvious in Figure 2(b), while the Figure 2(c) looks more natural and does not produce blocky effect.Image processed by FOTV withα=1 and α=1.5.
(a)
Noisy image
(b)
Processed by FOTV withα=1
(c)
Processed by FOTV withα=1.5In conclusion, the fractional-order TV algorithm can reduce blocky effect effectively comparing with the traditional first-order TV algorithm.


## 3.2. Analysis of Denoising Performance

The aim of this section is to analyze the denoising performance of the FOTV algorithm. For this purpose, three famous512×512 test images, called Barbara, Lena, and Peppers, are used in the experiments, as shown in Figure 3. In order to quantify the denoised image, we consider the peak-signal-to-noise ratio (PSNR), which has been largely used in the literature and commonly applied to determine the quality of a processed image. It can be calculated by the following formula:
(21)PSNR=10log10VL2MN∑i=1N∑j=1M(u(i,j)-u^(i,j))2,
where u is the original image, u^ is the denoised image, and VL is the maximal gray level of the image.Original image: Barbara, Lena, and Peppers.
(a)
Barbara
(b)
Lena
(c)
PeppersFirstly, in order to decide the value of the fractional-orderα, we study the relation between the PSNR and α. Figure 4 shows the relation between the PSNR and α on Lena and Peppers image corrupted by the additive white Gaussian noise with standard deviation (SD) of 20. The left image is the result of Lena image and the right one is the result of Peppers image. From the figure, we can obtain the following conclusions. Firstly, the PSNR reaches a maximum between α=1 and α=2. Secondly, the PSNR at α=1 is lower than PSNR at α>1, which is owing to the blocky effect. Finally, the PSNR decreases rapidly as α tends to zero, and when α=0, the PSNR reaches a minimum. According to these facts, we can select the fractional-order α between α=1.2 and α=2.Relation between the fractional-orderα and PSNR on Lena (a) and Peppers (b) image corrupted by the additive white Gaussian noise with standard deviation (SD) of 20.
(a)
Lena
(b)
PeppersIn order to analyze the denoising performance, the additive white Gaussian noise (AWGN) with standard deviation (SD) of10, 20, and 30 is, respectively, added in the three test images. Five denoising algorithms are used to process these given noisy images, which are the improved Perona and Malik (IP-M) model [33], fourth order (F-O) PDE model [7], improved fourth order (IF-O) PDE model [8], ROF model [1], and the proposed FOTV denoising algorithm. In this experiment, we set σ=5 and Δt=0.25 for IP-M, F-O-PDE, and IF-O-PDE models, and we set λ=0.07, τ=0.1, and K=20 for ROF model. These parameters are able to ensure the best denoising performance of the corresponding denoising model. Based on the conclusion of the foregoing experiment, we set α=1.8 for FOTV denoising algorithm, and the denoising results are shown in Table 1. In the table, the first column lists the given image and the second column lists the standard deviation of noisy image. The PSNRs of the denoised images processed by the different denoising algorithms are listed under the corresponding denoising model. The bigger the PSNR is, the better the denoising performance is. For easy observation, the biggest PSNR values are shown in boldface. From the table, it is obvious that the PSNR of our proposed algorithm is bigger than that of the other four algorithms, so we can conclude that our FOTV denoising algorithm outperforms the other algorithms.Table 1
PSNR quantitative comparison among five denoising models.
Image
SD
IP-M
F-O-PDE
IF-O-PDE
ROF
FOTV
Barbara
10
31.2427
29.3776
29.3777
31.0871
31.2850
20
26.6231
24.8155
24.8156
26.8212
27.0976
30
24.4496
22.9267
22.9273
24.7415
24.9479
Lena
10
33.6393
31.5440
31.5442
33.8378
34.4747
20
29.7665
27.8080
27.8098
30.4101
31.2153
30
27.4437
26.0883
26.1019
28.5889
29.3855
Peppers
10
33.6967
31.8199
31.8205
33.8715
34.0619
20
30.0275
28.0437
28.0457
30.1768
31.0285
30
27.4982
26.1613
26.1694
28.2689
29.0046In order to further verify the denoising performance of our proposed algorithm, Figure5 shows the denoised images processed by four denoising algorithms. The first figure is the famous image “Barbara” corrupted by the additive white Gaussian noise with standard deviation (SD) of 20. In order to show the comparison clearly, the partial enlarged view of the noisy image is shown in the second figure. The denoised result of IP-M algorithm is shown in the right-top figure. The second row shows the denoised images processed by IF-O PDE algorithm, ROF algorithm, and our FOTV algorithm. In Figures 5(c) and 5(d), although the detailed information is preserved, there are a lot of noises unremoved. The result in Figure 5(e) looks blocky. Only the result of our model, as shown in Figure 5(f), looks natural and does not produce false edges. So we conclude that our proposed FOTV algorithm is able to achieve a better tradeoff between edge preservation and noise removal.Comparison of four algorithms for the noisy image with additive Gaussian white noise ofσ=20.
(a)
Barbara with 20 SD Gaussian noise
(b)
Partial view of (a)
(c)
IP-M algorithm
(d)
IF-O PDE algorithm
(e)
ROF algorithm
(f)
FOTV algorithm


## 3.3. The Convergence and the Rate of Convergence

The convergence and the rate of convergence are two important factors for evaluating the performance of the denoising method. In this section, some experiments are given to show the convergence and convergence rate of our proposed fractional-order total variation (FOTV) denoising algorithm.For this purpose, the proposed FOTV algorithm is, respectively, used to process the image “Lena” corrupted by the additive white Gaussian noise with standard deviation (SD) of 10, 20, and 30, and the PSNR is recorded at each iteration. Figure6(a) shows the relation between the iteration and PSNR. The horizontal axis of this figure is iteration number and the vertical axis is PSNR. The different colorful curve denotes the result obtained from the different noisy image. The red dotted line is the result of the noisy image with σ=30, the green dashed line is the result of the noisy image with σ=20, and the red real line is the result of the noisy image with σ=10. From the figure, we can see that the PSNR values are stable after 50 iterations, which indicate that our proposed FOTV denoising algorithm is convergent. In addition, it can be seen that the bigger the variance of white noise is, the longer the stopping times are and the lower the PSNRs are.Relation among the iteration, PSNR, and primal-dual gap on the noisy Lena image withσ=10, σ=20, and σ=30.
(a)
Relation between the iteration and PSNR
(b)
Relation between the iteration and primal-dual gapFigure6 plots the convergence of the FOTV algorithm for the above experiment together with the theoretical O(1/N2) and O(1/N4) rate. For easy observation, the horizontal axis is the log curve of iteration number and the vertical axis is the log curve of primal-dual gap calculated in Step 5. The dashed line shows the theoretical O(1/N2) rate and the dotted line shows the theoretical O(1/N4) rate. The red line shows the convergence of FOTV algorithm in the noisy image with σ=30, the green line shows the convergence of FOTV algorithm in the noisy image with σ=20, and the blue line shows the convergence of FOTV algorithm in the noisy image with σ=10. From this figure, it can be seen that the smaller the variance of white noise is, the faster convergence is. In addition, we can see that the FOTV algorithm converges at O(1/N2) rate firstly and O(1/N4) rate in the end.


## 4. Conclusion

A fractional-order total variation image denoising algorithm based on primal-dual method was proposed in this paper. The main contributions are as follows:(1) the saddle point and dual formulation of the proposed fractional-order total variation model are constructed in theory; (2) the primal-dual algorithm was used to solve the constructed saddle-point problem, in which O(1/N2) convergence rate can be guaranteed; (3) the experiments demonstrate that the proposed methodology is able to avoid the blocky effect, achieve state-of-the-art performance, and guarantee O(1/N2) convergence rate.



---

*Source: 585310-2013-11-21.xml*