# The Restoration of Textured Images Using Fractional-Order Regularization

Authors: Ying Fu; Xiaohua Li; Lei Liang; Yi Zhang; Jiliu Zhou
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Image restoration problem is ill-posed, so most image restoration algorithms exploit sparse prior in gradient domain to regularize it to yield high-quality results, reconstructing an image with piecewise smooth characteristics. While sparse gradient prior has good performance in noise removal and edge preservation, it also tends to remove midfrequency component such as texture. In this paper, we introduce the sparse prior in fractional-order gradient domain as texture-preserving strategy to restore textured images degraded by blur and/or noise. And we solve the unknown variables in the proposed model using method based on half-quadratic splitting by minimizing the nonconvex energy functional. Numerical experiments show our algorithm's robust outperformance.


---

## Body

## 1. Introduction

Mathematically, the image degradation is modeled as(1)
y
=
x
⊗
h
+
n
,
where x is the original latent image and y is an observed image degraded by blur and/or noise, which is produced by convolving x with a blur point-spread-function (a.k.a. kernel) h and adding zero mean Gaussian noise n. Image restoration is recovering latent image from observed image.Image restoration is ill-posed problem, so many methods introducing priors based on natural image statistics can regularize it. Total variation regularization is originally used for noise reduction [1, 2] and has also been used for image deblurring [3]. Chan and Wong [4] introduced total variational blind deconvolution method for motion blur kernel and out-of-focus kernel. Heavy-tailed natural image priors [5, 6] and hyper-Laplacian priors [7–10] were also extensively introduced. Numerous regularization approaches have been proposed too. Wang et al. [7] presented a fast total variation deconvolution algorithm to compute TV image deconvolution. Krishnan and Fergus [8] take a novel approach to the image restoration problem arising from the use of a hyper-Laplacian prior. Xu and Jia [11] developed a fast TV-l
1 deconvolution method based on half-quadratic splitting.While image reconstructed by algorithms above suppresses noise and preserve edges, it has piecewise smooth characteristic that the mid-frequency components such as textures are removed too.In digital images, the gray values between neighboring pixels have high correlation. This highly self-similar fractal information of image fractal information is usually represented by complex textural features, and the works in [12–18] showed that fractional-order gradient is more suitable to deal with fractal-like textures. It has been proved in [12] that the fractional-order derivative satisfies the lateral inhibition principle of biologic visual system better than the integer-order derivative. The fractional-order derivative operators have been used in texture enhancement [13], image denoising [14, 15], and image inpainting [16, 17]. Jun and Zhihui [14] replaced the first-order derivative in the regularized term of ROF model with the fractional-order derivative. Bai and Feng [15] designed fractional-order anisotropic diffusion equation to remove noise. Zhang et al. [16] exploited fractional-order TV sinogram inpainting model to reduce metal artifacts for X-ray computed tomography. In [18], fractional total variation method was introduced to restore textured image. This work shows that the fractional-order derivative not only nonlinearly preserves the textural details but also eliminates the staircase effect caused by low integral-order derivative in image processing. Different from work in [18], the sparse prior in fractional-order gradient domain is considered in our work, which is more suitable for the texture of image. It is explained clearly in Figures 2 and 3.This paper presents fractional-order regularization for the restoration of textured image degraded by blur and/or additive noise. R. Tony uses the Laplacian prior in fractional-order gradient domain forα
=
1 to preserve the texture. According to our analysis in the next section, hyper-Laplacian image prior in fractional-order gradient domain for0
<
α
≤
1 is more suitable to keep different texture for different texture image.The outline of this paper is as follows. In Section2, we analyze the reason why integral-order regularization fails to restore image texture. In Section 3, our fractional-order regularization model is proposed and based on half-quadratic splitting, we solve model using efficient alternating minimization method. Numerical experiments and comments are provided in Section 4 and the paper is concluded in Section 5.


## 2. Motivation

