An adaptive model combining a total variation filter and a fractional-order filter for image restoration

Abstract

In this paper, we propose a model that combines a total variation filter with a fractional-order filter, which can unite the advantages of the two filters, and has a remarkable effect in the protection of image edges and texture details; simultaneously, the proposed model can eliminate the staircase effect. In addition, the model improves the PSNR compared with the total variation filter and the fractional-order filter when removing noise. Zhu and Chan presented the primal-dual hybrid gradient algorithm and proved that it is effective for the total variation filter. On the basis of their work, we employ the primal-dual hybrid gradient algorithm to solve the combined model in this article. The final experimental results show that the new model and algorithm are effective for image restoration.

Introduction

The essence of image restoration problem is to find an original image u by solving an inverse problem. The regularization method, proposed by Tikhonov and Arsenin in 1977, ¹ has been proven to be one of the efficient methods for the inverse problem. Since then, the image restoration technology based on variational calculus has been extremely successful. Tikhonov and Arsenin first introduced Sobolev space $W^{\text{1}�\text{2}}\left(\Omega \right)$ into the regularization term, and they proposed the restoration model as follows

$$\underset{u\in W^{\text{1},\text{2}}(\Omega )}{\min }\lambda {\displaystyle {\int }_{\Omega }|\nabla u{|}^2}\text{d}x+\frac{1}{2}\left|\right|u-u_0|{|}^2$$

(1)where $\Omega \in R^n$ be an open bounded domain with a Lipschitzian boundary, u₀ represents the degraded image. The first term is regularization term which can guarantee the well-posedness of solution while removing noise. The second term is fidelity term which can fit statistical distribution of the image. The regularization parameter $\lambda$ adjusts the weight of regularization term and fidelity term. In 1992, Rudin et al. ² regarded the image as a function of the bounded variation functions space BV (Ω) and they proposed the total variation model (also called ROF model) for image restoration, and it is depicted as

$$\underset{u\in \mathit{BV}(\Omega )\cap L^{\text{2}}(\Omega )}{\min }{\displaystyle {\int }_{\Omega }\left|\mathit{Du}\right|}+\frac{\lambda }{2}\left|\right|u-u_0|{|}^2$$

(2)where

$${\displaystyle {\int }_{\Omega }\left|\mathit{Du}\right|}=\sup \left\{{\displaystyle {\int }_{\Omega }u\cdot \mathit{div}\phi \mathit{dx};\phi =({\phi }_1,{\phi }_2,\cdots {\phi }_N)\in C_0^1{(\Omega )}^N,|\phi {|}_{L^{\infty }(\Omega )}\le 1}\right\}$$

is the total variation of u,

$$\mathit{BV}(\Omega )=\left\{u\in L^1(\Omega );{\displaystyle {\int }_{\Omega }\left|\mathit{Du}\right|}<\infty \right\}.$$

If $u\in W^{1,1}(\Omega )$ , the total variation model can be written as

$$\underset{u\in W^{1,1}(\Omega )\cap L^{\text{2}}(\Omega )}{\min }{\displaystyle {\int }_{\Omega }|\nabla u|}+\frac{\lambda }{2}\left|\right|u-u_0|{|}^2$$

(3)where $\nabla u$ is the gradient of $u$ . The Euler–Lagrange equation of equation (3) is an anisotropic diffusion second-order partial differential equation (PDE).

Owing to the bounded variation functions, space contains discontinuous functions which conform to the characteristics of image edges, and the total variation model performs very well for preserving image edges. However, it favors piecewise constant solution in image smooth regions, which will cause staircase effect. For the sake of solving these problems, it was suggested that we could improve the model on the basis of the original variational model, such as generalized total variation mode, ³ adaptive total variation model ⁴ and mixed model with L1/L2 fidelity terms. ⁵ To a certain extent, these models can suppress staircase effect, but it is difficult to eliminate it fundamentally because of the intrinsic characteristics of solution space BV (Ω). In order to overcome this drawback, high-order PDE filter technique was widely researched. In 2003, Lysaker et al. ⁶ proposed a classical LLT model, and its expression is as follows

$$\underset{u\in B{V}^{\text{2}}(\Omega )\cap L^{\text{2}}(\Omega )}{\min }{\displaystyle {\int }_{\Omega }\left|D^{\text{2}}u\right|}+\frac{\lambda }{2}\left|\right|u-u_0|{|}^2$$

(4)where

$${\displaystyle {\int }_{\Omega }\left|D^{\text{2}}u\right|}=\sup \left\{{\displaystyle {\int }_{\Omega }{\displaystyle \sum _{i,j=1}^{n}u{\partial }_j{\partial }_i{\phi }_{i,j}}dx;\phi \in C_0^2(\Omega ,R^{n\times n}),|\phi {|}_{L_{\infty }}\le 1}\right\}$$

