# Mixed Higher Order Variational Model for Image Recovery

Authors: Pengfei Liu; Liang Xiao; Liancun Xiu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A novel mixed higher order regularizer involving the first and second degree image derivatives is proposed in this paper. Using spectral decomposition, we reformulate the new regularizer as a weightedL
1-L
2 mixed norm of image derivatives. Due to the equivalent formulation of the proposed regularizer, an efficient fast projected gradient algorithm combined with monotone fast iterative shrinkage thresholding, called, FPG-MFISTA, is designed to solve the resulting variational image recovery problems under majorization-minimization framework. Finally, we demonstrate the effectiveness of the proposed regularization scheme by the experimental comparisons with total variation (TV) scheme, nonlocal TV scheme, and current second degree methods. Specifically, the proposed approach achieves better results than related state-of-the-art methods in terms of peak signal to ratio (PSNR) and restoration quality.


---

## Body

## 1. Introduction

During image acquisition, the observed images are mainly degraded by blurring and random noise, which are tightly linked to imaging system. In order to attenuate the undesired effects, image recovery plays an integral part of modern imaging sciences in such areas, including medical imaging, microscopy, astronomical imaging, remote sensing, and photography [1–5].Image recovery which amounts to estimating a desired image from its degraded measurements is mathematically an inverse problem [6]. To cope with the ill posedness of image recovery, the standard approach is to use image a priori to constrain the solutions so that many variational regularization approaches have been successfully proposed for image recovery. As we know, one of the most popular regularizers is the total variation (TV) seminorm which was firstly introduced by Rudin et al. for image denoising [7] and turned out to be very efficient for preserving image edges. Thus, TV has been widely and effectively used in many other image applications [8–11] and provides state-of-the-art result. However, TV-based reconstruction methods often give rise to the undesired staircase effect, since TV favors piecewise constant solutions. Thus, the main challenge for image recovery is to reduce the staircase effect while preserving edges sharpness as far as possible. As is known to us, piecewise-vanishing second degree derivatives yield piecewise-linear solutions which better fit smooth intensity changes. To alleviate the staircase effect, several regularizers involving second degree image derivatives were introduced mainly for image denoising over past few years [12–15] and provide better performance than TV. Recently, many novel higher order regularization methods which constitute valid extensions of TV are effectively introduced for more image processing tasks [16–18] so that they minimize the staircase effect and provide state-of-the-art results. However, the second-order regularizers [16–18] which are the mixed norms of second degree directional image derivatives tend to suppress larger second degree directional derivatives thus leading to blurring across the sharp edges. Thus, the reconstructed images with better preserved edges and less staircase effect are more desirable. To this end, as reported in [19–23], the regularizers that combine first and second degree image derivatives have been proposed so as to provide better edge preservation while reducing the staircase effect in smooth regions of image; thus, we are motivated to investigate the regularizers involving derivatives of different degrees in this paper.In this paper, our main contributions are twofolds. First, we reinterpret the isotropic second degree total variation (ISDTV) regularizer [16] which is theL
1-L
2mixed norms of the second degree directional derivatives and give its equivalent formulation. We then propose a novel mixed higher order regularizer involving the first and second degree derivatives of the image and reformulate the new regularizer as a weightedL
1-L
2mixed norm of image derivatives in spectral decomposition framework. The new regularizer which inherits the desirable properties of TV turns out to be well suited for image recovery. Second, for the proposed variational image recovery model, a fast projected gradient algorithm combined with monotone fast iterative shrinkage thresholding (FPG-MFISTA) which is derived from majorization-minimization (MM) framework has been utilized to solve the minimization problem under additional convex constraints. Finally, we compare our proposed schemes with TV [24], nonlocal TV [25], and current second degree methods [18, 26], which currently provide state-of-the-art results.The rest of the paper is organized as follows. In Section2, we first present the mathematic formulation of image recovery problem and review the ISDTV regularizer. We then propose a novel regularizer at the end of Section 2. In Section 3, we propose our variational regularization-based image recovery model and give a detailed description of the corresponding minimization algorithm. In Section 4, we assess the performance of our approach on both natural and biomedical images and provide systematic comparisons with TV, nonlocal TV, and current second degree methods. Finally, we conclude the paper in Section 5.


## 2. Mixed Higher Order Total Variation for Image Recovery

### 2.1. Regularization for Image Recovery Problem

We first consider the recovery of a continuously differentiable imagef
:
Ω
→
ℜfrom its degraded observationg, whereΩ
⊂
ℜ
2 is the spatial support of the image. The most commonly used image observation model can be linearly formulated as
(1)
g
=
A
f
+
e
,
whereAis a linear blur operator andeis Gaussian white noise with standard deviationσ. As we know, the recovery offfrom its observationgis an ill-posed problem, due to the fact that the blur operatorAis ill-conditioned. To stabilize the recovery process, some available a priori knowledge about the desired image is utilized. Thus, a common variational regularization approach is employed to formulate image recovery as an optimization problem:
(2)
f
^
=
argmin
f
{
C
(
f
)
=
1
2
∫
Ω
|
g
-
A
f
|
2
d
x
d
y
+
λ
J
(
f
)
}
,
where the first term is known as data fidelity term, which measures the consistency between estimation and observation, while the second one is called regularization term which provides some a priori knowledge of the desired image, and the regularization parameterλ
≥
0provides a balance between the two terms.


### 2.2. Isotropic Second Degree Total Variation and Equivalent Formulation

Hu and Jacob have introduced the isotropic higher degree total variation (I-HDTV) scheme [16], which inherits the desirable properties of TV, such as, convexity and invariance to rotations and translations. As detailed in [16], the first degree directional derivative offalong the unit vectors
0
(
θ
)
=
(
cos
⁡
θ
,
sin
θ
)
Tat coordinate(
x
,
y
)denoted byf
θ
,
1
(
x
,
y
)is defined as
(3)
f
θ
,
1
(
x
,
y
)
=
[
cos
⁡
θ
,
sin
θ
]
︸
s
0
T
(
θ
)
[
∂
f
(
x
,
y
)
∂
x
∂
f
(
x
,
y
)
∂
y
]
︸
g
0
(
x
,
y
)
,
whereg
0
(
x
,
y
)denotes the gradient offat coordinate(
x
,
y
)and(
·
)
T denotes transpose operator. Specially, as proved in [16], the standard TV is reinterpreted as theL
1-L
2mixed norms of the first degree directional derivative, formulated as
(4)
TV
(
f
)
=
∫
Ω
∥
g
0
(
x
,
y
)
∥
2
d
x
d
y
=
2
∫
Ω
∥
f
θ
,
1
(
x
,
y
)
∥
L
2
[
0,2
π
]
d
x
d
y
=
2
∫
Ω
1
2
π
∫
0
2
π
|
f
θ
,
1
(
x
,
y
)
|
2
d
θ
d
x
d
y
,
where∥
·
∥
2denotes the Euclidean norm.Now, we consider the isotropic second degree total variation (ISDTV) of imagef, which is defined as
(5)
J
2
(
f
)
=
∫
Ω
∥
f
θ
,
2
(
x
,
y
)
∥
L
2
[
0,2
π
]
d
x
d
y
=
∫
Ω
1
2
π
∫
0
2
π
|
f
θ
,
2
(
x
,
y
)
|
2
d
θ
d
x
d
y
,
wheref
θ
,
2
(
x
,
y
)is the second degree directional derivative of imagefat coordinate(
x
,
y
), defined as
(6)
f
θ
,
2
(
x
,
y
)
=
s
T
(
θ
)
g
(
x
,
y
)
,
with
(7)
s
(
θ
)
=
(
cos
⁡
2
θ
,
2
sin
θ
cos
⁡
θ
,
si
n
2
θ
)
T
,
(8)
g
(
x
,
y
)
=
(
∂
2
f
(
x
,
y
)
∂
x
2
,
∂
2
f
(
x
,
y
)
∂
x
∂
y
,
∂
2
f
(
x
,
y
)
∂
y
2
)
T
.Based on the work of [16], we derive an equivalent formulation of the integrand of (5) in following proposition, whose proof is provided in Appendix A.Proposition 1.
LetW
=
(
1
/
(
2
2
)
)
[
1
0
-
1
2
0
2
0
2
0
]; then theL
2
[
0,2
π
]norm of the second degree directional derivative offat coordinate(
x
,
y
)is equal to∥
W
g
(
x
,
y
)
∥
2, where∥
·
∥
2stands for the Euclidean norm.
Therefore, according to the differential operatorUdefined in Table 1, the ISDTV regularizer in (5) can be equivalently expressed as
(9)
J
2
(
f
)
=
∫
Ω
∥
W
U
f
(
x
,
y
)
∥
2
d
x
d
y
.Table 1
Definition of differential operators.
Vectorial  operators
U
=
(
∂
x
x
,
∂
x
y
,
∂
y
y
)
T
V
=
(
∂
x
,
∂
y
,
∂
x
x
,
∂
x
y
,
∂
y
y
)
T


### 2.3. Mixed Higher Order Total Variation and Equivalent Formulation

Based on the discussion in above Section2.2, TV is reinterpreted as theL
1-L
2mixed norms of first degree directional derivative of the image and can well preserve image edges, and ISDTV is defined as theL
1-L
2mixed norms of second degree directional derivatives of the image and can effectively preserve image features sharpness and minimize the staircase effect. As we know, the recovered image with well-preserved edges and features while not suffering from severely staircase effect is more desirable in image recovery. Furthermore, as detailed in [19–23], the regularizers that combine first and second degree derivatives of the image are reported to well preserve edges and reduce the staircase effect at the same time, and it enables us to extend TV and ISDTV to mixed higher order case. To this end, we thus propose a new regularizer which combines first and second degree derivatives of the image in this subsection.Firstly, with the definitions ofs
0
(
θ
),g
0
(
x
,
y
),s
(
θ
),  andg
(
x
,
y
)in (3), (7), and (8), we define
(10)
S
(
θ
)
=
(
cos
⁡
θ
,
sin
θ
,
cos
⁡
2
θ
,
2
sin
θ
cos
⁡
θ
,
sin
2
θ
)
T
,
G
(
x
,
y
)
=

(

∂
f
(
x
,
y
)
∂
x
,

∂
f
(
x
,
y
)
∂
y
,

∂
2
f
(
x
,
y
)
∂
x
2
,

∂
2
f
(
x
,
y
)
∂
x
∂
y
,

∂
2
f
(
x
,
y
)
∂
y
2

)