The priorp
(
x
) favors natural image, usually based on the observation that their heavy-tailed gradient distribution is sparse. For example, Figure 1 shows textured image and a histogram of its gradient magnitudes inx-direction andy-direction, respectively. The distribution shows that the image contains primarily small or zero gradients, but a few gradients have large magnitudes. A common measure [19] is
(2)
log
⁡
p
(
x
)
=
-
∑
i
‍
|
G
x
,
i
(
x
)
|
α
+
|
G
y
,
i
(
x
)
|
α
+
constant
,
where G
x
,
i
(
x
) and G
y
,
i
(
x
) denote the horizontal and vertical derivatives at pixel i (here, the simple [
-
1
1
] and [
-
1
1
]
T filters are used) and exponent value α
∈
(
0,2
]. α
<
1 leads to sparse prior and natural images usually correspond to α in the range of [
0.5
,
0.8
] [19]. α
=
1 and α
=
2 are Laplacian prior and Gaussian prior, respectively.(a) A textured image. (b) Thex-direction log distribution of gradient magnitudes. (c) The y-direction log distribution of gradient magnitudes.
(a)
(b)
(c)Analysis of restoration on 1D signal using gradient prior. ((a) and (c)) Sharp and blurred signal; ((b) and (d)) sum of gradients-
log
⁡
p
(
x
)
=
∑
i
‍
|
G
i
(
x
)
|
α as a function of α.
(a)
(b)
(c)
(d)Analysis of restoration on 1D signal: (a) sharp versus blurred signal; (b) sum of gradients-
log
⁡
p
(
x
)
=
∑
i
‍
|
G
i
(
x
)
|
α as a function of α; (c) sum of gradients -
log
⁡
p
(
x
)
=
∑
i
‍
|
G
i
v
(
x
)
|
α as a function of α (the value of v is 0.3).
(a)
(b)
(c)The image restoration methods use the sparse prior term as a regularized term of variational energy functional [19], which is
(3)
min
⁡
x
λ
(
x
⊗
h
-
y
)
+
∑
i
=
1
I
‍
(
|
G
x
,
i
x
|
α
+
|
G
y
,
i
x
|
α
)
.The failure of restoring texture with the sparse gradient prior depends on the fact that the value of energy does not always decrease during restoration process, so the no-blur explanation is usually favored. To understand this, consider the 1D signals in Figure2.For sharp edge in Figure2(a), while Gaussian prior favors the blurry explanation, the sparse prior (
α
<
1
) favors the correct sharp explanation in Figure 2(b). The signal considered in Figure 2(c) illustrates that natural image contains a lot of medium contrast textures, which dominate the statistics more than step edges. As a result, blurring natural image reduces the overall contrast which cannot be restored by Gaussian prior or even sparse priors as in Figure 2(d).The reason is that the gradient profile in fractal-like textures is close to Gaussian distribution and these small values are severely penalized by the sparse gradient prior.A fractional-order gradient log distribution can be expressed as follows [18]:
(4)
log
⁡
p
(
x
)
=
-
∑
i
‍
|
G
x
,
i
v
x
|
α
+
|
G
y
,
i
v
x
|
α
+
constant,
where G
x
,
i
v
x and G
y
,
i
v
x denote the horizontal and vertical fractional-order derivatives at pixel i and v is the fractional order v
∈
(
0,4
]. The exponent value is the same as α value in (2).Compared with result in Figure3(b), the sharp explanation in Figure 3(c) is favored by sparse prior even by Gaussian prior in fractional-order gradient domain.


## 3. The Proposed Model and Algorithm