And $B{V}^2(\Omega )=\left\{u\in L^1(\Omega );{\displaystyle {\int }_{\Omega }\left|D^{2}u\right|}<\infty \right\}$ . If $u\in W^{\text{2},1}(\Omega )$ , the LLT model becomes

$$\underset{u\in W^{2,1}(\Omega )\cap L^{\text{2}}(\Omega )}{\min }{\displaystyle {\int }_{\Omega }\left|{\nabla }^{\text{2}}u\right|}+\frac{\lambda }{2}\left|\right|u-u_0|{|}^2$$

(5)where $\text{|}{\nabla }^{2}u\text{|=}\sqrt{u_{\mathit{xx}}{}^2+u_{\mathit{xy}}{}^2+u_{\mathit{yx}}{}^2+u_{\mathit{yy}}{}^2}$ . The Euler–Lagrange equation of LLT model is a fourth-order PDE. The LLT model can remove the staircase effect, which satisfies the human visual effect. But because the diffusion speed of the fourth-order PDE is too fast, it causes image blurring. For the purpose of conquering the staircase effect and avoiding the blurring effect simultaneously, Li and Shen presented an approach of combining ROF model and LLT model. ⁷ The numerical results showed that the combined model not only overcame their deficiencies but also combined the advantages of each model.

In recent years, with the research and development of fractional calculus, the fractional calculus has been widely used in biology, physics, fluid mechanics and image processing field. The fractional calculus theory, especially the fractional differential, has achieved some results in image restoration and edges detection.^8–10 At present, the Grünewald-Letnikov fractional-order derivative definition was often used in image processing. ¹¹ Its definition is as follows

$${}^{\mathit{GL}}D\alpha u(t)=\underset{h\to 0}{\lim }\frac{{\displaystyle {\sum }_{k=0}^K{(-1)}^k{C}_k^{\alpha }}u(t-\mathit{kh})}{h^{\alpha }}$$

where $C_k^{\alpha }=\frac{\Gamma (\alpha +1)}{\Gamma (k+1)\Gamma (\alpha -k+1)}$ , $\Gamma (\cdot )$ is the gamma function. In Jun and Zhi, ¹² He et al. ¹³ and Chen and Sun ¹⁴ the researchers generalized the integer-order variation model to the fractional-order variation model for image restoration. The fractional-order variation model (abbreviated as FV model) is described as follows

$$\underset{u}{\min }{\displaystyle {\int }_{\Omega }\left|{\nabla }^{\alpha }u\right|}+\frac{\lambda }{2}\left|\right|u-u_0|{|}^2$$

(6)where $\alpha \in \text{(1,2)}$ , ${\nabla }^{\alpha }u$ is the fractional-order gradient of u and $\text{|}{\nabla }^{\alpha }u\text{|=}\sqrt{{\left({\nabla }_x^{\alpha }u\right)}^2+{\left({\nabla }_y^{\alpha }u\right)}^2}$ . The numerical results proved that the fractional-order variation model had higher PSNR and performed well in retaining texture details and eliminating the staircase effect.

In this paper, based on the previous work, we present a new model for image restoration and apply the primal-dual hybrid gradient algorithm to solve the model. ¹⁵ The rest of the paper is organized as follows. In the next section, we introduce the proposed model and its Euler–Lagrange equation in detail; then , we derive the primal-dual hybrid gradient algorithm for the proposed model and give its convergence theorem. In the subsequent section, we give the numerical experiments and results to demonstrate the effectiveness of the proposed model. Finally, we conclude this paper with conclusion.

The proposed model

In this section, we mainly study the adaptivity of the new model and the Euler–Lagrange equation of it. At first, we give a detailed introduction for the proposed model, as follows

$$\underset{u}{\min }{\displaystyle {\int }_{\Omega }\left(1-g\right)|\nabla u|}+{\displaystyle {\int }_{\Omega }g\left|{\nabla }^{\alpha }u\right|}+\frac{\lambda }{2}\left|\right|u-u_0|{|}^2$$

(7)where g and 1−g are the weight coefficients. Specially, it is the ROF model when g = 0 and it is the FV model when g = 1.

Strong et al. ¹⁶ proposed an adaptive total variation model for image restoration. The adaptive weight function $g(x)$ can control the diffusion speed, and its expression is given as follows

$$g(x)=\frac{1}{1+K|\nabla G_{\sigma }\ast u_0{|}^2}$$