T
.
Based on the definitions in (10), we define the new regularizer as
(11)
J
1,2
(
f
)
=
∫
Ω
∥
S
T
(
θ
)
G
(
x
,
y
)
∥
L
2
[
0,2
π
]
d
x
d
y
=
∫
Ω
1
2
π
∫
0
2
π
|
S
T
(
θ
)
G
(
x
,
y
)
|
2
d
θ
d
x
d
y
.Definition (11) leads to a convex regularizer which is also homogeneous, rotation, and translation invariant. Moreover, we have
(12)
S
T
(
θ
)
G
(
x
,
y
)
=
s
0
T
(
θ
)
g
0
(
x
,
y
)
+
s
T
(
θ
)
g
(
x
,
y
)
=
f
θ
,
1
(
x
,
y
)
+
f
θ
,
2
(
x
,
y
)
;
namely, the regularizer in (11) involves the first and second degree directional derivatives of the image which are related to TV and ISDTV, respectively, and takes full advantage of the ability to preserve edges and minimize the staircase effect. Thus, the resulting regularizer constitutes a valid extension of TV and ISDTV. For simplicity, we term the new regularizer in (11) as mixed higher order total variation (MHOTV).Based on the following proposition whose proof is given in AppendixB, we give an equivalent formulation of the MHOTV regularizer.Proposition 2.
LetR
=
(
1
/
(
2
2
)
)
[
0
0
1
0
-
1
2
0
0
0
0
0
2
0
0
0
0
0
2
0
2
0
0
0
2
0
]; then the integrand of (11) is equal to∥
R
G
(
x
,
y
)
∥
2, where∥
·
∥
2stands for the Euclidean norm.
Therefore, according to the differential operatorVdefined in Table 1, we further express our MHOTV regularizer as
(13)
J
1,2
(
f
)
=
∫
Ω
∥
R
V
f
(
x
,
y
)
∥
2
d
x
d
y
.
With the equivalent formulations of the ISDTV and MHOTV regularizers defined in (9) and (13), which can be also understood as weightedL
1-L
2mixed norm of image derivatives, thus the two equivalent formulations are more preferable for the minimization algorithm we will propose in Section 3.


## 2.1. Regularization for Image Recovery Problem

We first consider the recovery of a continuously differentiable imagef
:
Ω
→
ℜfrom its degraded observationg, whereΩ
⊂
ℜ
2 is the spatial support of the image. The most commonly used image observation model can be linearly formulated as
(1)
g
=
A
f
+
e
,
whereAis a linear blur operator andeis Gaussian white noise with standard deviationσ. As we know, the recovery offfrom its observationgis an ill-posed problem, due to the fact that the blur operatorAis ill-conditioned. To stabilize the recovery process, some available a priori knowledge about the desired image is utilized. Thus, a common variational regularization approach is employed to formulate image recovery as an optimization problem:
(2)
f
^
=
argmin
f
{
C
(
f
)
=
1
2
∫
Ω
|
g
-
A
f
|
2
d
x
d
y
+
λ
J
(
f
)
}
,
where the first term is known as data fidelity term, which measures the consistency between estimation and observation, while the second one is called regularization term which provides some a priori knowledge of the desired image, and the regularization parameterλ
≥
0provides a balance between the two terms.


## 2.2. Isotropic Second Degree Total Variation and Equivalent Formulation

Hu and Jacob have introduced the isotropic higher degree total variation (I-HDTV) scheme [16], which inherits the desirable properties of TV, such as, convexity and invariance to rotations and translations. As detailed in [16], the first degree directional derivative offalong the unit vectors
0
(
θ
)
=
(
cos
⁡
θ
,
sin
θ
)
Tat coordinate(
x
,
y
)denoted byf
θ
,
1
(
x
,
y
)is defined as
(3)
f
θ
,
1
(
x
,
y
)
=
[
cos
⁡
θ
,
sin
θ
]
︸
s
0
T
(
θ
)
[
∂
f
(
x
,
y
)
∂
x
∂
f
(
x
,
y
)
∂
y
]
︸
g
0
(
x
,
y
)
,
whereg
0
(
x
,
y
)denotes the gradient offat coordinate(
x
,
y
)and(
·
)
T denotes transpose operator. Specially, as proved in [16], the standard TV is reinterpreted as theL
1-L
2mixed norms of the first degree directional derivative, formulated as
(4)
TV
(
f
)
=
∫
Ω
∥
g
0
(
x
,
y
)
∥
2
d
x
d
y
=
2
∫
Ω
∥
f
θ
,
1
(
x
,
y
)
∥
L
2
[
0,2
π
]
d
x
d
y
=
2
∫
Ω
1
2
π
∫
0
2
π
|
f
θ
,
1
(
x
,
y
)
|
2
d
θ
d
x
d
y
,
where∥
·
∥
2denotes the Euclidean norm.Now, we consider the isotropic second degree total variation (ISDTV) of imagef, which is defined as
(5)
J
2
(
f
)
=
∫
Ω
∥
f
θ
,
2
(
x
,
y
)
∥
L
2
[
0,2
π
]
d
x
d
y
=
∫
Ω
1
2
π
∫
0
2
π
|
f
θ
,
2
(
x
,
y
)
|
2
d
θ
d
x
d
y
,
wheref
θ
,
2
(
x
,
y
)is the second degree directional derivative of imagefat coordinate(
x
,
y
), defined as
(6)
f
θ
,
2
(
x
,
y
)
=
s
T
(
θ
)
g
(
x
,
y
)
,
with
(7)
s
(
θ
)
=
(
cos
⁡
2
θ
,
2
sin
θ
cos
⁡
θ
,
si
n
2
θ
)
T
,
(8)
g
(
x
,
y
)
=
(
∂
2
f
(
x
,
y
)
∂
x
2
,
∂
2
f
(
x
,
y
)
∂
x
∂
y
,
∂
2
f
(
x
,
y
)
∂
y
2
)
T
.Based on the work of [16], we derive an equivalent formulation of the integrand of (5) in following proposition, whose proof is provided in Appendix A.Proposition 1.
LetW
=
(
1
/
(
2
2
)
)
[
1
0
-
1
2
0
2
0
2
0
]; then theL
2
[
0,2
π
]norm of the second degree directional derivative offat coordinate(
x
,
y
)is equal to∥
W
g
(
x
,
y
)
∥
2, where∥
·
∥
2stands for the Euclidean norm.
Therefore, according to the differential operatorUdefined in Table 1, the ISDTV regularizer in (5) can be equivalently expressed as
(9)
J
2
(
f
)
=
∫
Ω
∥
W
U
f
(
x
,
y
)
∥
2
d
x
d
y
.Table 1
Definition of differential operators.
Vectorial  operators
U
=
(
∂
x
x
,
∂
x
y
,
∂
y
y
)
T
V
=
(
∂
x
,
∂
y
,
∂
x
x
,
∂
x
y
,
∂
y
y
)
T


## 2.3. Mixed Higher Order Total Variation and Equivalent Formulation

Based on the discussion in above Section2.2, TV is reinterpreted as theL
1-L
2mixed norms of first degree directional derivative of the image and can well preserve image edges, and ISDTV is defined as theL
1-L
2mixed norms of second degree directional derivatives of the image and can effectively preserve image features sharpness and minimize the staircase effect. As we know, the recovered image with well-preserved edges and features while not suffering from severely staircase effect is more desirable in image recovery. Furthermore, as detailed in [19–23], the regularizers that combine first and second degree derivatives of the image are reported to well preserve edges and reduce the staircase effect at the same time, and it enables us to extend TV and ISDTV to mixed higher order case. To this end, we thus propose a new regularizer which combines first and second degree derivatives of the image in this subsection.Firstly, with the definitions ofs
0
(
θ
),g
0
(
x
,
y
),s
(
θ
),  andg
(
x
,
y
)in (3), (7), and (8), we define
(10)
S
(
θ
)
=
(
cos
⁡
θ
,
sin
θ
,
cos
⁡
2
θ
,
2
sin
θ
cos
⁡
θ
,
sin
2
θ
)
T
,
G
(
x
,
y
)
=

(

∂
f
(
x
,
y
)
∂
x
,

∂
f
(
x
,
y
)
∂
y
,

∂
2
f
(
x
,
y
)
∂
x
2
,

∂
2
f
(
x
,
y
)
∂
x
∂
y
,

∂
2
f
(
x
,
y
)
∂
y
2

)

T
.
Based on the definitions in (10), we define the new regularizer as
(11)
J
1,2
(
f
)
=
∫
Ω
∥
S
T
(
θ
)
G
(
x
,
y
)
∥
L
2
[
0,2
π
]
d
x
d
y
=
∫
Ω
1
2
π
∫
0
2
π
|
S
T
(
θ
)
G
(
x
,
y
)
|
2
d
θ
d
x
d
y
.Definition (11) leads to a convex regularizer which is also homogeneous, rotation, and translation invariant. Moreover, we have
(12)
S
T
(
θ
)
G
(
x
,
y
)
=
s
0
T
(
θ
)
g
0
(
x
,
y
)
+
s
T
(
θ
)
g
(
x
,
y
)
=
f
θ
,
1
(
x
,
y
)
+
f
θ
,
2
(
x
,
y
)
;
namely, the regularizer in (11) involves the first and second degree directional derivatives of the image which are related to TV and ISDTV, respectively, and takes full advantage of the ability to preserve edges and minimize the staircase effect. Thus, the resulting regularizer constitutes a valid extension of TV and ISDTV. For simplicity, we term the new regularizer in (11) as mixed higher order total variation (MHOTV).Based on the following proposition whose proof is given in AppendixB, we give an equivalent formulation of the MHOTV regularizer.Proposition 2.
LetR
=
(
1
/
(
2
2
)
)
[
0
0
1
0
-
1
2
0
0
0
0
0
2
0
0
0
0
0
2
0
2
0
0
0
2
0
]; then the integrand of (11) is equal to∥
R
G
(
x
,
y
)
∥
2, where∥
·
∥
2stands for the Euclidean norm.
Therefore, according to the differential operatorVdefined in Table 1, we further express our MHOTV regularizer as
(13)
J
1,2
(
f
)
=
∫
Ω
∥
R
V
f
(
x
,
y
)
∥
2
d
x
d
y
.
With the equivalent formulations of the ISDTV and MHOTV regularizers defined in (9) and (13), which can be also understood as weightedL
1-L
2mixed norm of image derivatives, thus the two equivalent formulations are more preferable for the minimization algorithm we will propose in Section 3.