The corresponding energy functional is as follows [18]:
(5)
min
⁡
x
λ
(
x
⊗
h
-
y
)
+
∑
i
=
1
‍
(
|
G
x
,
i
v
x
|
α
+
|
G
y
,
i
v
x
|
α
)
,
where i is the pixel index and ⊗ is the 2-dimensional convolution operator, and a weighting term λ
=
3
e
3 controls the strength of the regularization. G
x
,
i
v
x and G
y
,
i
v
x denote the horizontal and vertical fractional-order derivatives at pixel i defined by our coauthor as Tables 1 and 2 [16].Table 1
Operator ofx-direction: G
x
v
x.
(a)
⋮
⋮
⋮
⋮
⋮
⋮
⋮
0
0
0
⋯
0
⋮
0
0
0
C
s
-
1
C
s
0
C
s
1
⋯
C
s
k
⋯
C
s
n
-
2
C
s
n
-
1
C
s
n
0
0
0
⋯
0
⋮
0
0
0
⋮
⋮
⋮
⋮
⋮
⋮
⋮
(b)
⋮
⋮
⋮
⋮
⋮
⋮
⋮
0
0
0
⋯
0
⋮
0
0
0
C
s
n
C
s
n
-
1
C
s
n
-
2
⋯
C
s
k
⋯
C
s
1
C
s
0
C
s
-
1
0
0
0
⋯
0
⋮
0
0
0
⋮
⋮
⋮
⋮
⋮
⋮
⋮
7Table 2
Operator ofy-direction: G
y
v
x.
(a)
⋯
0
C
s
n
0
⋯
⋯
0
C
s
n
-
1
0
⋯
⋯
0
C
s
n
-
2
0
⋯
⋮
⋮
⋮
⋯
0
C
s
n
0
⋯
⋮
⋮
⋮
⋯
0
C
s
1
0
⋯
⋯
0
C
s
0
0
⋯
⋯
0
C
s
-
1
0
⋯
(b)
⋯
0
C
s
-
1
0
⋯
⋯
0
C
s
0
0
⋯
⋯
0
C
s
1
0
⋯
⋮
⋮
⋮
⋯
0
C
s
k
0
⋯
⋮
⋮
⋮
⋯
0
C
s
n
-
2
0
⋯
⋯
0
C
s
n
-
1
0
⋯
⋯
0
C
s
n
0
⋯The coefficients of the operator in Tables1 and 2 are
(6)
C
s
-
1
=
v
4
+
v
2
8
C
s
0
=
1
-
v
2
2
-
v
3
8
C
s
1
=
-
5
v
4
-
5
v
3
16
+
v
4
16
⋮
C
s
k
=
1
Γ
(
-
v
)
[
Γ
(
k
-
v
-
1
)
(
k
+
1
)
!
·
(
v
4
+
v
2
8
)
+
Γ
(
k
-
v
)
k
!
·
(
1
-
v
2
4
)
+
Γ
(
k
-
v
-
1
)
(
k
-
1
)
!
·
(
-
v
4
+
v
2
8
)
]
,
⋮
C
s
n
-
2
=
1
Γ
(
-
v
)
[
Γ
(
n
-
v
-
1
)
(
n
-
1
)
!
·
(
v
4
+
v
2
8
)
+
Γ
(
n
-
v
-
2
)
(
n
-
2
)
!
·
(
1
-
v
2
4
)
+
Γ
(
n
-
v
-
3
)
(
n
-
3
)
!
·
(
-
v
4
+
v
2
8
)
]
,
⋮
C
s
n
-
2
=
Γ
(
n
-
v
-
1
)
(
n
+
1
)
!
Γ
(
-
v
)
·
(
1
-
v
2
8
)
+
Γ
(
n
-
v
-
2
)
(
n
-
2
)
!
Γ
(
-
v
)
·
(
-
v
4
+
v
2
8
)
C
s
n
=
Γ
(
n
-
v
-
1
)
(
n
-
1
)
!
Γ
-
v
·
(
-
v
4
+
v
2
8
)
.
Equation (5) contains nonlinear penalties for regularization term, so we propose alternating minimization (AM) method, based on a half-quadratic splitting to solve it [18, 20]. We introduce auxiliary variables u and w
=
(
w
x
,
w
y
) at each pixel, so the energy functional in (5) can be modified as
(7)
min
⁡
x
,
w
θ
2
(
x
⊗
h
-
y
)
2
+
λ
|
u
|
+
∑
i
=
1
‍
(
β
2
(
∥
G
x
,
i
v
x
-
w
x
,
i
∥
2
2
+
∥
G
y
,
i
v
x
-
w
y
,
i
∥
2
2
)
+
|
w
x
,
i
|
α
+
|
w
y
,
i
|
α
β
2
)
,
where the first two terms are used to ensure the similarity between the measures and the corresponding auxiliary variables. As β
→
∞ and θ
→
∞ the solution of (6) converges to that of (5). Equation (7) can be solved by AM method through fixing other variables to solve x
,
w
, and u independently.