(8)where $K\ge \text{0}$ , the Gaussian kernel function $G_{\sigma }=\frac{1}{2\pi {\sigma }^2}\exp \left(-\frac{x^2}{2{\sigma }^2}\right)$ and $0<g(x)\le \text{1}$ . Here, the image can be preprocessed with Gaussian filter to remove a small part of the noise in advance, which can reduce the error of mistaking the noise as a false edge and decrease the staircase effect to some extent, thus making our image processing more accurate. In Liu et al. and Tian etal.,^17,18 in the same way, the researchers presented the adaptive LLT model and the adaptive FV model for image restoration. Their experimental results indicated that the model with adaptive weight function can improve image restoration quality more effectively than the original model. In equation (7), inspired by the above models, we use $1-g(x)$ and $g(x)$ as weight functions of the ROF model and the FV model. By analyzing equations (7) and (8), we can obtain: (1) For given K, $|\nabla G_{\sigma }*u_0|$ is small and $g(x)\to \text{1}$ in image flat regions, the FV model plays the main role; therefore, the proposed model can theoretically suppress the staircase effect and preserve more texture details. (2) $|\nabla G_{\sigma }*u_0|$ is large and $g\to \text{0}$ in image edges regions, the ROF model plays the main role, and hence the model theoretically can protect image edges well.

Mathematically, for the variational problem, the general approach is to solve the corresponding Euler–Lagrange equation^19,20 or apply optimization algorithm. From the perspective of variational calculus, we derive the Euler–Lagrange equation of equation (7) with the variation method.

Let

$$E(u)={\displaystyle {\int }_{\Omega }\left(1-g\right)|\nabla u|}+{\displaystyle {\int }_{\Omega }g\left|{\nabla }^{\alpha }u\right|}+\frac{\lambda }{2}\left|\right|u-u_0|{|}^2$$

(9)

Then we construct an auxiliary function $F(\epsilon )$ , for any real number $\epsilon$ , function $\phi \in C_0^{\infty }(\Omega )$

$$\begin{array}{l}F(\epsilon )=E(u+\epsilon \phi )={\displaystyle {\int }_{\Omega }\left(1-g\right)|\nabla (u+\epsilon \phi )|}\\ +{\displaystyle {\int }_{\Omega }g|{\nabla }^{\alpha }(u+\epsilon \phi )|}+\frac{\lambda }{2}\left|\right|u+\epsilon \phi -u_0|{|}^2\end{array}$$

(10)

By the function extremum theorem, we have $F\hspace{0.17em}\prime {\text{|}}_{\epsilon \text{=0}}=0$ , that is

$$\begin{array}{l}{\displaystyle {\int }_{\Omega }\left(1-g\right)\left(\frac{{\nabla }_{x}u}{|\nabla u|}\cdot {\nabla }_x\phi +\frac{{\nabla }_{y}u}{|\nabla u|}\cdot {\nabla }_y\phi \right)}\\ +g\left(\frac{{\nabla }_x^{\alpha }u}{\left|{\nabla }^{\alpha }u\right|}\cdot {\nabla }_x^{\alpha }\phi +\frac{{\nabla }_y^{\alpha }u}{\left|{\nabla }^{\alpha }u\right|}\cdot {\nabla }_y^{\alpha }\phi \right)+\lambda (u-u_0)\phi \text{d}x\text{d}y=0\end{array}$$

Introducing adjoint operator ${\nabla }_x^{*}$ , ${\nabla }_y^{*}$ , ${({\nabla }_x^{\alpha })}^{*}$ , ${({\nabla }_y^{\alpha })}^{*}$ , the equation becomes

$$\begin{array}{l}{\displaystyle {\int }_{\Omega }\left(1-g\right)\left({\nabla }_x^{*}\left(\frac{{\nabla }_{x}u}{|\nabla u|}\right)+{\nabla }_y^{*}\left(\frac{{\nabla }_{y}u}{|\nabla u|}\right)\right)}\\ \cdot \phi +g({({\nabla }_x^{\alpha })}^{*}\left(\frac{{\nabla }_x^{\alpha }u}{\left|{\nabla }^{\alpha }u\right|}\right)+{({\nabla }_y^{\alpha })}^{*}\left(\frac{{\nabla }_y^{\alpha }u}{\left|{\nabla }^{\alpha }u\right|}\right))\\ \cdot \phi +\lambda (u-u_0)\phi \text{d}x\text{d}y=0\end{array}$$

By the variation method preparatory theorem, we get the Euler–Lagrange equation of equation (7) as follows

$$\begin{array}{l}\left(1-g\right)\left({\nabla }_x^{*}\left(\frac{{\nabla }_{x}u}{|\nabla u|}\right)+{\nabla }_y^{*}\left(\frac{{\nabla }_{y}u}{|\nabla u|}\right)\right)\\ +g({({\nabla }_x^{\alpha })}^{*}\left(\frac{{\nabla }_x^{\alpha }u}{\left|{\nabla }^{\alpha }u\right|}\right)\\ +{({\nabla }_y^{\alpha })}^{*}\left(\frac{{\nabla }_y^{\alpha }u}{\left|{\nabla }^{\alpha }u\right|}\right))\\ +\lambda (u-u_0)=0\end{array}$$

(11)where ${\nabla }_x^{*}=-{\nabla }_x$ , ${\nabla }_y^{*}=-{\nabla }_y$ . For fraction-order adjoint operator ${({\nabla }_x^{\alpha })}^{*}$ and ${({\nabla }_y^{\alpha })}^{*}$ , by fractional Fourier Transform formulas