## 3. FPG-MFISTA: The Efficient Algorithm for the Proposed Model

Using the ISDTV and MHOTV regularizers described in (9) and (13), we present a unified variational regularization image recovery model to reconstruct the desired image from its noisy observation, which is formulated as
(14)
f
^
=
argmin
f
{
C
(
f
)
=
1
2
∫
Ω
|
g
-
A
f
|
2
d
x
d
y
+
λ
J
(
f
)
}
J
(
f
)
∈
{
J
2
(
f
)
,
J
1,2
(
f
)
}
,
where the last expression in (14) means that the regularization term can be either ISDTV or MHOTV regularizer.


### 3.1. Discrete Formulation for ISDTV- and MHOTV-Based Image Recovery Model

Now we consider the discrete formulation of proposed model (14). In order to simplify our analysis, we assume that the image with sizem
×
nis stacked in a vector of sizeN
=
m
×
n. Then, the objective function in (14) can be discretized as
(15)
C
(
f
)
=
1
2
∥
g
-
A
f
∥
2
2
+
λ
J
(
f
)
,
whereA
∈
ℜ
N
×
Nis the blurring matrix andg,f
∈
ℜ
Nare the vectorized versions of the observed image and the image to be recovered, respectively. To discretize the differential operators defined in Table 1, we assume reflexive boundary conditions for the images and use forward finite differences to approximate the partial derivatives [27].To simplify our analysis, we will introduce some auxiliary definitions and notations as follows.First, we define operatorsU
:
ℜ
N
→
ℜ
3
×
N andV
:
ℜ
N
→
ℜ
5
×
Nto refer to the discrete versions of the operatorsUandVdefined in Table 1, respectively. For simplicity, we denote spaceX
j
=
ℜ
j
×
N
=
{
u
=
(
u
1
,
u
2
,
…
,
u
N
)
∣
u
i
∈
ℜ
j
,
∀
i
=
1,2
,
…
,
N
} and equip the spaceX
jwith the inner product〈
·
,
·
〉
X
jand norm∥
·
∥
X
j. Letu
=
(
u
1
,
u
2
,
…
,
u
N
)
∈
X
jandv
=
(
v
1
,
v
2
,
…
,
v
N
)
∈
X
j; we define them as
(16)
〈
u
,
v
〉
X
j
=
∑
i
=
1
N
u
i
T
v
i
,
∥
u
∥
X
j
=
〈
u
,
u
〉
X
j
.Next, let us define the adjoint operators ofUandVwhich are the discrete operatorsU
T
:
X
3
→
ℜ
N andV
T
:
X
5
→
ℜ
N, respectively. Thus, letp
∈
X
3, q
∈
X
5andf
∈
ℜ
N; we have
(17)
〈
p
,
U
f
〉
X
3
=
〈
U
T
p
,
f
〉
2
,
〈
q
,
V
f
〉
X
5
=
〈
V
T
q
,
f
〉
2
,
where〈
·
,
·
〉
2stands for the inner product in Euclidean spaceℜ
N.Furthermore, with thel
1
-
l
2mixed norm of a vector fieldu
=
(
u
1
,
u
2
,
…
,
u
N
)
∈
X
jdefined as∥
u
∥
1
,
2
=
∑
i
=
1
N
∥
u
i
∥
2, we thus write the discrete versions of the ISDTV and MHOTV regularizers defined in (9) and (13) as
(18)
J
2
(
f
)
=
∑
i
=
1
N
∥
W
(
U
f
)
i
∥
2
=
∥
W
U
f
∥
1
,
2
,
(19)
J
1,2
(
f
)
=
∑
i
=
1
N
∥
R
(
V
f
)
i
∥
2
=
∥
R
V
f
∥
1
,
2
,
where(
·
)
idenotes theith column element of the corresponding vector filed.Therefore, we can write the discrete formulation of proposed model (14) as
(20)
f
^
=
arg
min
⁡
f
{
C
(
f
)
=
1
2
∥
g
-
A
f
∥
2
2
+
λ
∥
H
f
∥
1
,
2
}
H
∈
{
W
U
,
R
V
}
,
with the last expressionH
∈
{
W
U
,
R
V
}in (20) meaning that the regularization term can be either ISDTV or MHOTV regularizer.


### 3.2. Majorization-Minimization (MM)

In the following, we focus on the minimization of model (20). Due to the fact that the ISDTV and MHOTV regularizers are both nondifferentiable, we solve model (20) by a majorization-minimization (MM) approach [28–30]. As described in the MM framework, an iterative algorithm for solving the proposed model (20) atkth iteratef
(
k
)involves the following two steps:(1)
constructing the surrogate functionalL
(
f
;
f
(
k
)
)that satisfies the following properties:
(21)
C
(
f
)
≤
L
(
f
;
f
(
k
)
)
,
∀
f
;
C
(
f
(
k
)
)
=
L
(
f
(
k
)
;
f
(
k
)
)
,
(2)
solving the next iteratef
(
k
+
1
)
=
argmin
f
⁡
L
(
f
;
f
(
k
)
).Under the MM framework, we find the optimal solution of model (20) by iteratively minimizing a sequence of surrogate functionals{
L
(
f
;
f
(
k
)
)
}that upper bound the objective functionC
(
f
).Therefore, to obtain the surrogate functionals{
L
(
f
;
f
(
k
)
)
}, we first upper bound the data fidelity term ofC
(
f
)atkth iteratef
(
k
)using the following majorizer:
(22)
Q
(
f
;
f
(
k
)
)
=
α
2
∥
f
-
z
∥
2
2
+
c
,
wherez
=
f
(
k
)
+
(
1
/
α
)
A
T
(
g
-
A
f
(
k
)
)andcis a constant which is independent off. To make sure thatQ
(
f
;
f
(
k
)
)is valid, we chooseα
>
∥
A
T
A
∥to makeα
I
-
A
T
Apositive definite. Thus, the whole resulting majorizer ofC
(
f
)atkth iteratef
(
k
)can be written as
(23)
L
(
f
;
f
(
k
)
)
=
Q
(
f
;
f
(
k
)
)
+
λ
∥
H
f
∥
1
,
2
=
α
2
∥
f
-
z
∥
2
2
+
λ
∥
H
f
∥
1
,
2
+
c
.Then, the minimization problem (20) is converted to iteratively minimizing (23) with regard tof, which can be formulated as
(24)
f
(
k
+
1
)
=
arg
min
⁡
f
{
L
(
f
;
f
(
k
)
)
=
α
2
∥
f
-
z
∥
2
2
+
λ
∥
H
f
∥
1
,
2
+
c
}
H
∈
{
W
U
,
R
V
}
.As we can see, the resulting minimization problem (24) is much simpler and its solution can be alternatively interpreted as the solution of a denoising problem withzbeing the noisy observation. In the following, we will present a primal-dual approach based on the monotone fast iterative shrinkage thresholding algorithm (MFISTA) [24] to seek the optimal solution of (24).


### 3.3. FPG-MFISTA