### 3.1.<?bold?>
x Subsolution

Given fixed values ofu and w from the previous iteration, (7) is quadratic in x. So we compute x by minimizing
(8)
E
(
x
;
u
,
w
)
=
∥
x
⊗
h
-
y
-
u
∥
2
+
∑
i
=
1
‍
(
β
θ
(
∥
G
x
,
i
v
x
-
w
x
,
i
∥
2
2
+
∥
G
y
,
i
v
x
-
w
y
,
i
∥
2
2
)
)
.
The optimal x is
(9)
(
θ
β
H
T
H
+
G
x
v
T
G
x
v
+
G
y
v
T
G
y
v
)
x
=
G
x
v
T
w
x
+
G
y
v
T
w
y
+
θ
β
H
T
(
y
+
u
)
,
where H
x
=
h
⊗
x. According to Parseval’s theorem after the Fourier transform, (8) has the closed form solution in minimization, which enables us to find the optimal x directly:
(10)
x
=
F
-
1
(
+
θ
β
F
(
H
)
*
F
(
H
)
)
-
1
(
∑
F
(
G
x
v
)
*
F
(
w
x
)
+
F
(
G
y
v
)
*
F
(
w
y
)
+
θ
β
F
(
H
)
*
F
(
y
+
u
)
∑
)
×
(
∑
F
(
G
x
v
)
*
F
(
G
x
v
)
+
F
(
G
y
v
)
*
F
(
G
y
v
)
+
θ
β
F
(
H
)
*
F
(
H
)
∑
)
-
1
)
,
where F
(
·
) and F
(
·
)
-
1 denote the fast Fourier transform and inverse fast Fourier transform, respectively. * is the complex conjugate operator.


### 3.2.u Subsolution

Here,u and w belong to different terms. They are not coupled with each other in the functional, so their optimization is independent. Given fixed value of x, we compute u by minimizing
(11)
E
(
u
;
x
)
=
1
2
∥
u
-
(
x
⊗
h
-
y
)
∥
+
λ
θ
|
u
|
.
According to shrinkage formula [21], the optimal u is
(12)
u
=
sign
⁡
(
x
⊗
h
-
y
)
max
⁡
(
∥
x
⊗
h
-
y
∥
-
λ
θ
,
0
)
.


### 3.3.w Subsolution