$$\widehat{{\nabla }_x^{\alpha }h}={(j2\pi {\omega }_1)}^{\alpha }\cdot \widehat{h},\widehat{{\nabla }_y^{\alpha }h}={(j2\pi {\omega }_2)}^{\alpha }\cdot \widehat{h}$$

And Parseval identical formula

$${\displaystyle {\int }_{\Omega }f\cdot h\text{d}x\text{d}y}=\frac{1}{2\pi }{\displaystyle {\int }_{\Omega }\widehat{f}\cdot \overline{\widehat{h}}\text{d}{\omega }_1\text{d}{\omega }_2}$$

We have ${({\nabla }_x^{\alpha })}^{*}={(-1)}^{\alpha }\cdot {\nabla }_x^{\alpha }$ , ${({\nabla }_y^{\alpha })}^{*}={(-1)}^{\alpha }\cdot {\nabla }_y^{\alpha }$ .

Now, we plug the above deductions into equation (11) and get the Euler–Lagrange equation as follows

$$\begin{array}{l}-(1-g)\nabla \cdot \left(\frac{\nabla u}{|\nabla u|}\right)+g\cdot {(-1)}^{\alpha }\left({\nabla }_x^{\alpha }\left(\frac{{\nabla }_x^{\alpha }u}{\left|{\nabla }^{\alpha }u\right|}\right)+{\nabla }_y^{\alpha }\left(\frac{{\nabla }_y^{\alpha }u}{\left|{\nabla }^{\alpha }u\right|}\right)\right)\\ \text{+}\lambda (u-u_0)=0\end{array}$$

(12)where we replace $\text{|}{\nabla }^{\alpha }u\text{|}$ by $\text{|}{\nabla }^{\alpha }u{\text{|}}_{\epsilon }=\sqrt{\text{|}{\nabla }^{\alpha }u{\text{|}}^2+\epsilon }$ to avoid singularities of equation (12). For fractional PDE, we introduce artificial time $t$ and employ the steepest descent method to solve it.