Considering thatfis a digital image whose intensities are finite and bounded, we can use a convex closed set to constrain the range of image intensities. Based on that, the resulting denoising subproblem (24) can be equivalently formulated as
(25)
f
(
k
+
1
)
=
arg
min
⁡
f
α
2
∥
f
-
z
∥
2
2
+
λ
∥
H
f
∥
1
,
2
+
I
S
(
f
)
=
arg
min
⁡
f
∈
S
α
2
∥
f
-
z
∥
2
2
+
λ
∥
H
f
∥
1
,
2
,
whereH
∈
{
W
U
,
R
V
}means that the regularization term can be either ISDTV or MHOTV regularizer,S
=
{
f
∈
ℜ
N
∣
f
i
∈
[
l
,
u
]
,
∀
i
=
1
,
…
,
N
}is a convex closed set which enforces bounded constraints on the solution, andI
S
(
f
)is the indicator function ofSthat takes the value0forf
∈
Sand∞otherwise.Next, we present a primal-dual approach which can efficiently compute the solution of the constrained problem (25). Without loss of generality, here we just describe all the steps of the minimization algorithm for model (25) involving the proposed MHOTV regularizer, which corresponds to the case thatH
=
R
V.Therefore, by denotingμ
=
λ
/
α, the corresponding MHOTV based denoising subproblem can be further rewritten as
(26)
f
(
k
+
1
)
=
arg
min
⁡
f
∈
S
1
2
∥
f
-
z
∥
2
2
+
μ
∥
R
V
f
∥
1
,
2
.As referred in [31],∥
·
∥
∞
,
2is the dual norm of∥
·
∥
1
,
2, and then thel
1
-
l
2mixed norm of the vector fieldu
=
(
u
1
,
u
2
,
…
,
u
N
)
∈
X
jcan be rewritten as
(27)
∥
u
∥
1
,
2
=
max
⁡
ω
∈
X
〈
u
,
ω
〉
X
j
,
whereX
=
{
ω
=
(
ω
1
,
ω
2
,
…
,
ω
N
)
∈
X
j
∣
∥
ω
i
∥
2
≤
1
,
∀
i
=
1
,
…
,
N
}denotes thel
∞
-
l
2unit norm ball. Using (27), we rewrite model (26) as
(28)
f
(
k
+
1
)
=
arg
min
⁡
f
∈
S
1
2
∥
f
-
z
∥
2
2
+
μ
max
⁡
ω
∈
X
〈
R
V
f
,
ω
〉
X
5
=
arg
min
⁡
f
∈
S
1
2
∥
f
-
z
∥
2
2
+
μ
max
⁡
ω
∈
X
〈
f
,
V
T
R
T
ω
〉
2
;
thus, the resulting formulation (28) makes us to solve the following minimax problem:
(29)
min
⁡
f
∈
S
max
⁡
ω
∈
X
{
E
(
f
,
ω
)
=
1
2
∥
f
-
z
∥
2
2
+
μ
〈
f
,
V
T
R
T
ω
〉
2
}
.Since the cost functionE
(
f
,
ω
)is strictly convex infand concave inω, we can exchange the order of the minimum and maximum, and the saddle-point(
f
(
k
+
1
)
,
ω
(
k
)
)of (29) is attained as reported in [32]. Therefore,
(30)
min
⁡
f
∈
S
max
⁡
ω
∈
X
E
(
f
,
ω
)
=
E
(
f
(
k
+
1
)
,
ω
(
k
)
)
=
max
⁡
ω
∈
X
min
⁡
f
∈
S
E
(
f
,
ω
)
.As we know,E
(
f
,
ω
)can be equivalently written as
(31)
E
(
f
,
ω
)
=
1
2
∥
f
-
(
z
-
μ
V
T
R
T
ω
)
∥
2
2
+
1
2
∥
z
∥
2
2
-
1
2
∥
z
-
μ
V
T
R
T
ω
∥
2
2
.
Based on (30) and (31), we can obtain the saddle-point(
f
(
k
+
1
)
,
ω
(
k
)
)of (29) by solving
(32)
f
(
k
+
1
)
=
argmin
f
∈
S
{
max
⁡
ω
∈
X
1
2
∥
f
-
(
z
-
μ
V
T
R
T
ω
)
∥
2
2
+
1
2
∥
z
∥
2
2
-
1
2
∥
z
-
μ
V
T
R
T
ω
∥
2
2
}
=
P
S
(
z
-
μ
V
T
R
T
ω
(
k
)
)
,
whereP
Sdenotes the orthogonal projection onto the convex setS, whileω
(
k
)is the maximizer of the dual problem, formulated as
(33)
ω
(
k
)
=
argmax
ω
∈
X
{
h
(
ω
)
=
1
2
∥
P
S
(
z
-
μ
V
T
R
T
ω
)
-
(
z
-
μ
V
T
R
T
ω
)
∥
2
2
+
1
2
∥
z
∥
2
2
-
1
2
∥
z
-
μ
V
T
R
T
ω
∥
2
2
}
.Since the dual functionh
(
ω
)in (33) is differentiable and has Lipschitz continuous gradient, as referred in [24], then we can compute the gradient ofh
(
ω
)as
(34)
∇
h
(
ω
)
=
μ
R
V
P
S
(
z
-
μ
V
T
R
T
ω
)
.To accelerate the convergence, we use Nesterov’s iterative method [33] to solve model (33), which is a gradient-based approach that exhibits convergence rates of one order higher than the standard gradient ascent method. Furthermore, we employ the operatorP
Xto return the orthogonal projection onto thel
∞
-
l
2unit norm ballX, where this operation is performed by projecting independently each componentω
i,  i
=
1,2
,
…
,
N ofωonto thel
2unit norm ball.Therefore, as described above, we first find the maximizer of the dual problem (33) and then obtain the optimal solution of our primal problem (26) through (32).In the same way, by using (28)–(33) and replacingRand V withWandU, respectively, we can obtain the solution of model (25) involving the ISDTV regularizer.Finally, a detailed description of the fast projected gradient algorithm combined with monotone fast iterative shrinkage thresholding (FPG-MFISTA) for ISDTV and MHOTV regularization-based image recovery is given in Algorithm1, while the corresponding denoising subproblem (24) can be efficiently solved by a fast projected gradient algorithm described in Algorithm 2.Algorithm 1:FPG-MFISTA for ISDTV and MHOTV regularizer-based image recovery.
Step  1. Input g, A, λ
>
0, α
>
∥
A
T
A
∥, ε (ε
>
0 is a tolerance), M
a
x
i
t
e
r
(Number of iterations). Set  k
=
0. Initialize v
(
1
)
=
f
(
0
)
,
t
1
=
1
,
c
1
=
C
(
f
(
0
)
).
Step  2. Compute
s
(
k
)
⟵
FPG
(
v
(
k
)
+
1
α
A
T
(
g
-
A
v
(
k
)
)
,
λ
α
);
t
k
+
1
⟵
1
+
1
+
4
t
k
2
2
;
c
k
+
1
=
C
(
s
(
k
)
);
if  c
k
+
1
>
c
k  then
c
k
+
1
=
c
k
;
f
(
k
+
1
)
⟵
f
(
k
);
else
f
(
k
+
1
)
⟵
s
(
k
);
end
v
(
k
+
1
)
⟵
f
(
k
+
1
)
+
t
k
t
k
+
1
(
s
(
k
)
-
f
(
k
+
1
)
)
+
t
k
-
1
t
k
+
1
(
f
(
k
+
1
)
-
f
(
k
)
).
Step  3. Terminate if  ∥
f
(
k
+
1
)
-
f
(
k
)
∥
2
/
∥
f
(
k
)
∥
2
<
ε. Otherwise, set k
=
k
+
1 and go to Step  2
Step  4. Output f
^-An optimal solution of (20).Algorithm 2:FPG (
z
,
μ
) (fast projected gradient) algorithm for image denoising under ISDTV (Η
=
W
U) or MHOTV (Η
=
R
V) regularizer.
Step  1. Input z
,
Η
,
μ
>
0
,
γ
>
∥
Η
T
Η
∥, ε
0
(
ε
0
>
0 is a tolerance), I
n
n
e
r
i
t
e
r
(Number of iterations). Setk
=
0. Initialize u
(
1
)
=
ω
(
0
)
=
0, t
1
=
1.
Step  2. Compute
ω
(
k
)
⟵
P
X
(
u
(
k
)
+
1
μ
γ
Η
P
S
(
z
-
μ
Η
T
u
(
k
)
)
);
t
k
+
1
⟵
1
+
1
+
4
t
k
2
2;
u
(
k
+
1
)
⟵
ω
(
k
)
+
t
k
-
1
t
k
+
1
(
ω
(
k
)
-
ω
(
k
-
1
)
);
Step  3. Terminate if ∥
ω
(
k
+
1
)
-
ω
(
k
)
∥
2
/
∥
ω
(
k
)
∥
2
<
ε
0. Otherwise, set k
=
k
+
1 and go to Step  2
Step  4. Output f
(
k
)
=
P
S
(
z
-
μ
Η
T
ω
(
k
-
1
)
)—An optimal solution of (25).


## 3.1. Discrete Formulation for ISDTV- and MHOTV-Based Image Recovery Model

Now we consider the discrete formulation of proposed model (14). In order to simplify our analysis, we assume that the image with sizem
×
nis stacked in a vector of sizeN
=
m
×
n. Then, the objective function in (14) can be discretized as
(15)
C
(
f
)
=
1
2
∥
g
-
A
f
∥
2
2
+
λ
J
(
f
)
,
whereA
∈
ℜ
N
×
Nis the blurring matrix andg,f
∈
ℜ
Nare the vectorized versions of the observed image and the image to be recovered, respectively. To discretize the differential operators defined in Table 1, we assume reflexive boundary conditions for the images and use forward finite differences to approximate the partial derivatives [27].To simplify our analysis, we will introduce some auxiliary definitions and notations as follows.First, we define operatorsU
:
ℜ
N
→
ℜ
3
×
N andV
:
ℜ
N
→
ℜ
5
×
Nto refer to the discrete versions of the operatorsUandVdefined in Table 1, respectively. For simplicity, we denote spaceX
j
=
ℜ
j
×
N
=
{
u
=
(
u
1
,
u
2
,
…
,
u
N
)
∣
u
i
∈
ℜ
j
,
∀
i
=
1,2
,
…
,
N
} and equip the spaceX
jwith the inner product〈
·
,
·
〉
X
jand norm∥
·
∥
X
j. Letu
=
(
u
1
,
u
2
,
…
,
u
N
)
∈
X
jandv
=
(
v
1
,
v
2
,
…
,
v
N
)
∈
X
j; we define them as
(16)
〈
u
,
v
〉
X
j
=
∑
i
=
1
N
u
i
T
v
i
,
∥
u
∥
X
j
=
〈
u
,
u
〉
X
j
.Next, let us define the adjoint operators ofUandVwhich are the discrete operatorsU
T
:
X
3
→
ℜ
N andV
T
:
X
5
→
ℜ
N, respectively. Thus, letp
∈
X
3, q
∈
X
5andf
∈
ℜ
N; we have
(17)
〈
p
,
U
f
〉
X
3
=
〈
U
T
p
,
f
〉
2
,
〈
q
,
V
f
〉
X
5
=
〈
V
T
q
,
f
〉
2
,
where〈
·
,
·
〉
2stands for the inner product in Euclidean spaceℜ
N.Furthermore, with thel
1
-
l
2mixed norm of a vector fieldu
=
(
u
1
,
u
2
,
…
,
u
N
)
∈
X
jdefined as∥
u
∥
1
,
2
=
∑
i
=
1
N
∥
u
i
∥
2, we thus write the discrete versions of the ISDTV and MHOTV regularizers defined in (9) and (13) as
(18)
J
2
(
f
)
=
∑
i
=
1
N
∥
W
(
U
f
)
i
∥
2
=
∥
W
U
f
∥
1
,
2
,
(19)
J
1,2
(
f
)
=
∑
i
=
1
N
∥
R
(
V
f
)
i
∥
2
=
∥
R
V
f
∥
1
,
2
,
where(
·
)
idenotes theith column element of the corresponding vector filed.Therefore, we can write the discrete formulation of proposed model (14) as
(20)
f
^
=
arg
min
⁡
f
{
C
(
f
)
=
1
2
∥
g
-
A
f
∥
2
2
+
λ
∥
H
f
∥
1
,
2
}
H
∈
{
W
U
,
R
V
}
,
with the last expressionH
∈
{
W
U
,
R
V
}in (20) meaning that the regularization term can be either ISDTV or MHOTV regularizer.


## 3.2. Majorization-Minimization (MM)