We have the following:(13)
E
˙
(
w
x
;
x
)
=
|
w
x
|
α
+
β
2
(
∥
G
x
v
x
-
w
x
∥
2
2
)
,
E
(
w
y
;
x
)
=
|
w
y
|
α
+
β
2
(
∥
G
y
v
x
-
w
y
∥
2
2
)
.Forα
=
2 case, the subproblem about w
x and w
y is quadratic.Forα
=
1 case, the optimal solution for w
x and w
y can be derived by shrinkage formula too:
(14)
w
x
˙
=
sign
⁡
(
G
x
v
x
)
*
max
⁡
(
abs
(
G
x
v
x
)
-
1
β
,
0
)
,
w
y
=
sign
⁡
(
G
y
v
x
)
*
max
⁡
(
abs
(
G
y
v
x
)
-
1
β
,
0
)
.For the otherα case, setting the derivative of (13) with regard to w
x and w
y to zero gives
(15)
α
˙
|
w
x
|
α
-
1
sign
⁡
(
w
x
)
+
β
(
w
x
-
G
x
v
x
)
=
0
,
α
|
w
y
|
α
-
1
sign
⁡
(
w
y
)
+
β
(
w
y
-
G
y
v
x
)
=
0
.
Two special α cases are discussed here.Forα
=
1
/
2 case, about w
x, (15) becomes
(16)
1
2
|
w
x
|
-
1
/
2
sign
⁡
(
w
x
)
+
β
(
w
x
-
G
x
v
x
)
=
0
,
(17)
w
x
3
-
2
(
G
x
v
x
)
w
x
2
+
(
G
x
v
x
)
2
w
x
-
sign
⁡
(
w
x
)
4
β
2
=
0
.BecauseG
x
v
x is fixed and w
x lies between 0 and G
x
v
x, we can replace sign
⁡
(
w
x
) with sign
⁡
(
G
x
v
x
). Equation (17) can be rewritten as
(18)
w
x
3
-
2
(
G
x
v
x
)
w
x
2
+
(
G
x
v
x
)
2
w
x
-
sign
⁡
(
G
x
v
x
)
4
β
2
=
0
.So we can get the cubic polynomials aboutw
x and w
y:
(19)
w
˙
x
3
-
2
(
G
x
v
x
)
w
x
2
+
(
G
x
v
x
)
2
w
x
-
sign
⁡
(
G
x
v
x
)
4
β
2
=
0
,
w
y
3
-
2
(
G
y
v
x
)
w
y
2
+
(
G
y
v
x
)
2
w
y
-
sign
⁡
(
G
y
v
x
)
4
β
2
=
0
.The value ofw
x and w
y is either 0 or the root of cubic polynomial in (19).Forα
=
2
/
3 case, we can get the quartic polynomials about w
x and w
y:
(20)
w
˙
x
4
-
3
(
G
x
v
x
)
w
x
3
+
3
(
G
x
v
x
)
2
w
x
2
-
(
G
x
v
x
)
3
w
x
+
8
27
β
3
=
0
,
w
y
4
-
3
(
G
y
v
x
)
w
y
3
+
3
(
G
y
v
x
)
2
w
y
2
-
(
G
y
v
x
)
3
w
y
+
8
27
β
3
=
0
.The value ofw
x and w
y is either 0 or the root of cubic polynomial in (20).Given the roots of cubic and quartic polynomials and zero solution, we need to determine which one corresponds to the global minima of (12), which can be confirmed by the following scheme.Letr be the nonzero real root. If 0 value is the optimum solution of (13), for G
x
v
x
, this implies
(21)
|
r
|
α
+
β
2
(
r
-
G
x
v
x
)
2
>
β
(
G
x
v
x
)
2
2
,
(22)
sign
⁡
(
r
)
|
r
|
(
α
-
1
)
+
β
2
(
r
-
G
x
v
x
)
≶
0
,
r
≶
0
.
We can use (15) to eliminate sign
⁡
(
r
)
|
r
|
(
α
-
1
) from (15) and (22), yielding the condition
(23)
r
≶
2
G
x
v
x
α
-
1
α
-
2
,
since sign
⁡
(
r
)
=
sign
⁡
(
G
x
v
x
). So w
x
=
r if r is between 2
G
x
v
x
/
3 and G
x
v
x in the α
=
1
/
2 case or between G
x
v
x
/
2 and G
x
v
x in the α
=
2
/
3 case. Otherwise, w
x
=
0. The same scheme applies to G
y
v
x.For otherα cases, w
x and w
y can be computed by Newton method.


### 3.4. Algorithm

β and θ are positive values to enforce the similarity between the auxiliary variables and the respective terms. We empirically set β
0
=
1, β
Max
⁡
=
256 and θ
0
=
1, θ
Max
⁡
=
λ.The algorithm of this fractional-order regularization model is shown in Algorithm1.Algorithm 1:Fractional-order regularization.
Input: observed image y, PSF h, penalty parameters λ, exponent α and fractional-order v
Input: β
0, β
Max
⁡, θ
0, θ
Max
⁡
Initializex
=
y, θ
=
θ
0, β
=
β
0
whileθ
<
θ
Max
⁡, do
solve for u using (12)
whileβ
<
β
Max
⁡ do
Given x, solve for w according to our discussion
Given w, solve for x using (10)
β
=
2
β
end while
θ
=
2
θ
end while
Output: Estimated image x