$$\{\begin{array}{l}\frac{\partial u}{\partial t}=(1-g)\nabla \cdot \left(\frac{\nabla u}{|\nabla u{|}_{\epsilon }}\right)\hfill \\ -g\cdot {(-1)}^{\alpha }\left({\nabla }_x^{\alpha }\left(\frac{{\nabla }_x^{\alpha }u}{|{\nabla }^{\alpha }u{|}_{\epsilon }}\right)+{\nabla }_y^{\alpha }\left(\frac{{\nabla }_y^{\alpha }u}{|{\nabla }^{\alpha }u{|}_{\epsilon }}\right)\right)-\lambda (u-u_0),\\ u(x,y,0)=u_0\hfill \end{array}$$

(13)

Primal-dual hybrid gradient algorithm

In the last two decades, researchers have proposed many optimization algorithms for image restoration, such as augmented Lagrangian method, ²¹ alternating minimization algorithm, ²² duality-based algorithm, ²³ primal-dual algorithm, ²⁴ etc. In Bonettini and Ruggiero, ²⁵ the author proposed a primal-dual hybrid gradient algorithm for ROF model and proved its convergence. Zhu and Chan ¹⁵ simplified the primal-dual hybrid gradient algorithm by two ingenious treatments. Firstly, they converted the gradient of image $\nabla u$ into a multiplication of matrix and vector. Secondly, they introduced the dual variable through a conclusion of Cauchy–Schwarz inequality instead of the definition of convex conjugate operator. In this section, we will focus our study on the primal-dual hybrid gradient algorithm for the variational problem of equation (7).

Notation and preparation of algorithm

In this section, we consider the model in discrete setting. To simplify, we assume that the image is two-dimension matrices of size $n\times n$ . By discrete treatment, equation (7) becomes

$$\underset{u\in X}{\min }{\displaystyle \sum _{\text{1}\le i,j\le n}(1-g_{i,j})|\nabla u_{i,j}|}+{\displaystyle \sum _{\text{1}\le i,j\le n}{g}_{i,j}\left|{\nabla }^{\alpha }{u}_{i,j}\right|}+\frac{\lambda }{2}\left|\right|u-u_0|{|}^2$$

(14)where $X=R^{n\times n}$ , gradient $\nabla u$ and fractional-order gradient ${\nabla }^{\alpha }u$ are vectors in $X\times X$ , ${\nabla }^{\alpha }{u}_{i,j}$ is the $(i,j)$ element of ${\nabla }^{\alpha }u$ , ${\nabla }^{\alpha }{u}_{i,j}={\left({\nabla }_x^{\alpha }{u}_{i,j},{\nabla }_y^{\alpha }{u}_{i,j}\right)}^{\text{T}}$ , $\left|{\nabla }^{\alpha }{u}_{i,j}\right|=\sqrt{{\left({\nabla }_x^{\alpha }{u}_{i,j}\right)}^2+{\left({\nabla }_y^{\alpha }{u}_{i,j}\right)}^2}$ . $g_{i,j}$ is the $(i,j)$ element of the weight coefficient matrix $g$ , $g_{i,j}=\frac{1}{1+K|\nabla {(G_{\sigma }*u_0)}_{i,j}{|}^2}.$

The fractional-order gradient ${\nabla }^{\alpha }{u}_{i,j}$ $(i,j=1,2,\cdots ,n)$ difference form ( $h=1$ ) is as follows

$${\nabla }^{\alpha }{u}_{i,j}=\{\begin{array}{lllll}{\left({\displaystyle \sum _{k=0}^K{(-1)}^k{C}_k^{\alpha }{u}_{i-k,j}},{\displaystyle \sum _{k=0}^K{(-1)}^k{C}_k^{\alpha }{u}_{i,j-k}}\right)}^{\text{T}}\hfill & \text{ if }\hfill & i>1\hfill & \text{ and }\hfill & j>1\hfill \\ {\left(0,{\displaystyle \sum _{k=0}^K{(-1)}^k{C}_k^{\alpha }{u}_{i,j-k}}\right)}^{\text{T}}\hfill & \text{ if }\hfill & i=1\hfill & \text{ and }\hfill & j>1\hfill \\ {\left({\displaystyle \sum _{k=0}^K{(-1)}^k{C}_k^{\alpha }{u}_{i-k,j}},0\right)}^{\text{T}}\hfill & \text{ if }\hfill & i>1\hfill & \text{ and }\hfill & j=1\hfill \\ {\left(0,0\right)}^{\text{T}}\hfill & \text{ if }\hfill & i=1\hfill & \text{ and }\hfill & j=1\hfill \end{array}$$

Specially, it is the backward difference form of $\nabla u_{i,j}$ when $\alpha =\text{1}$ . We reorder the image matrices u and u₀ row-wisely into the vector $y$ and $f$ $(y,f\in R^N(N=n\times n))$ , and the $(i,j)$ element of the two-dimension matrices corresponds to the $(j-1)n+i$ element of the vectors, namely $u_{i,j}=y_{(j-1)n+i}$ , ${(u_{\text{0}})}_{i,j}=f_{(j-1)n+i}$ . The $(i,j)$ component of $\nabla u$ and ${\nabla }^{\alpha }u$ can be represented as a multiplication of vector $y$ and matrices $A_s$ , $B_s$ .

For $s=1,2,\cdots ,N$ , we have

$$\nabla u_{i,j}=A_{s}y=\{\begin{array}{lllll}{\left(y_s-y_{s-1},y_s-y_{s-n}\right)}^{\text{T}}\hfill & \text{ if }\hfill & i>1\hfill & \text{ and }\hfill & j>1\hfill \\ {\left(0,y_s-y_{s-n}\right)}^{\text{T}}\hfill & \text{ if }\hfill & i=1\hfill & \text{ and }\hfill & j>1\hfill \\ {\left(y_s-y_{s-1},0\right)}^{\text{T}}\hfill & \text{ if }\hfill & i>1\hfill & \text{ and }\hfill & j=1\hfill \\ {\left(0,0\right)}^{\text{T}}\hfill & \text{ if }\hfill & i=1\hfill & \text{ and }\hfill & j=1\hfill \end{array}$$

and

$${\nabla }^{\alpha }{u}_{i,j}=B_{s}y=\{\begin{array}{lllll}{\left({\displaystyle \sum _{k=0}^K{(-1)}^k{C}_k^{\alpha }{y}_{s-k}},{\displaystyle \sum _{k=0}^K{(-1)}^k{C}_k^{\alpha }{y}_{s-\mathit{nk}}}\right)}^{\text{T}}\hfill & \text{ if }\hfill & i>1\hfill & \text{ and }\hfill & j>1\hfill \\ {\left(0,{\displaystyle \sum _{k=0}^K{(-1)}^k{C}_k^{\alpha }{y}_{s-\mathit{nk}}}\right)}^{\text{T}}\hfill & \text{ if }\hfill & i=1\hfill & \text{ and }\hfill & j>1\hfill \\ {\left({\displaystyle \sum _{k=0}^K{(-1)}^k{C}_k^{\alpha }{y}_{s-k}},0\right)}^{\text{T}}\hfill & \text{ if }\hfill & i>1\hfill & \text{ and }\hfill & j=1\hfill \\ {\left(0,0\right)}^{\text{T}}\hfill & \text{ if }\hfill & i=1\hfill & \text{ and }\hfill & j=1\hfill \end{array}$$

where $s=(j-1)n+i$ , $\text{T}$ is the transpose symbol of matrix. Then equation (14) becomes

$$\underset{y}{\min }{\displaystyle \sum _{s=1}^N(1-g_s)\left|\right|A_{s}y\left|\right|}+{\displaystyle \sum _{s=1}^N{g}_s\left|\right|B_{s}y\left|\right|}+\frac{\lambda }{2}\left|\right|y-f|{|}^2$$

(15)where $A_s$ , $B_s\in R^{\text{2}\times N}$ . In order to simplify equation (15), let ${\overline{A}}_s=(1-g_s)A_s$ , ${\overline{B}}_s=g_s{B}_s$ , then equation (15) finally becomes

$$\underset{y}{\min }{\displaystyle \sum _{s=1}^N\left|\right|{\overline{A}}_{s}y\left|\right|}+{\displaystyle \sum _{s=1}^N\left|\right|{\overline{B}}_{s}y\left|\right|}+\frac{\lambda }{2}\left|\right|y-f|{|}^2$$

(16)

Equation (16) is the primal problem. Next we derive its primal-dual formulation. We first give a conclusion from the Cauchy–Schwarz inequality: for any vector $a$ , there is a vector $b$ satisfies the equation $\left|\right|a\left|\right|=\underset{\left|\right|b\left|\right|\le 1}{\max }{a}^{\text{T}}b.$

Applying the above conclusion to equation (16), we have

$${\displaystyle \sum _{s=1}^N\left|\right|{\overline{A}}_{s}y\left|\right|}=\underset{x}{\max }{y}^{\text{T}}\overline{A}x$$

(17)where $\overline{A}=\left[{\overline{A}}_1^{\text{T}},{\overline{A}}_2^{\text{T}},\cdots ,{\overline{A}}_N^{\text{T}}\right]\in R^{N\times 2N}$ , $x_s={\left[x_s^1,x_s^2\right]}^{\text{T}}\in R^{2\times 1}$ , $x={\left[x_1,x_2,\cdots ,x_N\right]}^{\text{T}}\in R^{2N}.$ In the same way, there is a matrix $\overline{B}$ satisfying the following formulation

$${\displaystyle \sum _{s=1}^N\left|\right|{\overline{B}}_{s}y\left|\right|}=\underset{z}{\max }{y}^{\text{T}}\overline{B}z$$

(18)where $\overline{B}=\left[{\overline{B}}_1^{\text{T}},{\overline{B}}_2^{\text{T}},\cdots ,{\overline{B}}_N^{\text{T}}\right]\in R^{N\times 2N}$ , $z_s={\left[z_s^1,z_s^2\right]}^{\text{T}}\in R^{2\times 1}$ , $z={\left[z_1,z_2,\cdots ,z_N\right]}^{\text{T}}\in R^{2N}.$ Here we define a set

$$C=\left\{c\right|c={[c_1,c_2,\cdots ,c_N]}^{\text{T}},c_i=\left(\begin{array}{c}{c}_i^1\\ c_i^2\end{array}\right),\left|\right|c_i\left|\right|\le 1(i=1,2,\cdots ,N)\}$$

then $x,z\in C$ .

We substitute equations (17) (18) in equation (16) and get the primal-dual formulation

$$\underset{y}{\min }\underset{x,z\in C}{\max }{y}^{\text{T}}\overline{A}x+y^{\text{T}}\overline{B}z+\frac{\lambda }{2}\left|\right|y-f|{|}^2$$

(19)

Its corresponding dual problem is

$$\underset{x,z\in C}{\max }\frac{\lambda }{2}\left(\left|\right|f|{|}^2-||\frac{1}{\lambda }(\overline{A}x+\overline{B}z)-f|{|}^2\right)$$

(20)

Let $\Phi (y,x,z)=y^{\text{T}}\overline{A}x+y^{\text{T}}\overline{B}z+\frac{\lambda }{2}\left|\right|y-f|{|}^2,$ for the multivariate function $\Phi (y,x,z)$ , suppose the primal variable $y$ is fixed and we take the derivative of the dual variable $x,z$

$${\Phi }_x(y,x,z)={\overline{A}}^{\text{T}}y,{\Phi }_z(y,x,z)={\overline{B}}^{\text{T}}y$$

Similarly, fixing the dual variable $x$ and $z$ , we obtain the derivation of the primal variable $y$

$${\Phi }_y(y,x,z)=\overline{B}z+\overline{A}x+\lambda (y-f)$$

According to the gradient descent method, we update the variable $x,z$ and $y$ in turn. To summarize, the primal-dual hybrid gradient algorithm is as follows:

1. Initialization. Pick $y^0$ and feasible $x^0,z^0\in C$ , choose the step size ${\delta }_k,{\tau }_k$ and ${\theta }_k$ .

2. Iteration. For $k=0,1,2,\cdots$ , compute as follows

$$\begin{array}{l}(a)x^{k+1}=P_C(x^k+{\delta }_k{\overline{A}}^{\text{T}}{y}^k)\\ (b)z^{k+1}=P_C(z^k+{\tau }_k{\overline{B}}^{\text{T}}{y}^k)\\ (c)y^{k+1}=y^k-{\theta }_k\left((\overline{A}{x}^{k+1}+\overline{B}{z}^{k+1})/\lambda +y^k-f\right)\end{array}$$

(21)

3. If $y^{k+1}$ satisfies the stopping criterion, then we terminate the iteration and output; otherwise, go to step 2.

In the iterative formulation equation (21), $P_c(v)$ denotes the orthogonal projection of the point $v$ onto the set $C$ , the projection computation ensures that each iteration $x^{k+\text{1}},z^{k+\text{1}}\in C$ expression is

$$P_C(v)=\frac{v}{\max (1,|\left|v\right||)}$$

Silvia Bonettini and Valeria Ruggiero ²⁵ proposed the primal-implicit scheme algorithm and gave its convergence theorem and proof. The primal-dual hybrid gradient algorithm for our proposed model is a special case of primal-implicit scheme algorithm. Compared to algorithm in Bonettini and Ruggiero, ²⁵ the difference of our paper is that we introduce two dual variables with the primal variable to update. The algorithm satisfies the convergence theorem in Bonettini and Ruggiero, ²⁵ theorem 3, and thus the following theorem holds.

Theorem: Let $F(y)={\displaystyle {\sum }_{s=1}^N\left|\right|{\overline{A}}_{s}y\left|\right|}+{\displaystyle {\sum }_{s=1}^N\left|\right|{\overline{B}}_{s}y\left|\right|}+\frac{\lambda }{2}\left|\right|y-f|{|}^2$ , $Y^{*}$ is the set of the $F(y)$ minimum points. If the step size sequences $\left\{{\theta }_k\right\}$ , $\left\{{\tau }_k\right\}$ and $\left\{{\delta }_k\right\}$ satisfy the following conditions: (1) $\underset{k\to \infty }{\lim }{\theta }_k=0$ , $\underset{k\to \infty }{\lim }{\tau }_k=\infty$ and $\underset{k\to \infty }{\lim }{\delta }_k=\infty$ . (2) ${\displaystyle {\sum }_{k=0}^{\infty }\frac{{\theta }_k}{{\tau }_k}<\infty }$ , ${\displaystyle {\sum }_{k=0}^{\infty }\frac{{\theta }_k}{{\delta }_k}<\infty }$ and $\underset{k\to \infty }{\lim }{\theta }_{{}_k}^2<\infty$ . Then $\{F(y^k)\}$ converges to a minimum of $F(y)$ over $R^N$ and $\mathit{dist}(y^k,Y^{*})$ converges to zero.

Numerical experiments

In this section, three experiments verify that our proposed model make better denoising effect compared to the total variation model and the fractional-order variation model. The gray images used for experiments are “Lena” of size 128 × 128, “Women” of size 256 × 256 and “Ellipsoid” of size 256 × 256. The numerical calculations are all accomplished by using MATLAB R2012b.lnk on a laptop with Intel in Core(TM) 22.00 GHZ, 2.00 GB RAM, windows XP.

In these experiments, we take the peak signal to noise ratio (PSNR) as an objective evaluation index of restoration image quality, and its definition is

$$\mathit{PSNR}=10{\log }_{10}\frac{{255}^2}{\mathit{MSE}}\mathrm{dB}$$

$$\mathit{MSE}=\frac{1}{M\times N}{\displaystyle \sum _{i=1}^M{\displaystyle \sum _{j=1}^N{[\overline{u}(i,j)-u(i,j)]}^2}}$$

where $\mathit{MSE}$ is the mean variance measure, $M\times N$ is image size, $\overline{u}$ and u are the original image without noise and the restoration image, respectively. In the implementation of experiments, we use the primal-dual hybrid gradient algorithm for numerical calculation, and the iterative terminating condition of the algorithm is $\frac{\text{||}{y}^{k+\text{1}}-y^k|{|}^2}{\left|\right|y^{k+1}|{|}^2}\le 0.0001$ . In addition, in the experiments, we find that $\lambda =\text{0}\text{.001}$ is the most suitable regularization parameter value.

Improving PSNR

Experiment 1

In this experiment, respectively, we use the total variation model, the fractional variation model and our proposed model to restore the noisy image and compare their PSNR. We add white Gaussian noise with stand variance $\sigma =\text{20}$ to “Lena” image, and the PSNR of noisy image is 22.1412. Figure 1 shows the original image, noisy image and the restoration images from three denoising model. In Table 1, we list the PSNR of restoration image from different models. In order to get a more intuitive result, we make a line chart in Figure 2.

Figure 1.

2D Image: (a) original image; (b) noisy image; (c) restoration image by ROF model; (d) restoration image by FV model with $\alpha =\text{1}\text{.6}$ ; (e) restoration image by our proposed model.

Figure 2.

Line graph of experiment data.

Table 1.

PSNR of restoration images.

$\alpha$	1ROF	1.2	1.4	1.6	1.8
FV model	26.4995	26.6154	26.9291	27.5018	27.1866
The proposed model	26.4995	27.7277	27.5739	27.7252	28.0134

PSNR: peak signal to noise ratio.

In Figure 2, the yellow dot presents the PSNR of the ROF model, the blue dot data come from FV model and the red dot data derive from our proposed model. The FV model line chart general trend rises at first and then falls, our proposed model falls first and then rises, but the PSNR of our proposed model is higher than the ROF model and the FV model, which prove that our proposed model has a better image restoration effect than ROF model and FV model.

Suppressing staircase effect and protecting edges

Experiment 2

In this experiment, we compare the ability of three image restoration models to suppress the staircase effect and protect the edges. In Figure 3, we show four locally enlarged images of “Lena,” including the original image and the restoration images. Comparing image smooth area, we can discover obvious staircase effect in (3b), but few in (3c) and (3d). The experimental results intuitively show that our proposed model can effectively suppress the staircase effect.

Figure 3.

2D Image: “Lena” local original image and local restoration image by different denoised model. (a) Local original image; (b) local restoration image by ROF model; (c) local restoration image by FV model with $\alpha =\text{1}\text{.6}$ ; (d) local restoration image by our proposed model.

In this paper, we use “canny” edge detector to detect image edges. In Figure 5, we show four image edges, including the original image edge and the restoration image edge. Comparing the original image edge (4a) and restoration image edge (4b), we discover that the ROF model has good edge protection capability, but it contains many false artificial edges in hat brim and shoulder area. Comparing the original image edge (4a) and restoration image edge (4c), we can see that (4c) lost a lot of edge detail in hat brim and hair area. It is found that (4d) is the closest to (4a) by comparing four images in Figure 4, which indicates that our proposed model has better edge protection ability than ROF model and FV mode.

Figure 4.

2D Image: “Lena” original image edge and restoration image edge by different denoised model. (a) original image edge; (b) restoration image edge by ROF model; (c) restoration image edge by FV model with $\alpha =\text{1}\text{.6}$ ; (d) restoration image edge by our proposed model.

Preserving texture details

Experiment 3

In this experiment, we add white Gaussian noise with stand variance $\sigma =\text{20}$ to “Women” and “Ellipsoid” image and compare the ability of three restoration models to preserve texture details. Figure 5 displays the original ‘Women’, ‘Ellipsoid’ image and corresponding noisy image. Figures 6 and 7 show the restoration images and corresponding residual images. In Figure 6, comparing (6d), (6e) and (6f), we can clearly discover that (6e) and (6f) contain less texture details than (6d), thus the restoration image (6b) and (6c) retains more texture details than (6a). The experimental results show that the FV model and our proposed model preserve more texture details than ROF model. Analyzing images in Figure 7, we obtain the same conclusion that our proposed model doing well in preserving texture details.

Figure 5.

2D Image: (a) “Women” original image; (b) “Women” noisy image. (c) “Ellipsoid” Original image; (d) “Ellipsoid” noisy image.

Figure 6.

2D Image: “Women” restoration image by different denoised model and corresponding residue image. (a) ROF (b) FV (c) Our proposed model (d) Residual image by ROF (e) Residual image by FV (f) Residual image by our proposed model.

Figure 7.

2–D Image: “Ellipsoid” restoration image by different denoised model and corresponding residue image. (a) ROF (b) FV (c) Our proposed model (d) Residual image by ROF (e) Residual image by FV (f) Residual image by our proposed model.

Scratch repair for damaged images

Experiment 4

In this experiment, we add different scratches noise with stand variance $\sigma =\text{20}$ to “Green peppers” image and compare the ability of three restoration models to preserve texture details. Figure 8 displays the original ‘Green peppers’ image and corresponding scratched image. Figure 8 shows the restoration images and corresponding restored images. In Figure 8, comparing (8d), (8e) and (8f), we can clearly discover that (8e) and (8f) contain less texture details than (6d), thus the restoration image (8d) and (8e) retains more texture details than (8c). The experimental results show that the FV model and our proposed model preserve more texture details than ROF model. Analyzing images in Figure 8, we obtain the same conclusion that our proposed model doing well in preserving texture details.

Figure 8.

2-D Image. (a) Original image (b) Scratched image (c) Scratched image (d) ROF (e) FV (f) Our proposed model.

Conclusions

In this paper, we propose a new model that adaptively combines the total variation model and the fractional-order variation model. We apply the primal-dual hybrid gradient algorithm to solve its energy functional minimization problem. The numerical results show that our proposed model effectively restrains the staircase effect and preserve the texture details. In addition, the model is better than the total variation model and the fractional-order variation model in improving the PSNR and protecting edges.