In the following, we focus on the minimization of model (20). Due to the fact that the ISDTV and MHOTV regularizers are both nondifferentiable, we solve model (20) by a majorization-minimization (MM) approach [28–30]. As described in the MM framework, an iterative algorithm for solving the proposed model (20) atkth iteratef
(
k
)involves the following two steps:(1)
constructing the surrogate functionalL
(
f
;
f
(
k
)
)that satisfies the following properties:
(21)
C
(
f
)
≤
L
(
f
;
f
(
k
)
)
,
∀
f
;
C
(
f
(
k
)
)
=
L
(
f
(
k
)
;
f
(
k
)
)
,
(2)
solving the next iteratef
(
k
+
1
)
=
argmin
f
⁡
L
(
f
;
f
(
k
)
).Under the MM framework, we find the optimal solution of model (20) by iteratively minimizing a sequence of surrogate functionals{
L
(
f
;
f
(
k
)
)
}that upper bound the objective functionC
(
f
).Therefore, to obtain the surrogate functionals{
L
(
f
;
f
(
k
)
)
}, we first upper bound the data fidelity term ofC
(
f
)atkth iteratef
(
k
)using the following majorizer:
(22)
Q
(
f
;
f
(
k
)
)
=
α
2
∥
f
-
z
∥
2
2
+
c
,
wherez
=
f
(
k
)
+
(
1
/
α
)
A
T
(
g
-
A
f
(
k
)
)andcis a constant which is independent off. To make sure thatQ
(
f
;
f
(
k
)
)is valid, we chooseα
>
∥
A
T
A
∥to makeα
I
-
A
T
Apositive definite. Thus, the whole resulting majorizer ofC
(
f
)atkth iteratef
(
k
)can be written as
(23)
L
(
f
;
f
(
k
)
)
=
Q
(
f
;
f
(
k
)
)
+
λ
∥
H
f
∥
1
,
2
=
α
2
∥
f
-
z
∥
2
2
+
λ
∥
H
f
∥
1
,
2
+
c
.Then, the minimization problem (20) is converted to iteratively minimizing (23) with regard tof, which can be formulated as
(24)
f
(
k
+
1
)
=
arg
min
⁡
f
{
L
(
f
;
f
(
k
)
)
=
α
2
∥
f
-
z
∥
2
2
+
λ
∥
H
f
∥
1
,
2
+
c
}
H
∈
{
W
U
,
R
V
}
.As we can see, the resulting minimization problem (24) is much simpler and its solution can be alternatively interpreted as the solution of a denoising problem withzbeing the noisy observation. In the following, we will present a primal-dual approach based on the monotone fast iterative shrinkage thresholding algorithm (MFISTA) [24] to seek the optimal solution of (24).


## 3.3. FPG-MFISTA

Considering thatfis a digital image whose intensities are finite and bounded, we can use a convex closed set to constrain the range of image intensities. Based on that, the resulting denoising subproblem (24) can be equivalently formulated as
(25)
f
(
k
+
1
)
=
arg
min
⁡
f
α
2
∥
f
-
z
∥
2
2
+
λ
∥
H
f
∥
1
,
2
+
I
S
(
f
)
=
arg
min
⁡
f
∈
S
α
2
∥
f
-
z
∥
2
2
+
λ
∥
H
f
∥
1
,
2
,
whereH
∈
{
W
U
,
R
V
}means that the regularization term can be either ISDTV or MHOTV regularizer,S
=
{
f
∈
ℜ
N
∣
f
i
∈
[
l
,
u
]
,
∀
i
=
1
,
…
,
N
}is a convex closed set which enforces bounded constraints on the solution, andI
S
(
f
)is the indicator function ofSthat takes the value0forf
∈
Sand∞otherwise.Next, we present a primal-dual approach which can efficiently compute the solution of the constrained problem (25). Without loss of generality, here we just describe all the steps of the minimization algorithm for model (25) involving the proposed MHOTV regularizer, which corresponds to the case thatH
=
R
V.Therefore, by denotingμ
=
λ
/
α, the corresponding MHOTV based denoising subproblem can be further rewritten as
(26)
f
(
k
+
1
)
=
arg
min
⁡
f
∈
S
1
2
∥
f
-
z
∥
2
2
+
μ
∥
R
V
f
∥
1
,
2
.As referred in [31],∥
·
∥
∞
,
2is the dual norm of∥
·
∥
1
,
2, and then thel
1
-
l
2mixed norm of the vector fieldu
=
(
u
1
,
u
2
,
…
,
u
N
)
∈
X
jcan be rewritten as
(27)
∥
u
∥
1
,
2
=
max
⁡
ω
∈
X
〈
u
,
ω
〉
X
j
,
whereX
=
{
ω
=
(
ω
1
,
ω
2
,
…
,
ω
N
)
∈
X
j
∣
∥
ω
i
∥
2
≤
1
,
∀
i
=
1
,
…
,
N
}denotes thel
∞
-
l
2unit norm ball. Using (27), we rewrite model (26) as
(28)
f
(
k
+
1
)
=
arg
min
⁡
f
∈
S
1
2
∥
f
-
z
∥
2
2
+
μ
max
⁡
ω
∈
X
〈
R
V
f
,
ω
〉
X
5
=
arg
min
⁡
f
∈
S
1
2
∥
f
-
z
∥
2
2
+
μ
max
⁡
ω
∈
X
〈
f
,
V
T
R
T
ω
〉
2
;
thus, the resulting formulation (28) makes us to solve the following minimax problem:
(29)
min
⁡
f
∈
S
max
⁡
ω
∈
X
{
E
(
f
,
ω
)
=
1
2
∥
f
-
z
∥
2
2
+
μ
〈
f
,
V
T
R
T
ω
〉
2
}
.Since the cost functionE
(
f
,
ω
)is strictly convex infand concave inω, we can exchange the order of the minimum and maximum, and the saddle-point(
f
(
k
+
1
)
,
ω
(
k
)
)of (29) is attained as reported in [32]. Therefore,
(30)
min
⁡
f
∈
S
max
⁡
ω
∈
X
E
(
f
,
ω
)
=
E
(
f
(
k
+
1
)
,
ω
(
k
)
)
=
max
⁡
ω
∈
X
min
⁡
f
∈
S
E
(
f
,
ω
)
.As we know,E
(
f
,
ω
)can be equivalently written as
(31)
E
(
f
,
ω
)
=
1
2
∥
f
-
(
z
-
μ
V
T
R
T
ω
)
∥
2
2
+
1
2
∥
z
∥
2
2
-
1
2
∥
z
-
μ
V
T
R
T
ω
∥
2
2
.
Based on (30) and (31), we can obtain the saddle-point(
f
(
k
+
1
)
,
ω
(
k
)
)of (29) by solving
(32)
f
(
k
+
1
)
=
argmin
f
∈
S
{
max
⁡
ω
∈
X
1
2
∥
f
-
(
z
-
μ
V
T
R
T
ω
)
∥
2
2
+
1
2
∥
z
∥
2
2
-
1
2
∥
z
-
μ
V
T
R
T
ω
∥
2
2
}
=
P
S
(
z
-
μ
V
T
R
T
ω
(
k
)
)
,
whereP
Sdenotes the orthogonal projection onto the convex setS, whileω
(
k
)is the maximizer of the dual problem, formulated as
(33)
ω
(
k
)
=
argmax
ω
∈
X
{
h
(
ω
)
=
1
2
∥
P
S
(
z
-
μ
V
T
R
T
ω
)
-
(
z
-
μ
V
T
R
T
ω
)
∥
2
2
+
1
2
∥
z
∥
2
2
-
1
2
∥
z
-
μ
V
T
R
T
ω
∥
2
2
}
.Since the dual functionh
(
ω
)in (33) is differentiable and has Lipschitz continuous gradient, as referred in [24], then we can compute the gradient ofh
(
ω
)as
(34)
∇
h
(
ω
)
=
μ
R
V
P
S
(
z
-
μ
V
T
R
T
ω
)
.To accelerate the convergence, we use Nesterov’s iterative method [33] to solve model (33), which is a gradient-based approach that exhibits convergence rates of one order higher than the standard gradient ascent method. Furthermore, we employ the operatorP
Xto return the orthogonal projection onto thel
∞
-
l
2unit norm ballX, where this operation is performed by projecting independently each componentω
i,  i
=
1,2
,
…
,
N ofωonto thel
2unit norm ball.Therefore, as described above, we first find the maximizer of the dual problem (33) and then obtain the optimal solution of our primal problem (26) through (32).In the same way, by using (28)–(33) and replacingRand V withWandU, respectively, we can obtain the solution of model (25) involving the ISDTV regularizer.Finally, a detailed description of the fast projected gradient algorithm combined with monotone fast iterative shrinkage thresholding (FPG-MFISTA) for ISDTV and MHOTV regularization-based image recovery is given in Algorithm1, while the corresponding denoising subproblem (24) can be efficiently solved by a fast projected gradient algorithm described in Algorithm 2.Algorithm 1:FPG-MFISTA for ISDTV and MHOTV regularizer-based image recovery.
Step  1. Input g, A, λ
>
0, α
>
∥
A
T
A
∥, ε (ε
>
0 is a tolerance), M
a
x
i
t
e
r
(Number of iterations). Set  k
=
0. Initialize v
(
1
)
=
f
(
0
)
,
t
1
=
1
,
c
1
=
C
(
f
(
0
)
).
Step  2. Compute
s
(
k
)
⟵
FPG
(
v
(
k
)
+
1
α
A
T
(
g
-
A
v
(
k
)
)
,
λ
α
);
t
k
+
1
⟵
1
+
1
+
4
t
k
2
2
;
c
k
+
1
=
C
(
s
(
k
)
);
if  c
k
+
1
>
c
k  then
c
k
+
1
=
c
k
;
f
(
k
+
1
)
⟵
f
(
k
);
else
f
(
k
+
1
)
⟵
s
(
k
);
end
v
(
k
+
1
)
⟵
f
(
k
+
1
)
+
t
k
t
k
+
1
(
s
(
k
)
-
f
(
k
+
1
)
)
+
t
k
-
1
t
k
+
1
(
f
(
k
+
1
)
-
f
(
k
)
).
Step  3. Terminate if  ∥
f
(
k
+
1
)
-
f
(
k
)
∥
2
/
∥
f
(
k
)
∥
2
<
ε. Otherwise, set k
=
k
+
1 and go to Step  2
Step  4. Output f
^-An optimal solution of (20).Algorithm 2:FPG (
z
,
μ
) (fast projected gradient) algorithm for image denoising under ISDTV (Η
=
W
U) or MHOTV (Η
=
R
V) regularizer.
Step  1. Input z
,
Η
,
μ
>
0
,
γ
>
∥
Η
T
Η
∥, ε
0
(
ε
0
>
0 is a tolerance), I
n
n
e
r
i
t
e
r
(Number of iterations). Setk
=
0. Initialize u
(
1
)
=
ω
(
0
)
=
0, t
1
=
1.
Step  2. Compute
ω
(
k
)
⟵
P
X
(
u
(
k
)
+
1
μ
γ
Η
P
S
(
z
-
μ
Η
T
u
(
k
)
)
);
t
k
+
1
⟵
1
+
1
+
4
t
k
2
2;
u
(
k
+
1
)
⟵
ω
(
k
)
+
t
k
-
1
t
k
+
1
(
ω
(
k
)
-
ω
(
k
-
1
)
);
Step  3. Terminate if ∥
ω
(
k
+
1
)
-
ω
(
k
)
∥
2
/
∥
ω
(
k
)
∥
2
<
ε
0. Otherwise, set k
=
k
+
1 and go to Step  2
Step  4. Output f
(
k
)
=
P
S
(
z
-
μ
Η
T
ω
(
k
-
1
)
)—An optimal solution of (25).