## 3.1.<?bold?>
x Subsolution

Given fixed values ofu and w from the previous iteration, (7) is quadratic in x. So we compute x by minimizing
(8)
E
(
x
;
u
,
w
)
=
∥
x
⊗
h
-
y
-
u
∥
2
+
∑
i
=
1
‍
(
β
θ
(
∥
G
x
,
i
v
x
-
w
x
,
i
∥
2
2
+
∥
G
y
,
i
v
x
-
w
y
,
i
∥
2
2
)
)
.
The optimal x is
(9)
(
θ
β
H
T
H
+
G
x
v
T
G
x
v
+
G
y
v
T
G
y
v
)
x
=
G
x
v
T
w
x
+
G
y
v
T
w
y
+
θ
β
H
T
(
y
+
u
)
,
where H
x
=
h
⊗
x. According to Parseval’s theorem after the Fourier transform, (8) has the closed form solution in minimization, which enables us to find the optimal x directly:
(10)
x
=
F
-
1
(
+
θ
β
F
(
H
)
*
F
(
H
)
)
-
1
(
∑
F
(
G
x
v
)
*
F
(
w
x
)
+
F
(
G
y
v
)
*
F
(
w
y
)
+
θ
β
F
(
H
)
*
F
(
y
+
u
)
∑
)
×
(
∑
F
(
G
x
v
)
*
F
(
G
x
v
)
+
F
(
G
y
v
)
*
F
(
G
y
v
)
+
θ
β
F
(
H
)
*
F
(
H
)
∑
)
-
1
)
,
where F
(
·
) and F
(
·
)
-
1 denote the fast Fourier transform and inverse fast Fourier transform, respectively. * is the complex conjugate operator.


## 3.2.u Subsolution

Here,u and w belong to different terms. They are not coupled with each other in the functional, so their optimization is independent. Given fixed value of x, we compute u by minimizing
(11)
E
(
u
;
x
)
=
1
2
∥
u
-
(
x
⊗
h
-
y
)
∥
+
λ
θ
|
u
|
.
According to shrinkage formula [21], the optimal u is
(12)
u
=
sign
⁡
(
x
⊗
h
-
y
)
max
⁡
(
∥
x
⊗
h
-
y
∥
-
λ
θ
,
0
)
.


## 3.3.w Subsolution

We have the following:(13)
E
˙
(
w
x
;
x
)
=
|
w
x
|
α
+
β
2
(
∥
G
x
v
x
-
w
x
∥
2
2
)
,
E
(
w
y
;
x
)
=
|
w
y
|
α
+
β
2
(
∥
G
y
v
x
-
w
y
∥
2
2
)
.Forα
=
2 case, the subproblem about w
x and w
y is quadratic.Forα
=
1 case, the optimal solution for w
x and w
y can be derived by shrinkage formula too:
(14)
w
x
˙
=
sign
⁡
(
G
x
v
x
)
*
max
⁡
(
abs
(
G
x
v
x
)
-
1
β
,
0
)
,
w
y
=
sign
⁡
(
G
y
v
x
)
*
max
⁡
(
abs
(
G
y
v
x
)
-
1
β
,
0
)
.For the otherα case, setting the derivative of (13) with regard to w
x and w
y to zero gives
(15)
α
˙
|
w
x
|
α
-
1
sign
⁡
(
w
x
)
+
β
(
w
x
-
G
x
v
x
)
=
0
,
α
|
w
y
|
α
-
1
sign
⁡
(
w
y
)
+
β
(
w
y
-
G
y
v
x
)
=
0
.
Two special α cases are discussed here.Forα
=
1
/
2 case, about w
x, (15) becomes
(16)
1
2
|
w
x
|
-
1
/
2
sign
⁡
(
w
x
)
+
β
(
w
x
-
G
x
v
x
)
=
0
,
(17)
w
x
3
-
2
(
G
x
v
x
)
w
x
2
+
(
G
x
v
x
)
2
w
x
-
sign
⁡
(
w
x
)
4
β
2
=
0
.BecauseG
x
v
x is fixed and w
x lies between 0 and G
x
v
x, we can replace sign
⁡
(
w
x
) with sign
⁡
(
G
x
v
x
). Equation (17) can be rewritten as
(18)
w
x
3
-
2
(
G
x
v
x
)
w
x
2
+
(
G
x
v
x
)
2
w
x
-
sign
⁡
(
G
x
v
x
)
4
β
2
=
0
.So we can get the cubic polynomials aboutw
x and w
y:
(19)
w
˙
x
3
-
2
(
G
x
v
x
)
w
x
2
+
(
G
x
v
x
)
2
w
x
-
sign
⁡
(
G
x
v
x
)
4
β
2
=
0
,
w
y
3
-
2
(
G
y
v
x
)
w
y
2
+
(
G
y
v
x
)
2
w
y
-
sign
⁡
(
G
y
v
x
)
4
β
2
=
0
.The value ofw
x and w
y is either 0 or the root of cubic polynomial in (19).Forα
=
2
/
3 case, we can get the quartic polynomials about w
x and w
y:
(20)
w
˙
x
4
-
3
(
G
x
v
x
)
w
x
3
+
3
(
G
x
v
x
)
2
w
x
2
-
(
G
x
v
x
)
3
w
x
+
8
27
β
3
=
0
,
w
y
4
-
3
(
G
y
v
x
)
w
y
3
+
3
(
G
y
v
x
)
2
w
y
2
-
(
G
y
v
x
)
3
w
y
+
8
27
β
3
=
0
.The value ofw
x and w
y is either 0 or the root of cubic polynomial in (20).Given the roots of cubic and quartic polynomials and zero solution, we need to determine which one corresponds to the global minima of (12), which can be confirmed by the following scheme.Letr be the nonzero real root. If 0 value is the optimum solution of (13), for G
x
v
x
, this implies
(21)
|
r
|
α
+
β
2
(
r
-
G
x
v
x
)
2
>
β
(
G
x
v
x
)
2
2
,
(22)
sign
⁡
(
r
)
|
r
|
(
α
-
1
)
+
β
2
(
r
-
G
x
v
x
)
≶
0
,
r
≶
0
.
We can use (15) to eliminate sign
⁡
(
r
)
|
r
|
(
α
-
1
) from (15) and (22), yielding the condition
(23)
r
≶
2
G
x
v
x
α
-
1
α
-
2
,
since sign
⁡
(
r
)
=
sign
⁡
(
G
x
v
x
). So w
x
=
r if r is between 2
G
x
v
x
/
3 and G
x
v
x in the α
=
1
/
2 case or between G
x
v
x
/
2 and G
x
v
x in the α
=
2
/
3 case. Otherwise, w
x
=
0. The same scheme applies to G
y
v
x.For otherα cases, w
x and w
y can be computed by Newton method.


## 3.4. Algorithm

β and θ are positive values to enforce the similarity between the auxiliary variables and the respective terms. We empirically set β
0
=
1, β
Max
⁡
=
256 and θ
0
=
1, θ
Max
⁡
=
λ.The algorithm of this fractional-order regularization model is shown in Algorithm1.Algorithm 1:Fractional-order regularization.
Input: observed image y, PSF h, penalty parameters λ, exponent α and fractional-order v
Input: β
0, β
Max
⁡, θ
0, θ
Max
⁡
Initializex
=
y, θ
=
θ
0, β
=
β
0
whileθ
<
θ
Max
⁡, do
solve for u using (12)
whileβ
<
β
Max
⁡ do
Given x, solve for w according to our discussion
Given w, solve for x using (10)
β
=
2
β
end while
θ
=
2
θ
end while
Output: Estimated image x


## 4. Numerical Experiments