## 4. Experimental Results

To evaluate the performance of our proposed ISDTV and MHOTV regularization schemes in the context of image deblurring, we provide systematical experimental comparisons with five state-of-the-art methods, namely, TV regularization method [24], nonlocal TV regularization method [25], Hessian Frobenius-norm and Hessian Spectral-norm regularization methods [18], and the Hessian Nuclear-norm regularization method [26]. Finally, the quality of the reconstructed image is evaluated in terms of peak signal-to-noise ratio (PSNR), which is measured in decibels (dB) and defined as
(35)
PSNR
=
10
lo
g
10
(
max
⁡
1
≤
i
≤
N
⁡
|
f
i
(
0
)
|
2
(
1
/
N
)
∥
f
(
0
)
-
f
^
∥
2
2
)
,
wheref
^
∈
ℜ
Nandf
(
0
)
∈
ℜ
Nare the vectorized versions of the reconstructed image and the original image, respectively.


### 4.1. Experiment Setting

For the image deblurring experiments, our test images contain six natural images and six biomedical cell images which are shown in Figures1 and 2, respectively, where the cell images are part of the biomedical image database [34], and they are all converted to grayscale. In our experiments, the intensities of all the test images are initially normalized to the range[
0,1
], and we use three kinds of point spread function (PSF) to produce the blurred images, which are a9
×
9Gaussian PSF of standard deviationσ
b
=
6, a9
×
9uniform one, and a19
×
19motion one, respectively. Furthermore, the Gaussian noise of two different levels corresponding to a blurred SNR (BSNR) of{
20,30
} dB is added to generate the eventually degraded images, and the BSNR is defined as
(36)
BSNR
=
var
(
A
f
)
σ
n
2
,
wherevar
(
A
f
)is the variance of blurred image andσ
nis the standard deviation of Gaussian noise. In addition, we constrain the intensities of image recovered by all the seven methods to lie in the convex closed setS
=
{
f
∈
ℜ
N
∣
f
i
∈
[
0
,
1
]
,
∀
i
=
1
,
…
,
N
}. For the sake of fair comparisons, the regularization parameterλfor each regularization method is chosen to give the best PSNR performance. Finally, we setε
=
1
0
-
5,M
a
x
i
t
e
r
=
100,ε
0
=
1
0
-
3
,andI
n
n
e
r
i
t
e
r
=
10in our experiments.Figure 1
Set of natural images (left to right): Barbara, Boat, Face, Fingerprint, Lena, and Peppers.Figure 2
Set of biomedical cell images (left to right): Fluorescent cell 1, Fluorescent cell 2, Fluorescent cell 3, CIL 248, CIL 7437, and CIL 12293.


### 4.2. Recovery on Natural Images

In Table2, we provide systematic comparative results for all the regularization schemes on the natural images under different degradation conditions. As we can see from Table 2, regarding comparisons among the seven regularization schemes, the best PSNR results, on average, are achieved for the proposed MHOTV regularizer, the nonlocal TV regularizer provides better results than TV regularizer for most of the cases, and ISDTV regularizer outperforms TV and Hessian Spectral-norm regularizers for most of the cases, while the Hessian Frobenius-norm regularizer always performs similar results with ISDTV scheme, and the PSNR results of Hessian Nuclear-norm regularizer are almost always a little better than that of Hessian Spectral-norm and Frobenius-norm regularizers. Thus, we can conclude that our MHOTV-based regularization scheme is proven to be a valid approach for image recovery.Table 2
PSNR (dB) results of the seven regularization methods on natural images for three blurring kernels and two different noise levels.
Blur
Gaussian
Uniform
Motion
Gaussian
Uniform
Motion
BSNR (dB)
20
30
20
30
20
30
20
30
20
30
20
30
Barbara
Boat
TV
23.22
24.12
23.26
24.15
23.34
24.72
26.83
29.12
26.97
29.35
26.23
28.79
Nonlocal TV
23.42
24.21
23.45
24.22
23.28
25.02
26.94
29.25
27.07
29.46
25.87
28.67
Hessian Spectral
23.32
24.18
23.35
24.22
23.32
24.68
26.91
29.11
27.00
29.30
26.05
28.47
Hessian Frobenius
23.36
24.23
23.38
24.27
23.37
24.78
27.01
29.27
27.12
29.47
26.15
28.60
Hessian Nuclear
23.38
24.26
23.40
24.29
23.41
24.85
27.07
29.34
27.19
29.57
26.19
28.65
ISDTV
23.35
24.22
23.37
24.26
23.36
24.76
26.98
29.23
27.09
29.43
26.13
28.58
MHOTV
23.42
24.28
23.45
24.32
23.50
24.88
27.15
29.46
27.27
29.67
26.48
28.96
Face
Fingerprint
TV
26.13
28.83
26.11
28.73
24.97
28.43
23.47
26.37
23.31
26.26
22.43
25.30
Nonlocal TV
26.90
29.28
26.73
29.32
24.44
27.77
24.86
26.88
24.57
26.72
23.03
25.84
Hessian Spectral
26.84
29.61
26.74
29.49
25.70
29.23
24.68
27.07
24.53
27.04
23.75
26.30
Hessian Frobenius
26.98
29.77
26.89
29.69
25.79
29.39
24.79
27.13
24.62
27.11
23.79
26.36
Hessian Nuclear
27.03
29.79
26.96
29.69
25.83
29.46
24.75
27.07
24.60
27.06
23.80
26.32
ISDTV
26.95
29.82
26.87
29.67
25.78
29.44
24.74
27.09
24.58
27.06
23.76
26.32
MHOTV
27.12
30.01
27.00
29.84
25.86
29.58
24.69
27.27
24.43
27.16
23.45
26.23
Lena
Peppers
TV
25.49
27.40
25.54
27.56
25.59
28.37
26.21
29.35
26.54
29.49
24.77
27.19
Nonlocal TV
25.57
27.49
25.67
27.57
25.07
28.17
26.53
28.64
26.82
28.98
24.24
27.64
Hessian Spectral
25.34
27.10
25.41
27.20
25.25
28.16
26.19
28.98
26.46
29.08
25.12
27.97
Hessian Frobenius
25.44
27.22
25.50
27.32
25.39
28.33
26.35
29.14
26.64
29.25
25.23
28.15
Hessian Nuclear
25.51
27.26
25.56
27.37
25.51
28.43
26.48
29.23
26.77
29.37
25.31
28.16
ISDTV
25.42
27.19
25.48
27.29
25.36
28.29
26.29
29.09
26.59
29.21
25.19
28.11
MHOTV
25.78
27.60
25.83
27.73
25.81
28.73
26.69
29.65
27.01
29.75
25.56
28.66Furthermore, we show the representative Fingerprint, Face, and Peppers deblurring examples in Figures3, 4, and 5, which are degraded by Gaussian blurring, uniform blurring, and motion blurring with different Gaussian noise levels, respectively.Results of different methods on the Fingerprint image degraded by Gaussian blurring and noise corresponding to a BSNR level of 30 dB: (a) original image, (b) degraded image (PSNR = 18.56 dB), (c) TV solution (PSNR = 26.37 dB), (d) nonlocal TV solution (PSNR = 26.88 dB), (e) Hessian Frobenius-norm solution (PSNR = 27.13 dB), (f) Hessian Spectral-norm solution (PSNR = 27.07 dB), (g) Hessian Nuclear-norm solution (PSNR = 27.07 dB), (h) ISDTV solution (PSNR = 27.09 dB), and (i) MHOTV solution (PSNR = 27.27 dB).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Results of different methods on the Face image degraded by uniform blurring and noise corresponding to a BSNR level of 20 dB: (a) original image, (b) degraded image (PSNR = 21.77 dB), (c) TV solution (PSNR = 26.11 dB), (d) nonlocal TV solution (PSNR = 26.73 dB), (e) Hessian Frobenius-norm solution (PSNR = 26.89 dB), (f) Hessian Spectral-norm solution (PSNR = 26.74 dB), (g) Hessian Nuclear-norm solution (PSNR = 26.96 dB), (h) ISDTV solution (PSNR = 26.87 dB), and (i) MHOTV solution (PSNR = 27.00 dB).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Results of different methods on the Peppers image degraded by motion blurring and noise corresponding to a BSNR level of 30 dB: (a) original image, (b) degraded image (PSNR = 18.98 dB), (c) TV solution (PSNR = 27.19 dB), (d) nonlocal TV solution (PSNR = 27.64 dB), (e) Hessian Frobenius-norm solution (PSNR = 28.15 dB), (f) Hessian Spectral-norm solution (PSNR = 27.97 dB), (g) Hessian Nuclear-norm solution (PSNR = 28.16 dB), (h) ISDTV solution (PSNR = 28.11 dB), and (i) MHOTV solution (PSNR = 28.66 dB).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)From these figures, we can clearly see that the images recovered by TV regularization method suffer from heavy blocky artifacts appearing in the smoother regions of the image and result in some loss of details, the nonlocal TV regularization method preserves the sharp edges but smoothes some moderate gradient areas, the Hessian Spectral-norm, Hessian Frobenius-norm, Hessian Nuclear-norm, and ISDTV regularization methods are very competitive, and they can also preserve some sharp edges and minimize the staircase effect, while our MHOTV regularization method provides the best image quality which can effectively minimize the staircase effect and better preserve edges sharpness. For example, in Figure3, since the Fingerprint image is mainly composed of smooth transitions and ridge features, in order to highlight the differences, we present the details of the corresponding zooming versions of the region (marked by the black box in Figure 3(a)) for all the seven methods in Figures 3(c)–3(i). Comparing the results of these seven methods, we observe that the TV solution shown in Figure 3(c) results in some loss of fine ridge details and produces the staircase effect, the nonlocal TV solution shown in Figure 3(d) smoothes some ridge details, and the remaining five methods all preserve fine ridge details and edges better than TV and nonlocal TV methods, while our MHOTV solution provides the best performance in terms of both PSNR and image quality.