We consider the restoration of a blur- and noise-contaminated test image represented by255
×
255 pixels. In order to compare the accuracy of FOR (fractional order regularization) and IOR (integer order regularization) more precisely, we list in Table 3 the peak signal-to-noise ratio (PSNR) and gray-scale structural similarity (SSIM) as quality metric. PSNR is most easily defined via the mean squared error (MSE). Given a noise-free M by N image I and its noisy approximation I
^, MSE is defined as
(24)
MSE
=
1
M
N
∑
m
=
1
,
n
=
1
M
,
N
‍
[
I
(
m
,
n
)
-
I
^
(
m
,
n
)
]
.
PSNR is defined as
(25)
PSNR
=
1
10
log
⁡
(
25
5
2
MSE
)
(
dB
)
,
and SSIM is defined as
(26)
SSIM
=
(
2
μ
I
μ
K
+
c
1
)
(
2
σ
I
K
+
c
2
)
(
μ
I
2
+
μ
K
2
+
c
1
)
(
σ
I
2
+
σ
K
2
+
c
2
)
,
where I and K are different images, μ
I and μ
K are the average of I and K, respectively, σ
I
2 and σ
K
2 are the variance of I and K, respectively, and σ
I
K is the covariance of I and K. c
1 and c
2 are constant.Table 3
PSNR and SSIM of image restoration by IOR and FOR.
Image
α
PSNRIOR
SSIMIOR
v
α
PSNRFOR
SSIMFOR
Barbara (256 * 256)
0.8
—
0.8411
0.4
0.3
—
0.8580
Bubble (512 * 512)
0.4
2.5249
0.6949
1.8
0.4
2.7343
0.6987The desired blur- and noise-free image is depicted in Figure4. The image is contaminated by motion blur generated by Matlab function (fspecial(“motion”,10,20)). The resulting image is displayed in Figure 4(f). In Table 3 the second column, with header PSNR and SSIM values for images that have been corrupted by motion blur, is characterized by v
=
0.4 and α
=
0.3.Deblurring: (a) clear image; (b) synthesized blurred image with PSF (fspecial(“motion”,10,20)); (c) image restoration by Lucy-Richardson algorithm; (d) image restoration by IOR; (e) image restoration by [18]; (f) image restoration by FOR (g) closeups of (d), (e), and (f).
(a)
(b)
(c)
(d)
(e)
(f)
(g)The desired blur- and noise-contaminated image is depicted in Figure5. The image is contaminated by motion blur, adding white Gaussian noise (its standard variance is 0.003). The resulting image is displayed in Figure 5(d). In Table 3 the third column, with header PSNR and SSIM values for images that have been corrupted by motion blur, is characterized by v
=
1.8 and α
=
0.4.Deblurring and denoising: (a) clear image; (b) synthesized blurred image and adding white Gaussian noise (its standard variance is0.003); (c) image restoration by IOR; (d) image restoration by FOR.
(a)
(b)
(c)
(d)Figure6 shows the result of deconvolving a real blurry image. We estimate the blur kernel using the algorithm in [9]. Again, textured regions are better reconstructed using our method in visual quality. Figure 6(b) is restored by total variation. Figure 6(c) is restored by fractional-order total variation. Figure 6(d) shows the details in Figure 6(b) and Figure 6(c).Testing our algorithm with real-life blurry images. (a) Blurry image. (b) Restored image by using the algorithm in [9]. (c) Restored image by our algorithm. (d) Comparison of details of image window. Left: details in (a), middle: details in (b), and right: details in (c).
(a)
(b)
(c)
(d)


## 5. Conclusion

By introducing sparse prior in fractional-order gradient domain, we propose a fractional-order regularization method for the restoration of textured image degraded by blur and/or noise. The regularizer is constructed by using fractional-order derivatives, where the choice of the fractional-order is driven by different textured image. This makes the proposed model an efficient tool to preserve texture well. Numerical results show that the proposed model yields better SSIM and PSNR value and visual effects than using integral-order regularization method.Our following work is to use an automatic texture detection procedure for textured image restoration. Different parameters are applied for different textures.



---

*Source: 356906-2014-04-14.xml*