### 4.3. Recovery on Biomedical Cell Images

In this section, we mainly talk about the performance of the proposed MHOTV regularization scheme on biomedical images, which are closely related to the practical applications. To evaluate the effectiveness of our approach, we provide deblurring experiments for all the methods on the set of biomedical cell images, and their comparative results are given in Table3.Table 3
PSNR (dB) results of the seven regularization methods on biomedical cell images for two blurring kernels and two different noise levels.
Blur
Gaussian
Uniform
Gaussian
Uniform
BSNR (dB)
20
30
20
30
20
30
20
30
Fluorescent cell 1
Fluorescent cell 2
TV
27.52
29.20
27.62
29.29
28.21
30.12
28.28
30.35
Nonlocal TV
27.44
29.22
27.46
29.23
28.40
30.12
28.41
30.28
Hessian Spectral
27.92
29.63
28.00
29.74
28.61
30.38
28.63
30.61
Hessian Frobenius
27.97
29.73
28.03
29.83
28.72
30.51
28.77
30.69
Hessian Nuclear
27.97
29.75
28.04
29.88
28.74
30.57
28.82
30.77
ISDTV
27.96
29.70
28.01
29.80
28.69
30.48
28.74
30.70
MHOTV
27.97
29.67
28.02
29.76
28.72
30.55
28.78
30.78
Fluorescent cell 3
CIL 248
TV
27.15
29.01
27.15
29.24
27.70
29.65
27.69
29.79
Nonlocal TV
27.47
29.18
27.48
29.25
27.89
29.68
27.80
29.76
Hessian Spectral
27.47
29.25
27.49
29.47
28.23
30.00
28.19
30.10
Hessian Frobenius
27.53
29.35
27.56
29.55
28.26
30.10
28.19
30.18
Hessian Nuclear
27.53
29.37
27.60
29.56
28.23
30.08
28.17
30.18
ISDTV
27.52
29.32
27.55
29.53
28.25
30.11
28.20
30.21
MHOTV
27.67
29.50
27.68
29.73
28.35
30.25
28.30
30.38
CIL 7437
CIL 12293
TV
34.61
36.31
34.66
36.42
24.95
26.81
25.04
27.26
Nonlocal TV
34.63
36.43
34.57
36.43
24.86
26.58
24.95
26.88
Hessian Spectral
35.10
36.79
35.08
36.88
25.31
26.97
25.33
27.37
Hessian Frobenius
35.13
36.85
35.11
36.96
25.40
27.03
25.43
27.44
Hessian Nuclear
35.10
36.83
35.09
36.92
25.42
27.00
25.44
27.40
ISDTV
35.12
36.83
35.10
36.94
25.44
27.04
25.48
27.46
MHOTV
35.19
36.91
35.19
37.01
25.49
27.08
25.56
27.52With the PSNR comparisons in Table3, we can clearly observe that the proposed MHOTV regularization approach provides the best PSNR results for most of the cases, the ISDTV method always performs better than nonlocal TV, TV, and Hessian Spectral-norm methods, and the nonlocal TV method provides the better results than TV method for most of the cases, while the Hessian Frobenius-norm method always provides the comparatively similar results to ISDTV method, and the PSNR results of Hessian Nuclear-norm method are almost always a little better than that of Hessian Spectral-norm and Frobenius-norm methods. Meanwhile, we show the representative Fluorescent cell 3 and CIL 12293 deblurring examples in Figures 6 and 7, respectively. While carefully inspecting the recovered images shown in Figures 6 and 7, we obviously find that the advantage of our MHOTV approach consists in better preserving the ridges and filament-like features sharpness without introducing severe block artifacts, and the nonlocal TV, ISDTV, and the three Hessian-norm regularization methods can preserve some fine ridge features and alleviate the staircase effect, as opposed to TV which introduces the heavy block artifacts which tends to produce the oversmoothed image ridge and filament-like features. Therefore, we conclude that our MHOTV regularization approach can be effectively and successfully applied in the process of biomedical imaging.Results of different methods on the Fluorescent cell 3 image degraded by Gaussian blurring and noise corresponding to a BSNR level of 20 dB: (a) original image, (b) degraded image (PSNR = 24.40 dB), (c) TV solution (PSNR = 27.15 dB), (d) nonlocal TV solution (PSNR = 27.47 dB), (e) Hessian Frobenius-norm solution (PSNR = 27.53 dB), (f) Hessian Spectral-norm solution (PSNR = 27.47 dB), (g) Hessian Nuclear-norm solution (PSNR = 27.53 dB), (h) ISDTV solution (PSNR = 27.52 dB), and (i) MHOTV solution (PSNR = 27.67 dB).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Results of different methods on the CIL 12293 image degraded by uniform blurring and noise corresponding to a BSNR level of 20 dB: (a) original image, (b) degraded image (PSNR = 21.00 dB), (c) TV solution (PSNR = 25.04 dB), (d) nonlocal TV solution (PSNR = 24.95 dB), (e) Hessian Frobenius-norm solution (PSNR = 25.43 dB), (f) Hessian Spectral-norm solution (PSNR = 25.33 dB), (g) Hessian Nuclear-norm solution (PSNR = 25.44 dB), (h) ISDTV solution (PSNR = 25.48 dB), and (i) MHOTV solution (PSNR = 25.56 dB).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)


## 4.1. Experiment Setting

For the image deblurring experiments, our test images contain six natural images and six biomedical cell images which are shown in Figures1 and 2, respectively, where the cell images are part of the biomedical image database [34], and they are all converted to grayscale. In our experiments, the intensities of all the test images are initially normalized to the range[
0,1
], and we use three kinds of point spread function (PSF) to produce the blurred images, which are a9
×
9Gaussian PSF of standard deviationσ
b
=
6, a9
×
9uniform one, and a19
×
19motion one, respectively. Furthermore, the Gaussian noise of two different levels corresponding to a blurred SNR (BSNR) of{
20,30
} dB is added to generate the eventually degraded images, and the BSNR is defined as
(36)
BSNR
=
var
(
A
f
)
σ
n
2
,
wherevar
(
A
f
)is the variance of blurred image andσ
nis the standard deviation of Gaussian noise. In addition, we constrain the intensities of image recovered by all the seven methods to lie in the convex closed setS
=
{
f
∈
ℜ
N
∣
f
i
∈
[
0
,
1
]
,
∀
i
=
1
,
…
,
N
}. For the sake of fair comparisons, the regularization parameterλfor each regularization method is chosen to give the best PSNR performance. Finally, we setε
=
1
0
-
5,M
a
x
i
t
e
r
=
100,ε
0
=
1
0
-
3
,andI
n
n
e
r
i
t
e
r
=
10in our experiments.Figure 1
Set of natural images (left to right): Barbara, Boat, Face, Fingerprint, Lena, and Peppers.Figure 2
Set of biomedical cell images (left to right): Fluorescent cell 1, Fluorescent cell 2, Fluorescent cell 3, CIL 248, CIL 7437, and CIL 12293.


## 4.2. Recovery on Natural Images

In Table2, we provide systematic comparative results for all the regularization schemes on the natural images under different degradation conditions. As we can see from Table 2, regarding comparisons among the seven regularization schemes, the best PSNR results, on average, are achieved for the proposed MHOTV regularizer, the nonlocal TV regularizer provides better results than TV regularizer for most of the cases, and ISDTV regularizer outperforms TV and Hessian Spectral-norm regularizers for most of the cases, while the Hessian Frobenius-norm regularizer always performs similar results with ISDTV scheme, and the PSNR results of Hessian Nuclear-norm regularizer are almost always a little better than that of Hessian Spectral-norm and Frobenius-norm regularizers. Thus, we can conclude that our MHOTV-based regularization scheme is proven to be a valid approach for image recovery.Table 2
PSNR (dB) results of the seven regularization methods on natural images for three blurring kernels and two different noise levels.
Blur
Gaussian
Uniform
Motion
Gaussian
Uniform
Motion
BSNR (dB)
20
30
20
30
20
30
20
30
20
30
20
30
Barbara
Boat
TV
23.22
24.12
23.26
24.15
23.34
24.72
26.83
29.12
26.97
29.35
26.23
28.79
Nonlocal TV
23.42
24.21
23.45
24.22
23.28
25.02
26.94
29.25
27.07
29.46
25.87
28.67
Hessian Spectral
23.32
24.18
23.35
24.22
23.32
24.68
26.91
29.11
27.00
29.30
26.05
28.47
Hessian Frobenius
23.36
24.23
23.38
24.27
23.37
24.78
27.01
29.27
27.12
29.47
26.15
28.60
Hessian Nuclear
23.38
24.26
23.40
24.29
23.41
24.85
27.07
29.34
27.19
29.57
26.19
28.65
ISDTV
23.35
24.22
23.37
24.26
23.36
24.76
26.98
29.23
27.09
29.43
26.13
28.58
MHOTV
23.42
24.28
23.45
24.32
23.50
24.88
27.15
29.46
27.27
29.67
26.48
28.96
Face
Fingerprint
TV
26.13
28.83
26.11
28.73
24.97
28.43
23.47
26.37
23.31
26.26
22.43
25.30
Nonlocal TV
26.90
29.28
26.73
29.32
24.44
27.77
24.86
26.88
24.57
26.72
23.03
25.84
Hessian Spectral
26.84
29.61
26.74
29.49
25.70
29.23
24.68
27.07
24.53
27.04
23.75
26.30
Hessian Frobenius
26.98
29.77
26.89
29.69
25.79
29.39
24.79
27.13
24.62
27.11
23.79
26.36
Hessian Nuclear
27.03
29.79
26.96
29.69
25.83
29.46
24.75
27.07
24.60
27.06
23.80
26.32
ISDTV
26.95
29.82
26.87
29.67
25.78
29.44
24.74
27.09
24.58
27.06
23.76
26.32
MHOTV
27.12
30.01
27.00
29.84
25.86
29.58
24.69
27.27
24.43
27.16
23.45
26.23
Lena
Peppers
TV
25.49
27.40
25.54
27.56
25.59
28.37
26.21
29.35
26.54
29.49
24.77
27.19
Nonlocal TV
25.57
27.49
25.67
27.57
25.07
28.17
26.53
28.64
26.82
28.98
24.24
27.64
Hessian Spectral
25.34
27.10
25.41
27.20
25.25
28.16
26.19
28.98
26.46
29.08
25.12
27.97
Hessian Frobenius
25.44
27.22
25.50
27.32
25.39
28.33
26.35
29.14
26.64
29.25
25.23
28.15
Hessian Nuclear
25.51
27.26
25.56
27.37
25.51
28.43
26.48
29.23
26.77
29.37
25.31
28.16
ISDTV
25.42
27.19
25.48
27.29
25.36
28.29
26.29
29.09
26.59
29.21
25.19
28.11
MHOTV
25.78
27.60
25.83
27.73
25.81
28.73
26.69
29.65
27.01
29.75
25.56
28.66Furthermore, we show the representative Fingerprint, Face, and Peppers deblurring examples in Figures3, 4, and 5, which are degraded by Gaussian blurring, uniform blurring, and motion blurring with different Gaussian noise levels, respectively.Results of different methods on the Fingerprint image degraded by Gaussian blurring and noise corresponding to a BSNR level of 30 dB: (a) original image, (b) degraded image (PSNR = 18.56 dB), (c) TV solution (PSNR = 26.37 dB), (d) nonlocal TV solution (PSNR = 26.88 dB), (e) Hessian Frobenius-norm solution (PSNR = 27.13 dB), (f) Hessian Spectral-norm solution (PSNR = 27.07 dB), (g) Hessian Nuclear-norm solution (PSNR = 27.07 dB), (h) ISDTV solution (PSNR = 27.09 dB), and (i) MHOTV solution (PSNR = 27.27 dB).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Results of different methods on the Face image degraded by uniform blurring and noise corresponding to a BSNR level of 20 dB: (a) original image, (b) degraded image (PSNR = 21.77 dB), (c) TV solution (PSNR = 26.11 dB), (d) nonlocal TV solution (PSNR = 26.73 dB), (e) Hessian Frobenius-norm solution (PSNR = 26.89 dB), (f) Hessian Spectral-norm solution (PSNR = 26.74 dB), (g) Hessian Nuclear-norm solution (PSNR = 26.96 dB), (h) ISDTV solution (PSNR = 26.87 dB), and (i) MHOTV solution (PSNR = 27.00 dB).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Results of different methods on the Peppers image degraded by motion blurring and noise corresponding to a BSNR level of 30 dB: (a) original image, (b) degraded image (PSNR = 18.98 dB), (c) TV solution (PSNR = 27.19 dB), (d) nonlocal TV solution (PSNR = 27.64 dB), (e) Hessian Frobenius-norm solution (PSNR = 28.15 dB), (f) Hessian Spectral-norm solution (PSNR = 27.97 dB), (g) Hessian Nuclear-norm solution (PSNR = 28.16 dB), (h) ISDTV solution (PSNR = 28.11 dB), and (i) MHOTV solution (PSNR = 28.66 dB).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)From these figures, we can clearly see that the images recovered by TV regularization method suffer from heavy blocky artifacts appearing in the smoother regions of the image and result in some loss of details, the nonlocal TV regularization method preserves the sharp edges but smoothes some moderate gradient areas, the Hessian Spectral-norm, Hessian Frobenius-norm, Hessian Nuclear-norm, and ISDTV regularization methods are very competitive, and they can also preserve some sharp edges and minimize the staircase effect, while our MHOTV regularization method provides the best image quality which can effectively minimize the staircase effect and better preserve edges sharpness. For example, in Figure3, since the Fingerprint image is mainly composed of smooth transitions and ridge features, in order to highlight the differences, we present the details of the corresponding zooming versions of the region (marked by the black box in Figure 3(a)) for all the seven methods in Figures 3(c)–3(i). Comparing the results of these seven methods, we observe that the TV solution shown in Figure 3(c) results in some loss of fine ridge details and produces the staircase effect, the nonlocal TV solution shown in Figure 3(d) smoothes some ridge details, and the remaining five methods all preserve fine ridge details and edges better than TV and nonlocal TV methods, while our MHOTV solution provides the best performance in terms of both PSNR and image quality.


## 4.3. Recovery on Biomedical Cell Images

In this section, we mainly talk about the performance of the proposed MHOTV regularization scheme on biomedical images, which are closely related to the practical applications. To evaluate the effectiveness of our approach, we provide deblurring experiments for all the methods on the set of biomedical cell images, and their comparative results are given in Table3.Table 3
PSNR (dB) results of the seven regularization methods on biomedical cell images for two blurring kernels and two different noise levels.
Blur
Gaussian
Uniform
Gaussian
Uniform
BSNR (dB)
20
30
20
30
20
30
20
30
Fluorescent cell 1
Fluorescent cell 2
TV
27.52
29.20
27.62
29.29
28.21
30.12
28.28
30.35
Nonlocal TV
27.44
29.22
27.46
29.23
28.40
30.12
28.41
30.28
Hessian Spectral
27.92
29.63
28.00
29.74
28.61
30.38
28.63
30.61
Hessian Frobenius
27.97
29.73
28.03
29.83
28.72
30.51
28.77
30.69
Hessian Nuclear
27.97
29.75
28.04
29.88
28.74
30.57
28.82
30.77
ISDTV
27.96
29.70
28.01
29.80
28.69
30.48
28.74
30.70
MHOTV
27.97
29.67
28.02
29.76
28.72
30.55
28.78
30.78
Fluorescent cell 3
CIL 248
TV
27.15
29.01
27.15
29.24
27.70
29.65
27.69
29.79
Nonlocal TV
27.47
29.18
27.48
29.25
27.89
29.68
27.80
29.76
Hessian Spectral
27.47
29.25
27.49
29.47
28.23
30.00
28.19
30.10
Hessian Frobenius
27.53
29.35
27.56
29.55
28.26
30.10
28.19
30.18
Hessian Nuclear
27.53
29.37
27.60
29.56
28.23
30.08
28.17
30.18
ISDTV
27.52
29.32
27.55
29.53
28.25
30.11
28.20
30.21
MHOTV
27.67
29.50
27.68
29.73
28.35
30.25
28.30
30.38
CIL 7437
CIL 12293
TV
34.61
36.31
34.66
36.42
24.95
26.81
25.04
27.26
Nonlocal TV
34.63
36.43
34.57
36.43
24.86
26.58
24.95
26.88
Hessian Spectral
35.10
36.79
35.08
36.88
25.31
26.97
25.33
27.37
Hessian Frobenius
35.13
36.85
35.11
36.96
25.40
27.03
25.43
27.44
Hessian Nuclear
35.10
36.83
35.09
36.92
25.42
27.00
25.44
27.40
ISDTV
35.12
36.83
35.10
36.94
25.44
27.04
25.48
27.46
MHOTV
35.19
36.91
35.19
37.01
25.49
27.08
25.56
27.52With the PSNR comparisons in Table3, we can clearly observe that the proposed MHOTV regularization approach provides the best PSNR results for most of the cases, the ISDTV method always performs better than nonlocal TV, TV, and Hessian Spectral-norm methods, and the nonlocal TV method provides the better results than TV method for most of the cases, while the Hessian Frobenius-norm method always provides the comparatively similar results to ISDTV method, and the PSNR results of Hessian Nuclear-norm method are almost always a little better than that of Hessian Spectral-norm and Frobenius-norm methods. Meanwhile, we show the representative Fluorescent cell 3 and CIL 12293 deblurring examples in Figures 6 and 7, respectively. While carefully inspecting the recovered images shown in Figures 6 and 7, we obviously find that the advantage of our MHOTV approach consists in better preserving the ridges and filament-like features sharpness without introducing severe block artifacts, and the nonlocal TV, ISDTV, and the three Hessian-norm regularization methods can preserve some fine ridge features and alleviate the staircase effect, as opposed to TV which introduces the heavy block artifacts which tends to produce the oversmoothed image ridge and filament-like features. Therefore, we conclude that our MHOTV regularization approach can be effectively and successfully applied in the process of biomedical imaging.Results of different methods on the Fluorescent cell 3 image degraded by Gaussian blurring and noise corresponding to a BSNR level of 20 dB: (a) original image, (b) degraded image (PSNR = 24.40 dB), (c) TV solution (PSNR = 27.15 dB), (d) nonlocal TV solution (PSNR = 27.47 dB), (e) Hessian Frobenius-norm solution (PSNR = 27.53 dB), (f) Hessian Spectral-norm solution (PSNR = 27.47 dB), (g) Hessian Nuclear-norm solution (PSNR = 27.53 dB), (h) ISDTV solution (PSNR = 27.52 dB), and (i) MHOTV solution (PSNR = 27.67 dB).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Results of different methods on the CIL 12293 image degraded by uniform blurring and noise corresponding to a BSNR level of 20 dB: (a) original image, (b) degraded image (PSNR = 21.00 dB), (c) TV solution (PSNR = 25.04 dB), (d) nonlocal TV solution (PSNR = 24.95 dB), (e) Hessian Frobenius-norm solution (PSNR = 25.43 dB), (f) Hessian Spectral-norm solution (PSNR = 25.33 dB), (g) Hessian Nuclear-norm solution (PSNR = 25.44 dB), (h) ISDTV solution (PSNR = 25.48 dB), and (i) MHOTV solution (PSNR = 25.56 dB).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)


## 5. Conclusion

In this paper, we have proposed a novel MHOTV regularizer which combines the first and second degree derivatives of the image. Within spectral decomposition framework, we reformulate both the ISDTV and MHOTV regularizers as weightedL
1-L
2mixed norms of image derivatives. Furthermore, we present a unified variational model based on the ISDTV and MHOTV regularizers for image recovery. Based on the equivalent formulations of ISDTV and MHOTV regularizers, we then introduce an efficient fast projected gradient algorithm combined with monotone fast iterative shrinkage thresholding (FPG-MFISTA) to solve the corresponding optimization problems under majorization-minimization framework. Finally, the performance of the proposed MHOTV regularization scheme was assessed by comparing the experimental results on both natural images and biomedical cell images against various state-of-the-art methods. Specifically, the proposed MHOTV regularization scheme performs improved results and can effectively minimize the staircase effect that is commonly met in TV scheme while better preserving edges sharpness.



---

*Source: 924686-2014-02-27.xml*