Mixed Fractional-Order and High-Order Adaptive Image Denoising Algorithm Based on Weight Selection Function

Abstract

In this paper, a mixed-order image denoising algorithm containing fractional-order and high-order regularization terms is proposed, which effectively suppresses the staircase effect generated by the TV model and its variants while better preserving the edges and details of the image. Adding different regularization penalties in different regions is fundamental to improving the denoising performance of the model. Therefore, a weight selection function is designed using the structure tensor to achieve a more effective selection of regularization terms in different regions. In each iteration, the regularization parameters are adaptively adjusted according to the Morozov discrepancy principle to promote the performance of the algorithm. Based on the primal–dual theory, the original algorithm is improved by using the predictor–corrector scheme to obtain a more accurate approximate solution while ensuring the convergence of the algorithm. The effectiveness of the proposed algorithm is demonstrated through simulation experiments.

1. Introduction

The field of digital image processing requires high-resolution and clean images. However, in most cases, images are often degraded due to noise and other factors, resulting in the loss of valuable image information. Therefore, it poses a significant challenge to model the noise accurately, employ efficient algorithms to remove it from the image, and achieve algorithm convergence.

Image denoising can be regarded as a reversible problem, but it is ill-posed, which means that the solution to this problem does not exist or is not unique. To overcome this problem and achieve better denoising performance, several methods have been proposed, including regularization methods [1], wavelet transforms [2,3], non-local methods such as NLM [4,5] and BM3D [6,7], deep learning methods such as DnCNN [8], DIVA [9], and BM3D-net [10], and quantum-based methods such as DeQuIP [11] and QAB [12]. The role of these methods in image denoising is acknowledged, but certain limitations are also present in them. For instance, the appearance of new unrealistic structures or textures is prone to occur in denoised images using the wavelet transform. Non-local methods are characterized by high computational complexity and are easily influenced by outliers. Deep learning methods necessitate a substantial dataset for model training, leading to longer modeling times. Quantum-based methods are still in the early stages of research and have many technical challenges and limitations. Although there are some limitations to regularization methods, they are supported by a solid understanding of mathematical theory, and the logical reasoning of the computational derivation process is also strong. Being easy to code and having a shorter modeling time, they are still widely utilized at present.

In models associated with regularization techniques, there are two prevalent issues: how to suppress the staircase effect caused by the total variation (TV) regularization model proposed by Rudin et al. [13], and how to adaptively adjust the regularization parameters based on image features. Currently, a considerable amount of research has been conducted on these problems. Regarding the first issue, researchers have presented using higher-order regularizers [14,15,16] or fractional-order regularizers [17,18,19] to modify the TV regularization term and suppress the staircase effect. However, studies have found that the higher-order TV model (HOTV) does not smooth noise as effectively as the TV model. Although the fractional-order TV model (FOTV) has better noise filtering and detail preservation effects, its solving process is challenging. For the second issue, researchers have designed different adaptive regularization parameters based on prior information about the image [20,21]. However, there are few algorithms that simultaneously address both of these problems. Therefore, the advantages of the HOTV and FOTV models are inherited in this paper, and a mixed-order TV (MOTV) image denoising algorithm is proposed. Moreover, adaptive regularization parameters are designed based on image features. To efficiently solve the model while minimizing the number of parameters, an improved primal–dual algorithm based on the predictor–corrector scheme is proposed, and the convergence of this algorithm is theoretically proven.

The main contributions of this paper are summarized below: Firstly, a new MOTV image denoising algorithm is proposed to suppress the staircase artifacts while filtering out more noise. The approximate optimal solution of the model is obtained by using the improved primal–dual algorithm. Secondly, an edge-weighted selection function is constructed by using the structure tensor, which preserves more edge details of the image. Finally, an adaptive regularization parameter is designed using the Morozov discrepancy principle to facilitate the performance of the algorithm.

The rest of this article is organized as follows: Section 2 briefly introduces the relevant knowledge background and model. In Section 3, a new model is proposed, and the derivation process of the model solution is given. Then an adaptive regularization parameter is set, and the convergence of the algorithm is evidenced. In Section 4, simulation experiments are carried out to prove the effectiveness of the proposed algorithm. Section 5 provides a simple summary.

2. Related Works

2.1. Fractional-Order Variation Model

The fractional-order derivative, as a generalization of the integer derivative, can preserve more image details. Therefore, it has been extensively studied by a quantity of scholars in recent decades [22,23]. Firstly, the definition of fractional-order derivatives and their related properties are reviewed here.

The Euclidean space ${\mathbb{R}}^{m\times n}$ of an image u with size $m\times n$ is represented as V. To facilitate interpretation, the image domain is discretized into a matrix grid $\{(x_i,y_i):1\le i\le m,1\le j\le n\}$, where $(i,j)$ represents the position of the pixel in the matrix. Consequently, the image $u\subset V$ can be defined as $u_{i,j}=u(x_i,y_j)$. Various definitions of fractional-order derivatives exist in the literature [24], including the Grünwald–Letnikov (G-L) definition, the Riemann–Liouville (R-L) definition, the Caputo definition, and others. In digital image processing, the G-L definition is the most popular. The G-L definition is also employed in this paper, and the relevant fractional-order gradient is defined as follows:

$${\nabla }^{\alpha }u={({\nabla }_1^{\alpha }u,{\nabla }_2^{\alpha }u)}^T,\alpha >0,$$

(1)

where ${\nabla }_1^{\alpha }u,{\nabla }_2^{\alpha }u\in V$ denote the derivatives of the image u in the horizontal and vertical directions, respectively, which can be discretized and approximated as

$${({\nabla }_1^{\alpha }u)}_{i,j}=\sum _{k=0}^{K-1}{(-1)}^k{C}_k^{\alpha }{u}_{i-k,j},{({\nabla }_2^{\alpha }u)}_{i,j}=\sum _{k=0}^{K-1}{(-1)}^k{C}_k^{\alpha }{u}_{i,j-k}.$$

(2)

Here $K\ge 3$ is the number of adjacent pixels used to approximate the fractional-order derivative at each pixel point. The coefficient ${\left\{C_k^{\alpha }\right\}}_{k=0}^{K-1}$ is defined as $C_k^{\alpha }=\frac{\Gamma (\alpha +1)}{\Gamma (k+1)\Gamma (\alpha +1-k)}$, where $\Gamma (·)$ is the Gamma function, $U=\left\{u\in V:{∥u∥}_{\infty }={max}_{ij}\left|u_{i,j}\right|\le 1\right\}$. Then the FOTV model can be defined as

$$\underset{u}{min}\frac{\lambda }{2}{\parallel Hu-f\parallel }_2^2+{\parallel {\nabla }^{\alpha }u\parallel }_1,$$

(3)

where

$$\left|\right|{\nabla }^{\alpha }u{\left|\right|}_1:=\sum _{ij}|{({\nabla }^{\alpha }u)}_{i,j}|,\left|{({\nabla }^{\alpha }u)}_{i,j}\right|=\sqrt{{({\nabla }_1^{\alpha }u)}_{i,j}^2+{({\nabla }_2^{\alpha }u)}_{i,j}^2},$$

(4)

Here, H is a blur operator, u is the original image, and f is the noisy image. According to Ref. [25], an important relation can be derived as $\langle {\nabla }^{\alpha }u,p\rangle =\langle u,\overline{{(-1)}^{\alpha }}{\mathrm{div}}^{\alpha }p\rangle$. Here, $\forall p=(p^1,p^2)\in Z$ where $Z=V\times V$, and p represents a dual variable. Utilizing the aforementioned relationship, the adjoint divergence of the discrete gradient (1) can be obtained as follows:

$${({\mathrm{div}}^{\alpha }p)}_{i,j}={(-1)}^{\alpha }\sum _{k=0}^{K-1}{(-1)}^k{C}_k^{\alpha }(p_{i+k,j}^1+p_{i,j+k}^2).$$

(5)

2.2. High-Order Variation Model

If the TV model is defined in the bounded variation (BV) function space, then the reconstruction of discontinuous points can be achieved. Accordingly, the HOTV models, which can better model noise, have been proposed by plenty of scholars to suppress the staircase effect and achieve improved results [14,15,16,26]. The HOTV model is generally defined as follows:

$$\underset{u}{min}{∥Hu-f∥}_2^2+\beta R(u),$$

(6)

where $\beta$ is a predefined parameter and $R(u)$ is a high-order regularization term.

2.3. Structural Tensor

In mathematics, the structure tensor, also known as the second moment matrix, is a matrix derived from the gradient of a function. It can summarize the main directions of the gradients in a specified neighborhood of a given point, as well as the degree of coherence in those directions. As a result, the structure tensor finds extensive application in the fields of computer vision and image processing [27,28]. In image processing, it excels at distinguishing between edge regions, flat regions, and corner regions. The structure tensor for the image u can be expressed as

$$ST=\left[\begin{array}{cc}{u}_x^2& u_x{u}_y\\ u_y{u}_x& u_y^2\end{array}\right].$$

(7)

Here, $u_x$, $u_y$ denote the gradients of the image u in the x and y directions, respectively. If $SK$ and $SH$ are used to represent the determinant and trace of matrix $ST$, respectively, then the judgment conditions for each region are summarized as follows:

• flat regions: $SH=0$;
• edge regions: $SH>0$ and $SK=0$;
• corner regions: $SH>0$ and $SK>0.$

In Figure 1, the white portions of subgraphs (a–c), respectively, represent the flat regions, edge regions, and corner regions of the Lena image. However, in practice, the values of $SK$ and $SH$ are not ideal, so approximate judgment is usually used.

2.4. Primal–Dual Algorithm with Guaranteed Convergence

The standard Euclidean inner product and norm of the finite dimensional vector spaces X and Y are expressed as $\langle ․,․\rangle$ and $\left|\right|․\left|\right|={\langle ․,․\rangle }^{\frac{1}{2}}$, respectively. If a continuous linear operator A is a map on $X\to Y$, then its induced norm can be expressed as

$$L=\left|\right|A\left|\right|=max\left\{\right|\left|Ax\right||:x\in X,|\left|x\right||\le 1\}.$$

(8)

According to Fenchel duality theory [29], the generalized saddle point problem can be expressed as

$$\underset{x\in X}{min}\underset{y\in Y}{max}\langle Ax,y\rangle +G(x)-F^{*}(y),$$

(9)

where G and F are lower semi-continuous functions, and $F^{*}$ is the conjugate function of F.

It is well known that the primal–dual problem of the following primal optimization problem (10) corresponds to this saddle point problem (9) mentioned above. The primal problem can be expressed as follows:

$$\underset{x\in X}{min}F(Ax)+G(x).$$

(10)

Then, according to the primal–dual problem, the dual problem related to the primal problem can be calculated, which is expressed as

$$\underset{y\in Y}{max}-(G^{*}(-A^{*}y))+F^{*}(y).$$

(11)

Chambler and Pock proposed a primal–dual algorithm for solving the saddle point structure problem (9), and the convergence of the algorithm was proven [30]. By fixing the primal variable x of the primal–dual problem and solving the first-order derivative function about the dual variable y, the expression for y can be obtained:

$$y={(I+\partial F^{*})}^{-1}(y+Ax).$$

(12)

Similarly, it can be obtained that the expression about the primal variable x is:

$$x={(I+\partial G)}^{-1}(x-A^{*}y)$$

(13)

where $\partial F^{*}$ and $\partial G$ represent the subgradients of $F^{*}$ and G, respectively. The operators ${(I+\partial F^{*})}^{-1}$ and ${(I+\partial G)}^{-1}$ are known as resolvent operators. When the concept of minimizing an energy functional is applied to define the resolvent operator, it can be expressed as

$$x={(I+\tau \partial F)}^{-1}(y)=arg\underset{x}{min}\left\{\frac{{∥x-y∥}^2}{2\tau }+F(x)\right\}.$$

See [31] for details of the subgradient and resolvent operators. Then, the primal–dual algorithm proposed in the literature [30] can ultimately be summarized as

• Initialization: Choose $t_0,s_0>0$ with $t_0{s}_0{L}^2\le 1,(x^0,y^0)\in X\times Y$, and ${\overline{x}}^0=x^0$.
• Iteration $(n\ge 0)$: update as follows

  $$\left\{\begin{array}{cc}& y^{n+1}={(I+s_n\partial F^{*})}^{-1}(y^n+s_{n}A{\overline{x}}^n),\hfill \\ & x^{n+1}={(I+t_n\partial G)}^{-1}(x^n-t_n{A}^{*}{y}^{n+1}),\hfill \\ & {\theta }_n=1/\sqrt{1+2\gamma t_n},t_{n+1}={\theta }_n{t}_n,s_{n+1}=s_n/{\theta }_n,\hfill \\ & {\overline{x}}^{n+1}=x^{n+1}+{\theta }_n(x^{n+1}-x^n).\hfill \end{array}\right.$$

3. Model Proposal and Solution

3.1. Proposed Model

In this study, the following image denoising problems with additive noise are considered:

$$f=Hu+\eta ,$$

(14)

where $\eta$ represents the Gaussian noise. TV regularization is a classical method for restoring high-quality images, guaranteeing the stability of the solution. Models (3) and (6) are classical FOTV and HOTV models, respectively, used for removing Gaussian noise. Building upon the strengths of both models, a mixed-order image denoising model is proposed in this paper, described as follows:

$$\underset{u}{min}\left\{\frac{\lambda }{2}{\parallel Hu-f\parallel }_2^2+{\xi }_i{\parallel {\nabla }^{\alpha }u\parallel }_1+{\xi }_i{\parallel {\nabla }^{2}u\parallel }_2+\beta {\parallel u\parallel }_1\right\}.$$

(15)

Here, $\lambda$ represents the regularization parameter, $\beta$ is the balance parameter, and ${\xi }_i(i=1,2)$ is a weight selection parameter. The purpose of introducing this weight selection parameter is to select different regions of the image by the structure tensor and to process the image using different regularization terms in different regions. To achieve this purpose, ${\xi }_i$ is defined as a function of pixel point characteristics. In other words, the intensity of regularization varies for pixel points in the flat regions, edge regions, and corner regions. At present, the weight function, which is the inverse function of the image gradient, is used by most researchers [29,32,33], but this paper takes a different approach.

Next, a matrix related to the structure tensor is defined:

$$K_{\sigma }(u(x,y))=G_{\sigma }*(\nabla u\nabla u^T)=\left[\begin{array}{cc}{G}_{\sigma }*u_x^2& G_{\sigma }*u_x{u}_y\\ G_{\sigma }*u_y{u}_x& G_{\sigma }*u_y^2\end{array}\right]=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right],$$

where $G_{\sigma }={(\sigma \sqrt{2\pi })}^{-1}exp(-{\left|x\right|}^2/2{\sigma }^2)$ is a Gaussian filter function with standard deviation $\sigma$.

The matrix K contains two eigenvectors that are standard orthogonal. Different eigenvalues ${\lambda }_1,{\lambda }_2$ describe the average contrast in the corresponding eigenvector direction. Different regions exhibit different eigenvalue relationships. In the edge regions, ${\lambda }_1\gg {\lambda }_2\approx 0$; in the flat regions, ${\lambda }_1,{\lambda }_2\approx 0$; and in the corner regions, ${\lambda }_1>{\lambda }_2\gg 0$ [34]. Hence, the characteristics of each pixel can be measured using the difference between the eigenvalues. The edge coherence [35] based on the eigenvalue is defined as

$$C(u)={({\lambda }_1-{\lambda }_2)}^2={(a_{11}-a_{22})}^2+4a_{12}^2.$$

(16)

Furthermore, the weight selection function ${\xi }_i$ can be defined as

$${\xi }_i=\left\{\begin{array}{cc}1-exp[-\frac{1}{{(C(u)/\varpi )}^2}],i\hfill & =1\\ exp[-\frac{1}{{(C(u)/\varpi )}^2}],i\hfill & =2\end{array}\right.$$

(17)

where $\varpi \in [{10}^{-4}\sim {10}^{-2}]$ is the comparison parameter. The role of $\varpi$ is equivalent to a threshold, which can be used to distinguish various regions. Note that in order to prevent the impact of abnormal values, $C(u)$ should be normalized before calculating ${\xi }_i$.

Obviously, ${\xi }_i$ depends on the edge consistency of pixel features. In the edge and corner regions, $C\to 1$, ${\xi }_1\to 0$, ${\xi }_2\to 1$, at this time, ${\xi }_2{∥{\nabla }^{2}u∥}_2+\beta {∥u∥}_1$ is selected as a regularization term of the model, which can better preserve the edge details of the image. In the flat regions, $C\to 0$, ${\xi }_2\to 0$, ${\xi }_1\to 1$, at this time, ${\xi }_1{∥{\nabla }^{\alpha }u∥}_1+\beta {∥u∥}_1$ is selected as a regularization term, which can reduce noise and artifacts in the image while also reinforcing the regularization process for the flat regions. By designing the selection function ${\xi }_i$ and performing correspondence matching within the model, the image edges and fine details can be effectively preserved, ultimately resulting in the recovery of high-quality images. Below is the corresponding matching model:

$$\underset{u}{min}\Phi (u,\lambda )\equiv \left\{\frac{\lambda }{2}{\parallel Hu-f\parallel }_2^2+{\xi }_1{\parallel {\nabla }^{\alpha }u\parallel }_1+{\xi }_2{\parallel {\nabla }^{2}u\parallel }_2+\beta {\parallel u\parallel }_1\right\}.$$

(18)

Solving the MOTV model is equivalent to solving the saddle point of the primal–dual problem. The iterative steps and convergence analysis of the model are introduced in detail below.

3.2. Numerical Algorithm

The MOTV model is a variant of the TV model. The famous Euler–Lagrange method can be utilized to solve both the TV model and its variant. However, the model contains a non-differentiable derivative, making its calculation process complex. To overcome the difficulties, the primal–dual algorithm is employed. The conjugate sets of all the dual variables $p,y$, and w are $P,Y$, and W, respectively. Let $\delta (u)={\parallel {\nabla }^{\alpha }u\parallel }_1$, then its conjugate function satisfies the following properties:

$${\delta }_P^{*}(p)=\left\{\begin{array}{cc}0,\hfill & p\in P\hfill \\ +\infty ,\hfill & p\notin P\hfill \end{array}\right.$$

(19)

where $P=\left\{p\in Z:{\parallel p\parallel }_{\infty }={max}_{ij}\left|p_{i,j}\right|\le 1\right\}$ denotes the conjugate set. Similarly, let $\zeta (u)={\parallel {\nabla }^{2}u\parallel }_2$, $\varsigma (u)={\parallel u\parallel }_1$, then the associated conjugate functions ${\zeta }_Y^{*}(y)$, ${\varsigma }_W^{*}(w)$ and conjugate sets Y, W also have the same properties.

Applying Fenchel duality theory to the primal problem (18), a general primal–dual model of the MOTV model can be obtained, which is expressed as

$$\begin{array}{cc}\hfill \underset{u}{min}\underset{p,y,w}{max}\frac{\lambda }{2}{\parallel Hu-f\parallel }_2^2& +{\xi }_1\langle p,{\nabla }^{\alpha }u\rangle +{\xi }_2\langle y,{\nabla }^{2}u\rangle +\beta \langle w,u\rangle \hfill \\ \hfill & -{\xi }_1{\delta }_P^{*}(p)-{\xi }_2{\zeta }_Y^{*}(y)-\beta {\varsigma }_W^{*}(w).\hfill \end{array}$$

(20)

Of these, $y,w$ has the same properties as p. Similarly, the corresponding dual problem can be obtained as follows:

$$\begin{array}{cc}\hfill \underset{p,y,w}{max}\langle {\xi }_1{({\nabla }^{\alpha })}^{T}p+{\xi }_2{({\nabla }^2)}^{T}y& +\beta w,f\rangle -\frac{\lambda }{2}{\parallel {\xi }_1{({\nabla }^{\alpha })}^{T}p+{\xi }_2{({\nabla }^2)}^{T}y+\beta w\parallel }_2^2\hfill \\ \hfill & -{\xi }_1{\delta }_P^{*}(p)-{\xi }_2{\zeta }_Y^{*}(y)-\beta {\varsigma }_W^{*}(w).\hfill \end{array}$$

(21)

The discretized $MTV$ of the image u is also the Legendre–Fenchel conjugate of $\delta$, $\zeta$, and $\varsigma$ [36,37]. That is:

$$\begin{array}{cc}\hfill MTV(u)& =\underset{p,y,w}{max}\{{\xi }_1\langle p,{\nabla }^{\alpha }u\rangle +{\xi }_2\langle y,{\nabla }^{2}u\rangle +\beta \langle w,u\rangle \hfill \\ & -{\xi }_1{\delta }_P^{*}(p)-{\xi }_2{\zeta }_Y^{*}(y)-\beta {\varsigma }_W^{*}(w)\}\hfill \\ & =\underset{p,y,w}{max}\{{\xi }_1\langle p,{\nabla }^{\alpha }u\rangle +{\xi }_2\langle y,{\nabla }^{2}u\rangle +\beta \langle w,u\rangle \}.\hfill \end{array}$$

Legendre–Fenchel’s duality is used to represent the $MTV(u)$ norm, (18) can be represented as

$$\underset{u}{min}\underset{p,y,w}{max}L(u,p,y,w;\lambda ).$$

(22)

where $L(u,p,y,w;\lambda )=\frac{\lambda }{2}{\parallel Hu-f\parallel }_2^2+{\xi }_1\langle p,{\nabla }^{\alpha }u\rangle +{\xi }_2\langle y,{\nabla }^{2}u\rangle +\beta \langle w,u\rangle$. Let $G(u)=\frac{\lambda }{2}{\parallel Hu-f\parallel }_2^2$, then the following equality holds [38].

$$p^{k+1}={(s{\xi }_1\partial {\delta }^{*}+I)}^{-1}\tilde{p},\tilde{p}=p^k+s{\xi }_1{\nabla }^{\alpha }{u}^k$$

(23)

$$y^{k+1}={(s{\xi }_2\partial {\zeta }^{*}+I)}^{-1}\tilde{y},\tilde{y}=y^k+s{\xi }_2{\nabla }^2{u}^k$$

(24)

$$w^{k+1}={(s\beta \partial {\varsigma }^{*}+I)}^{-1}\tilde{w},\tilde{w}=w^k+s\beta u^k$$

(25)

$$u^{k+1}={(t\partial G+I)}^{-1}\tilde{u},\tilde{u}=u^k-t{\xi }_1{({\nabla }^{\alpha })}^T{p}^{k+1}-t{\xi }_2{({\nabla }^2)}^T{y}^{k+1}-t\beta w^{k+1}$$

(26)

Among them, four resolvent operators are involved, namely ${(\partial {\delta }^{*}+I)}^{-1}$, ${(\partial {\zeta }^{*}+I)}^{-1}$, ${(\partial {\varsigma }^{*}+I)}^{-1}$, ${(\partial G+I)}^{-1}$, the solution of the resolvent operators is very difficult. Therefore, the method of solving saddle points can be employed to equivalently solve the above equations and reduce computational complexity. According to the literature [39], the function $L(u,p,y,w;\lambda )$ has a saddle point, which can be obtained by Lemma 1 as $(u^{*},p^{*},y^{*},w^{*})$.

Lemma 1.

If $(u^{*},p^{*},y^{*},w^{*})$ is a point of the function $L(u,p,y,w;\lambda )$, and the point satisfies the inequality

$$L(u^{*},p,y,w;\lambda )\le L(u^{*},p^{*},y^{*},w^{*};\lambda )\le L(u,p^{*},y^{*},w^{*};\lambda ),$$

then, this point can be represented as a saddle point of the function $L(u,p,y,w;\lambda )$.

In connection with the literature [40], it is clear that the solution procedure for the maximum and minimum of (22) is interchangeable, i.e.,

$$\underset{u}{min}\underset{p,y,w}{max}L(u,p,y,w;\lambda )=L(u^{*},p^{*},y^{*},w^{*};\lambda )=\underset{p,y,w}{max}\underset{u}{min}L(u,p,y,w;\lambda ).$$

This further reflects that the minimalized primal problem can be solved by finding the saddle point of the primal–dual problem, where the dual variable is updated alternately with the primal variable [41,42]. If $u^{*}$, $p^{*}$, $y^{*}$, $w^{*}$ is the optimal solution of (22), then $(u^{*},p^{*},y^{*},w^{*})$ is a saddle point of (22) [39], i.e.,

$$u^{*}=\underset{u}{argmin}L(u,p^{*},y^{*},w^{*};\lambda ),$$

(27)

$$p^{*}=\underset{p}{argmin}L(u^{*},p,y^{*},w^{*};\lambda ),$$

(28)

$$y^{*}=\underset{y}{argmin}L(u^{*},p^{*},y,w^{*};\lambda ),$$

(29)

$$w^{*}=\underset{w}{argmin}L(u^{*},p^{*},y^{*},w;\lambda ).$$

(30)

Combining (22) and (27) yields

$$\lambda H^T(H{u}^{*}-f)+{\xi }_1{({\nabla }^{\alpha })}^T{p}^{*}+{\xi }_2{({\nabla }^2)}^T{y}^{*}+\beta w^{*}=0.$$

(31)

According to the Karush–Kuhn–Tucker (KKT) necessary conditions of the dual optimality, the Lagrange multipliers ${{\mu }^1}_{i,j},{{\mu }^2}_{i,j},{{\mu }^3}_{i,j}\ge 0$ can be yielded, such that

$$\left\{\begin{array}{cc}& {({\nabla }^{\alpha }{u}^{*})}_{i,j}+{{\mu }^1}_{i,j}{p}_{i,j}^{*}=0,\hfill \\ & {({\nabla }^2{u}^{*})}_{i,j}+{{\mu }^2}_{i,j}{y}_{i,j}^{*}=0,\hfill \\ & u_{i,j}^{*}+{{\mu }^3}_{i,j}{w}_{i,j}^{*}=0.\hfill \end{array}\right.$$

where either ${{\mu }^1}_{i,j}>0$ with $|p_{i,j}|=1$ or ${{\mu }^1}_{i,j}=0$ with $|p_{i,j}|<1$. At this time, ${{\mu }^1}_{i,j}=\left|{({\nabla }^{\alpha }{u}^{*})}_{i,j}\right|$. Similarly, the other two formulas can be obtained. Thus, if the solutions of (28)–(30) are $p^{*}$, $y^{*}$, $w^{*}$, respectively, then

$$\left\{\begin{array}{c}{({\nabla }^{\alpha }{u}^{*})}_{i,j}+\left|{({\nabla }^{\alpha }{u}^{*})}_{i,j}\right|p_{i,j}^{*}=0,\hfill \\ {({\nabla }^2{u}^{*})}_{i,j}+\left|{({\nabla }^2{u}^{*})}_{i,j}\right|y_{i,j}^{*}=0,\hfill \\ u_{i,j}^{*}+\left|u_{i,j}^{*}\right|w_{i,j}^{*}=0.\hfill \end{array}\right.$$

(32)

Therefore, the following lemma can be obtained:

Lemma 2.

Suppose p, y, w are dual variables, and if a saddle point of the function $L(u,p,y,w;\lambda )$ is $(u^{*},p^{*},y^{*},w^{*})$, then (31) and (32) hold.

There exist plenty of methods available for computing saddle points through maximization and minimization techniques [37,43,44]. However, the hybrid gradient primal–dual algorithm is utilized by this paper to solve saddle points, which improves the original primal–dual algorithm. By employing this method to calculate the primal and dual variables alternately, a solution can be obtained as follows:

$$p^{k+1}=\underset{p}{argmax}\left\{L(u,p,y,w;\lambda )-\frac{1}{2s}{\parallel p-p^k\parallel }_2^2\right\},$$

(33)

$$y^{k+1}=\underset{y}{argmax}\left\{L(u,p,y,w;\lambda )-\frac{1}{2s}{\parallel y-y^k\parallel }_2^2\right\},$$

(34)

$$w^{k+1}=\underset{w}{argmax}\left\{L(u,p,y,w;\lambda )-\frac{1}{2s}{\parallel w-w^k\parallel }_2^2\right\},$$

(35)

$$u^{k+1}=\underset{u}{argmin}\left\{L(u,p,y,w;\lambda )+\frac{1}{2t}{\parallel u-u^k\parallel }_2^2\right\}.$$

(36)

Here, $s,t>0$ represent the update steps for the dual and primal variables, respectively, which satisfies $s·t\le 1/16$. Note that the regularization parameter $\lambda$ is not involved in the solution of the dual variables $p,y,w$. Before discussing the solution of the subproblem (33)–(36), the projection of q on P is first defined as follows:

$${\rho }_P(q)=\underset{p\in P}{argmin}{\parallel p-q\parallel }_2^2.$$

(37)

The Lagrange method can be used to solve (37). The corresponding Lagrangian function is expressed as

$${\parallel p-q\parallel }_2^2+\sum _{i,j}{\gamma }_{i,j}({\parallel p_{i,j}\parallel }^2-1).$$

Under the constraint ${|p_{i,j}|}^2\le 1$, the Lagrange multiplier ${\gamma }_{i,j}\ge 0$ can be obtained. It follows from the complementarity condition that either ${\gamma }_{i,j}=0$ with $|p_{i,j}|<1$ and $|q_{i,j}|<1$ or ${\gamma }_{i,j}>0$ with $|p_{i,j}|=1$ and $|q_{i,j}|\ge 1$. In the first case, $p_{i,j}=q_{i,j}$. In the other case, the KKT condition yields

$$p_{i,j}-q_{i,j}+{\gamma }_{i,j}{p}_{i,j}=0,\forall i,j.$$

Furthermore, ${\gamma }_{i,j}=\left|q_{i,j}\right|-1$. Thus, $p_{i,j}=q_{i,j}/|q_{i,j}|$. Then, the expression of ${\rho }_P(q)$ is:

$$({\rho }_P(q))=\frac{q}{max(1,|q|)}.$$

(38)

In what follows, the solution to each subproblem is discussed in detail. For the p-subproblem in (33), the solution of

$$\begin{array}{cc}\hfill p^{k+1}& =arg\underset{p}{max}{\xi }_1\langle p-p^k,{\nabla }^{\alpha }u\rangle -\frac{1}{2s}{\parallel p-p^k\parallel }_2^2\hfill \\ \hfill & =arg\underset{p}{max}\frac{1}{2s}\{2\langle p-p^k,s{\xi }_1{\nabla }^{\alpha }u\rangle -\parallel p-p^k{\parallel }_2^2\}\hfill \\ \hfill & =arg\underset{p}{min}{\parallel p-(p^k+s{\xi }_1{\nabla }^{\alpha }u)\parallel }_2^2.\hfill \end{array}$$

(39)

The minimization problem is equivalent to computing the projection of $(p^k+s{\xi }_1{\nabla }^{\alpha }u)$ onto P set. Therefore

$$\begin{array}{c}\hfill p^{k+1}={\rho }_P(p^k+s{\xi }_1{\nabla }^{\alpha }u).\end{array}$$

(40)

Similarly, $y^{k+1}$ in (34) and $w^{k+1}$ in (35) can be computed by

$$\begin{array}{cc}\hfill y^{k+1}& ={\rho }_Y{(}^{k}y+s{\xi }_2{\nabla }^{2}u),\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill w^{k+1}& ={\rho }_W(w^k+s\beta u).\hfill \end{array}$$

(42)

For the u-subproblem in (36), the subproblem (36) is a quadratic function of image u, which has a closed-form solution. The solution process can be expressed as follows:

$$\begin{array}{cc}\hfill u^{k+1}=& arg\underset{u}{min}{\xi }_1\langle p,{\nabla }^{\alpha }(u-u^k)\rangle +\frac{\lambda }{2}{\parallel Hu-f\parallel }_2^2+{\xi }_2\langle y,{\nabla }^2(u-u^k)\rangle \hfill \\ \hfill & +\beta \langle w,u-u^k\rangle +\frac{1}{2t}{\parallel u-u^k\parallel }_2^2,\hfill \end{array}$$

(43)

hence, $u^{k+1}$ can be easily computed by

$$\begin{array}{c}\hfill u^{k+1}={(t\lambda H^{T}H+I)}^{-1}(t\lambda H^{T}f+\tilde{u}).\end{array}$$

(44)

Here, $\tilde{u}=(u^k-t\xi {({\nabla }^{\alpha })}^{T}p-t\xi {({\nabla }^2)}^{T}y-t\beta w)$. Then, the solutions to each subproblem are summarized as follows:

$$\begin{array}{cc}\hfill p^{k+1}& ={(s{\xi }_1\partial {\delta }^{*}+I)}^{-1}\tilde{p}\iff p^{k+1}=\frac{p^k+s{\xi }_1{\nabla }^{\alpha }{u}^k}{max(1,|p^k+s{\xi }_1{\nabla }^{\alpha }{u}^k|)},\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill y^{k+1}& ={(s{\xi }_2\partial {\zeta }^{*}+I)}^{-1}\tilde{y}\iff y^{k+1}=\frac{y^k+s{\xi }_2{\nabla }^2{u}^k}{max(1,|y^k+s{\xi }_2{\nabla }^2{u}^k|)},\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill w^{k+1}& ={(s\beta \partial {\varsigma }^{*}+I)}^{-1}\tilde{w}\iff w^{k+1}=\frac{w^k+s\beta u^k}{max(1,|w^k+s\beta u^k|)},\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill u^{k+1}& ={(t\partial E+I)}^{-1}\tilde{u}\iff u^{k+1}={(t\lambda H^{T}H+I)}^{-1}(t\lambda H^{T}f+\tilde{u}).\hfill \end{array}$$

(48)

To obtain more accurate computational results, the primal–dual algorithm can be further enhanced using the predictive correction method [43]. This method also facilitates the convergence analysis of the algorithm.

In summary, the solution process can be summarized as the following Algorithm 1.

Algorithm 1: Primal–Dual Algorithm for MOTV Image Restoration

Require:  f, H, $\beta$. 1: Initialize u, p, y and w. Set the size of s and t. 2: while stopping criterion is not satisfied do 3:     $C(u)={({\lambda }_1-{\lambda }_2)}^2={(a_{11}-a_{22})}^2+4a_{12}^2$; 4:     ${\xi }_1(u)=1-exp[-\frac{1}{{(C(u)/\gamma )}^2}]$; 5:     ${\xi }_2(u)=exp[-\frac{1}{{(C(u)/\gamma )}^2}]$; 6:     $p^{k+\frac{1}{2}}=\frac{p^k+s{\xi }_1{\nabla }^{\alpha }{u}^k}{max(1,|p^k+s{\xi }_1{\nabla }^{\alpha }{u}^k|)}$; 7:     $y^{k+\frac{1}{2}}=\frac{y^k+s{\xi }_2{\nabla }^2{u}^k}{max(1,|y^k+s{\xi }_2{\nabla }^2{u}^k|)}$; 8:     $w^{k+\frac{1}{2}}=\frac{w^k+s\beta u^k}{max(1,|w^k+s\beta u^k|)}$; 9:     $\tilde{u}=u^k-t{\xi }_1{({\nabla }^{\alpha })}^T{p}^{k+\frac{1}{2}}-t{\xi }_2{({\nabla }^2)}^T{y}^{k+\frac{1}{2}}-t\beta w^{k+\frac{1}{2}}$; 10:     $u^{k+1}={(t\lambda H^{T}H+I)}^{-1}(t\lambda H^{T}f+\tilde{u})$; 11:     $p^{k+1}=\frac{p^{k+\frac{1}{2}}+s{\xi }_1{\nabla }^{\alpha }{u}^{k+1}}{max(1,|p^{k+\frac{1}{2}}+s{\xi }_1{\nabla }^{\alpha }{u}^{k+1}|)}$; 12:     $y^{k+1}=\frac{y^{k+\frac{1}{2}}+s{\xi }_2{\nabla }^2{u}^{k+1}}{max(1,|y^{k+\frac{1}{2}}+s{\xi }_2{\nabla }^2{u}^{k+1}|)}$; 13:     $w^{k+1}=\frac{w^{k+\frac{1}{2}}+s\beta u^{k+1}}{max(1,|w^{k+\frac{1}{2}}+s\beta u^{k+1}|)}$; 14: end while 15: return $\mathbf{u}={\mathbf{u}}^{k+1}$.

3.3. Adaptive Regularization Parameter

Algorithm 1 involves an important regularization parameter $\lambda$ when solving the subproblem u. In fact, $\lambda$ has a tremendous influence on the denoising performance of the algorithm. In general, a multitude of researchers have manually set the regularization parameter $\lambda$, and several researchers have adaptively selected the parameter $\lambda$ according to the prior information of the image [45,46]. In order to simplify the solution process, this paper designs an adaptive parameter $\lambda$ by using the structural features of the data fidelity term while combining the Morozov deviation principle [40], so that the denoised image is always in the feasible set D. The feasible set D is defined as

$$D=\left\{{u:\parallel Hu-f\parallel }_2^2\le c^2\right\}.$$

(49)

Here c is the positive number that depends on the noise variance ${\sigma }^2$[17,40]. If the noise variance ${\sigma }^2$ is known, then the upper bound of the feasible set D can be expressed as $c^2=\tau mn{\sigma }^2$. If the noise variance ${\sigma }^2$ is unknown, then ${\sigma }^2$ can be estimated by referring to Reference [47]. $\tau$ is a predefined parameter, and $\tau =1$ is chosen empirically.

For the denoised image to retain more image detail, an error associated with the parameter $\lambda$ is set, which in turn allows the algorithm to adaptively select the parameter $\lambda$. The error is defined as

$$e_{k+1}=H{u}^{k+1}({\lambda }_{k+1})-f={({\lambda }_{k+1}t{H}^{T}H+I)}^{-1}(H\tilde{u}-f).$$

(50)

Whether $H^{T}H$ is a singular matrix or not, $({\lambda }_{k+1}t{H}^{T}H+I)$ is always invertible. According to the Lemma 3, the updated rules for $\lambda$ in the iterative process are as follows:

$${\lambda }^k=\left\{\begin{array}{cc}0\hfill & ,if\parallel H\tilde{u}{-f\parallel }_2^2\le c^2\hfill \\ therootof\parallel e^{k+1}{(\lambda )\parallel }_2^2=c^2\hfill & ,if\parallel H\tilde{u}{-f\parallel }_2^2>c^2\hfill \end{array}\right.$$

(51)

Lemma 3.

Let

$$\varphi (u,\lambda )={∥e∥}_2^2.$$

(52)

Then $\varphi (u,\lambda )$ is a convex function and strictly monotone decreasing. This function makes the following equation have a unique solution.

$$\varphi (u,\lambda )=c^2$$

Proof. 

The first-order derivative of the function $\varphi (u,\lambda )$ is $\frac{\partial \varphi (u,\lambda )}{\partial \lambda }=\frac{-2t{(H\tilde{u}-f)}^2}{{(t\lambda +I)}^3}\le 0$. From this, it can be seen that the function $\varphi (u,\lambda )$ is a monotonically decreasing function concerning the regularization parameter $\lambda$. The second-order derivative of the function $\varphi (u,\lambda )$ is $\frac{{\partial }^2\varphi (u,\lambda )}{\partial {\lambda }^2}=\frac{6t^2{(H\tilde{u}-f)}^2}{{(t\lambda +I)}^4}\ge 0$. It is evident that the function $\varphi (u,\lambda )$ is a strictly positive convex function about the regularization parameter $\lambda$. Therefore, the above equation has a unique solution. □

Contacting Algorithm 1, a new MOTV denoising algorithm can be obtained, see Algorithm 2.

Algorithm 2: Primal–Dual Algorithm for New MOTV Image Restoration

Require: f, H, $\beta$, $c{}^2$. 1: Initialize u, p, y and w. Set the size of s and t. 2: while stopping criterion is not satisfied do 3:     $C(u)={({\lambda }_1-{\lambda }_2)}^2={(k_{11}-k_{22})}^2+4k_{12}^2$; 4:     ${\xi }_1(u)=1-exp[-\frac{1}{{(C(u)/\gamma )}^2}]$; 5:     ${\xi }_2(u)=exp[-\frac{1}{{(C(u)/\gamma )}^2}]$; 6:     $p^{k+\frac{1}{2}}=\frac{p^k+s{\xi }_1{\nabla }^{\alpha }{u}^k}{max(1,|p^k+s{\xi }_1{\nabla }^{\alpha }{u}^k|)}$; 7:     $y^{k+\frac{1}{2}}=\frac{y^k+s{\xi }_2{\nabla }^2{u}^k}{max(1,|y^k+s{\xi }_2{\nabla }^2{u}^k|)}$; 8:     $w^{k+\frac{1}{2}}=\frac{w^k+s\beta u^k}{max(1,|w^k+s\beta u^k|)}$; 9:     $\tilde{u}=u^k-t{\xi }_1{({\nabla }^{\alpha })}^T{p}^{k+\frac{1}{2}}-t{\xi }_2{({\nabla }^2)}^T{y}^{k+\frac{1}{2}}-t\beta w^{k+\frac{1}{2}}$; 10:     if ${∥H\tilde{u}-f∥}_2^2<c^2$; then 11:         ${\lambda }_{k+1}=0$; 12:     else 13:         ${\lambda }_{k+1}$ is set to the root of $\varphi ({\lambda }^{k+1},u^k)=c^2$; 14:     end if 15:     $u^{k+1}={(t\lambda H^{T}H+I)}^{-1}(t\lambda H^{T}f+\tilde{u})$; 16:     $p^{k+1}=\frac{p^{k+\frac{1}{2}}+s{\xi }_1{\nabla }^{\alpha }{u}^{k+1}}{max(1,|p^{k+\frac{1}{2}}+s{\xi }_1{\nabla }^{\alpha }{u}^{k+1}|)}$; 17:     $y^{k+1}=\frac{y^{k+\frac{1}{2}}+s{\xi }_2{\nabla }^2{u}^{k+1}}{max(1,|y^{k+\frac{1}{2}}+s{\xi }_2{\nabla }^2{u}^{k+1}|)}$; 18:     $w^{k+1}=\frac{w^{k+\frac{1}{2}}+s\beta u^{k+1}}{max(1,|w^{k+\frac{1}{2}}+s\beta u^{k+1}|)}$; 19: end while 20: return $\mathbf{u}={\mathbf{u}}^{k+1}$.

3.4. Convergence Analysis

The sequence $u^{k+1}$ generated by Algorithm 2 converges to the minimum $u^{*}$ of $MTV(u)$. When u belongs to the feasible set D, the sequence ${\lambda }^k$ converges to the Lagrange multiplier ${\lambda }^{*}$. Through the primal problem (18) and the dual problem (22), $\Phi (u;{\lambda }^{*})={max}_{p\in P,y\in Y,w\in W}L(u,p,y,w;{\lambda }^{*})$ can be obtained. Therefore, when the dual variables satisfy the constraints $p\in P$, $y\in Y$, $w\in W$, it is only necessary to prove that the sequence $(u^k,p^k,y^k,w^k)$ generated by Algorithm 2 converges to a saddle point of the function $L(u,p,y,w;{\lambda }^{*})$. Firstly, two lemmas are introduced to aid in proving the convergence of the algorithm.

Lemma 4.

Let the sequence generated by Algorithm 2 be $(u^k,p^k,y^k,w^k,{\lambda }^k)$, then

$$\begin{array}{c}\hfill \parallel u^{k+1}-u^k{\parallel }_2^2\ge max(\frac{1}{2s{\xi }_{1}k(v_0^2+v_1^2+\cdots +v_{k-1}^2)}\langle u^{k+1}-u^k,\overline{{(-1)}^{\alpha }}{\mathrm{div}}^{\alpha }(p^{k+1}-p^{k+\frac{1}{2}})\rangle ,\\ \hfill \frac{1}{24s{\xi }_2}\langle u^{k+1}-u^k,{\mathrm{div}}^2(y^{k+1}-y^{k+\frac{1}{2}})\rangle ,\frac{1}{s\beta }\langle u^{k+1}-u^k,w^{k+1}-w^{k+\frac{1}{2}}\rangle ).\end{array}$$

(53)

Proof. 

Let $q_1=p^k+s{\xi }_1{\nabla }^{\alpha }{u}^{k+1}$, $q_2=p^k+s{\xi }_1{\nabla }^{\alpha }{u}^k$. According to the properties of $\mathrm{div}$ and ∇, the following can be obtained:

$$\begin{array}{cc}\hfill s{\xi }_1& \langle u^{k+1}-u^k,\overline{{(-1)}^{\alpha }}{\mathrm{div}}^{\alpha }(p^{k+1}-p^{k+\frac{1}{2}})\rangle \hfill \\ & =s{\xi }_1\langle {\nabla }^{\alpha }(u^{k+1}-u^k),{\rho }_p(q_1)-{\rho }_p(q_2)\rangle \hfill \\ & =\langle q_1-q_2,{\rho }_p(q_1)-{\rho }_p(q_2)\rangle .\hfill \end{array}$$

Using the classical inequality for any vector: $2a^{T}b\le {\parallel a\parallel }_2^2+{\parallel b\parallel }_2^2$, it can be obtained that

$$\begin{array}{cc}\hfill \langle & q_1-q_2,{\rho }_p(q_1)-{\rho }_p(q_2)\rangle \hfill \\ & \le \frac{1}{2}({∥q_1-q_2∥}_2^2+{∥{\rho }_p(q_1)-{\rho }_p(q_2)∥}_2^2)\hfill \\ & \le {∥q_1-q_2∥}_2^2={(s{\xi }_1)}^2{\parallel {\nabla }^{\alpha }(u^{k+1}-u^k)\parallel }_2^2.\hfill \end{array}$$

Let $v_i={(-1)}^i{C}_i^{\alpha }$, it can be obtained that

$$\begin{array}{cc}\hfill {\parallel {\nabla }^{\alpha }u\parallel }^2& =\sum _{ij}[{({\nabla }^{\alpha }{u}^1)}_{i,j}^2+{({\nabla }^{\alpha }{u}^2)}_{i,j}^2]\hfill \\ & =\sum _{ij}[{(v_0{u}_{i,j}^1+\cdots +v_{k-1}{u}_{i+k-1,j}^1)}^2+{(v_0{u}_{i,j}^2+\cdots +v_{k-1}{u}_{i,j+k-1}^2)}^2]\hfill \\ & \le 2k\sum _{ij}[{(v_0{u}_{i,j}^1)}^2+{(v_0{u}_{i,j}^2)}^2+\cdots +{(v_{k-1}{u}_{i+k-1,j}^1)}^2+{(v_{k-1}{u}_{i,j+k-1}^2)}^2]\hfill \\ & \le 2k(v_0^2+v_1^2+\cdots +v_{k-1}^2){\parallel u\parallel }^2.\hfill \end{array}$$

Consequently,

$$\begin{array}{c}\hfill \parallel u^{k+1}-u^k{\parallel }_2^2\ge \frac{1}{2s{\xi }_{1}k(v_0^2+v_1^2+\cdots +v_{k-1}^2)}\langle u^{k+1}-u^k,\overline{{(-1)}^{\alpha }}{\mathrm{div}}^{\alpha }(p^{k+1}-p^{k+\frac{1}{2}})\rangle ,\end{array}$$

where ${\xi }_1\ne 0$. Similarly, it follows that

$$\begin{array}{cc}\hfill \parallel u^{k+1}-u^k{\parallel }_2^2& \ge \frac{1}{24s{\xi }_2}\langle u^{k+1}-u^k,{\mathrm{div}}^2(y^{k+1}-y^{k+\frac{1}{2}})\rangle ,{\xi }_2\ne 0,\hfill \\ \hfill \parallel u^{k+1}-u^k{\parallel }_2^2& \ge \frac{1}{s\beta }\langle u^{k+1}-u^k,(w^{k+1}-w^{k+\frac{1}{2}})\rangle .\hfill \end{array}$$

The $\parallel {\nabla }^2{u\parallel }_2^2$ appearing in the solution process is expressed as

$$\begin{array}{cc}\hfill \parallel {\nabla }^2{u\parallel }_2^2& =\sum _{ij}{({\nabla }_1^2{u}_{i,j})}^2+{({\nabla }_2^2{u}_{i,j})}^2\hfill \\ \hfill & =\sum _{ij}{(u_{i+1,j}-2u_{i,j}+u_{i-1,j})}^2+{(u_{i,j+1}-2u_{i,j}+u_{i,j-1})}^2\hfill \\ \hfill & \le 2\sum _{ij}(u_{i+1,j}^2+u_{i-1,j}^2+8u_{i,j}^2+u_{i,j+1}^2+u_{i,j-1}^2)\hfill \\ \hfill & \le {24\parallel u\parallel }_2^2.\hfill \end{array}$$

In summary,

$$\begin{array}{c}\hfill \parallel u^{k+1}-u^k{\parallel }_2^2\ge max(\frac{1}{2s{\xi }_{1}k(v_0^2+v_1^2+\cdots +v_{k-1}^2)}\langle u^{k+1}-u^k,\overline{{(-1)}^{\alpha }}{\mathrm{div}}^{\alpha }(p^{k+1}-p^{k+\frac{1}{2}})\rangle ,\\ \hfill \frac{1}{24s{\xi }_2}\langle u^{k+1}-u^k,{\mathrm{div}}^2(y^{k+1}-y^{k+\frac{1}{2}})\rangle ,\frac{1}{s\beta }\langle u^{k+1}-u^k,w^{k+1}-w^{k+\frac{1}{2}}\rangle ),\end{array}$$

which is proved. □

Lemma 5.

Let the sequence generated by Algorithm 2 be $(u^k,p^k,y^k,w^k,{\lambda }^k)$, then

$$\begin{array}{cc}\hfill \frac{1}{t}& \parallel u^{*}-u^k{\parallel }_2^2+\frac{1}{s}\parallel p^{*}-p^k{\parallel }_2^2+\frac{1}{s}\parallel y^{*}-y^k{\parallel }_2^2+\frac{1}{s}{\parallel w^{*}-w^k\parallel }_2^2\hfill \\ \hfill & \ge \frac{1}{t}\parallel u^{*}-u^{k+1}{\parallel }_2^2+\frac{1}{s}\parallel p^{*}-p^{k+1}{\parallel }_2^2+\frac{1}{s}\parallel y^{*}-y^{k+1}{\parallel }_2^2+\frac{1}{s}{\parallel w^{*}-w^{k+1}\parallel }_2^2.\hfill \end{array}$$

(54)

In particular, the sequence $(u^k,p^k,y^k,w^k,{\lambda }^k)$ converges to the limit point $(u^{★},p^{★},y^{★},w^{★},{\lambda }^{★})$, at which point there are two cases, if $u^{★}$ belongs to the feasible set D, then ${\lambda }^{★}=0$. Conversely, ${\lambda }^{★}$ takes the solution of the equation $\parallel H{u}^{★}{(\lambda )-f\parallel }_2^2=c^2$.

See [40] for a related proof.

Next, prove that a saddle point of $L(u,p,y,w;{\lambda }^{★})$ is $(u^{★},p^{★},y^{★},w^{★})$, when $p^{★}\in P$, $y^{★}\in Y$, $w^{★}\in W$, it follows from Algorithm 2 that $(u^{★},p^{★},y^{★},w^{★})$ satisfies the following equation:

$$u^{★}={(t{\lambda }^{★}{H}^{T}H+I)}^{-1}(t{\lambda }^{★}{H}^{T}f+u^{★}-t{\xi }_1{({\nabla }^{\alpha })}^T{p}^{★}-t{\xi }_2{({\nabla }^2)}^T{y}^{★}-t\beta w^{★}).$$

(55)

Further collation yields,

$${\lambda }^{★}{H}^T(H{u}^{★}-f)+{\xi }_1{({\nabla }^{\alpha })}^T{p}^{★}+{\xi }_2{({\nabla }^2)}^T{y}^{★}+\beta w^{★}=0.$$

(56)

This process is also the process of replacing $\lambda$ with ${\lambda }^{★}$. Starting from Algorithm 2, it is known that $(u^{★},p^{★},y^{★},w^{★})$ also satisfies the following equation system:

$$\left\{\begin{array}{c}{p}^{★}={\rho }_P(p^{★}+s{\xi }_1{\nabla }^{\alpha }{u}^{★}),\hfill \\ y^{★}={\rho }_Y(y^{★}+s{\xi }_2{\nabla }^2{u}^{★}),\hfill \\ w^{★}={\rho }_W(w^{★}+s\beta u^{★}).\hfill \end{array}\right.$$

(57)

The solution of the projection operator is shown in (38), which thus yields

$$\left\{\begin{array}{cc}\hfill p_{i,j}^{★}& =\frac{p_{i,j}^{★}+s{\xi }_1{({\nabla }^{\alpha }{u}^{★})}_{i,j}}{max(1,|p_{i,j}^{★}+s{\xi }_1{({\nabla }^{\alpha }{u}^{★})}_{i,j}|)},\hfill \\ \hfill y_{i,j}^{★}& =\frac{y_{i,j}^{★}+s{\xi }_2{({\nabla }^2{u}^{★})}_{i,j}}{max(1,|y_{i,j}^{★}+s{\xi }_2{({\nabla }^2{u}^{★})}_{i,j}|)},\hfill \\ \hfill w_{i,j}^{★}& =\frac{w_{i,j}^{★}+s\beta u_{i,j}^{★}}{max(1,|w_{i,j}^{★}+s\beta u_{i,j}^{★}|)}.\hfill \end{array}\right.$$

(58)

If $|p_{i,j}^{★}+s{\xi }_1{({\nabla }^{\alpha }{u}^{★})}_{i,j}|\le 1$, then ${({\nabla }^{\alpha }{u}^{★})}_{i,j}=0$ and $|p_{i,j}^{★}|\le 1$. If $|p_{i,j}^{★}+s{\xi }_1{({\nabla }^{\alpha }{u}^{★})}_{i,j}|>1$, then $|p_{i,j}^{★}|=1$ and $s{\xi }_1{({\nabla }^{\alpha }{u}^{★})}_{i,j}=(1-|p_{i,j}^{★}+s{\xi }_1{({\nabla }^{\alpha }{u}^{★})}_{i,j}|)p_{i,j}^{★}$. Similarly, $y_{i,j}^{★}$, $w_{i,j}^{★}$ has a similar property. Therefore, the $p^{★}\in P,y^{★}\in Y,w^{★}\in W$ and the KKT optimality condition of (32) hold.

Thus, the following theorem can be obtained.

Theorem 1.

When ${\lambda }^{*}$ is the Lagrange multiplier on the constraint $u\in D$, assuming that the saddle point of $L(u,p,y,w;{\lambda }^{*})$ is $(u^{*},p^{*},y^{*},w^{*})$, then the sequence $(u^k,p^k,y^k,w^k;{\lambda }^k)$ generated by Algorithm 2 converges to $(u^{*},p^{*},y^{*},w^{*};{\lambda }^{*})$. At this time, $s·t\le 1/16$. In particular, $u^k$ converges to the minimum, ${\lambda }^k$ converges to the Lagrange multiplier satisfying the constraint $u\in D$.

4. Simulation Results

4.1. Quantitative Evaluation Index Description

In this paper, the peak signal-to-noise ratio (PSNR) and structural self-similarity (SSIM) are used to assess the quality of denoised images. PSNR is used to evaluate the quality similarity between two images, while SSIM measures the structural similarity between two images. Here, the PSNR is defined as

$$\mathrm{PSNR}=10{log}_{10}(\frac{{\mathrm{b}}^2}{\mathrm{MSE}}).$$

(59)

where $b=max(m,n)$ and $m,n$ denotes the size of the image. MSE denotes the mean square error between the two images, which is defined as

$$\mathrm{MSE}=\frac{1}{\mathrm{mn}}\sum _{\mathrm{i}=0}^{\mathrm{m}-1}\sum _{\mathrm{j}=0}^{\mathrm{n}-1}({[{\mathrm{u}}_0(\mathrm{i},\mathrm{j})-\mathrm{u}(\mathrm{i},\mathrm{j})]}^2).$$

(60)

Here, $u_0,u$ represents the original image and the denoised image, respectively. The definition of SSIM is as follows:

$$\mathrm{SSIM}({\mathrm{u}}_0,\mathrm{u})=\frac{(2{\mathcal{A}}_{{\mathrm{u}}_0}{\mathcal{A}}_{\mathrm{u}}+{\mathrm{d}}_1)(2cov({\mathrm{u}}_0,\mathrm{u})+{\mathrm{d}}_2)}{({\mathcal{A}}_{{\mathrm{u}}_0}^2+{\mathcal{A}}_{\mathrm{u}}^2+{\mathrm{d}}_1)({\mathcal{B}}_{{\mathrm{u}}_0}^2+{\mathcal{B}}_{\mathrm{u}}^2+{\mathrm{d}}_2)}.$$

(61)

where $\mathcal{A}$ denotes the mean, $cov(u{}_0,u)$ denotes the covariance, $\mathcal{B}$ denotes the standard deviation, and $d_1$ and $d_2$ are two constants, and it is commonly recommended based on experience to set $d_1=0.01$ and $d_1=0.03$. To verify the feasibility of the algorithm, a set of images is first introduced as test images, as depicted in Figure 2.

4.2. The Influence of Fractional-Order $\alpha$

The regularization parameter $\lambda$ and the fractional-order $\alpha$ are both indefinite parameters in image denoising algorithms, but they both play a significant role. The effect of the fractional-order $\alpha$ on the performance of the algorithm is investigated by adding Gaussian noise with the noise standard deviation of $\sigma =10,20$ to the test images Lena, Boats, and Pepper. From Figure 3, it can be observed that when $\sigma =10$, with the increase of $\alpha$, PSNR shows a decreasing trend. When $\sigma =20$, with the increase of $\alpha$, PSNR shows an increasing trend. The value near $\alpha =1.5$ can be considered to ensure that the changing trend of PSNR has certain adaptability to different noise standard deviations. From Figure 4, it can be observed that if $\alpha$ is small, the image is prone to a blocky effect, and if $\alpha$ is large, the image is not well denoised. To effectively suppress the phenomenon observed in Figure 4, it is also possible to consider selecting values around $\alpha =1.5$. In conclusion, it is preferable to select $\alpha$ from the range of $\alpha \in [1.4,1.6]$.

4.3. The Influence of $\lambda$ on Algorithm Performance

The regularization parameter $\lambda$ is also known to affect the denoising performance of the algorithm. Hence, Gaussian noise with $\sigma =30,40$ is added to the Lena, Boats, and Pepper images to verify the impact of the regularization parameter $\lambda$ on the performance of the algorithm. From Figure 5, it can be observed that when manually selecting the parameter $\lambda$, the PSNR and SSIM of each algorithm are not satisfactory. However, when using the adaptive parameter $\lambda$ selection method proposed in Section 3.3, the algorithms achieve higher PSNR and SSIM. Therefore, it can be concluded that the adaptive parameter $\lambda$ proposed can improve the performance of the algorithm to a certain extent.

4.4. Comparison with Other Methods

To further prove the denoising performance of the new algorithm, a first set of comparative experiments was conducted using test images Parrots, Lena, Barbara, Cameraman, Pepper, Foot, and Head, in comparison with some state-of-the-art methods [48,49,50]. Gaussian noise with $\sigma =5,10,20,30$ was added to the images, and the PSNR and SSIM of each algorithm after denoising are shown in

Table 1. The PSNR and SSIM of each algorithm under the first set of comparative experiments.

Image	$\mathit{\sigma }$	PSNR/SSIM [48]	PSNR/SSIM [49]	PSNR/SSIM [50]	PSNR/SSIM Proposed
Parrots	5	32.8362/0.9071	32.3292/0.9050	34.3292/0.9058	34.8288/0.9080
	10	27.7650/0.9085	31.4979/0.9005	33.4346/0.9089	33.6375/0.9101
	20	23.9519/0.9146	28.2574/0.9162	28.5716/0.9218	28.5996/0.9241
	30	22.7536/0.92018	23.1230/0.9204	24.2050/0.9277	26.3916/0.9279
Lena	5	33.1229/0.9023	32.2596/0.8938	33.1639/0.8992	33.4520/0.9000
	10	27.7962/0.8963	31.4268/0.8974	32.4873/0.9023	32.5379/0.9011
	20	24.0047/0.9144	28.0919/0.9050	28.1660/0.9141	28.4511/0.9153
	30	23.0918/0.8942	23.9741/0.8973	24.2123/0.9075	25.8134/0.8896
Barbara	5	33.1016/0.8908	29.2736/0.8904	31.0013/0.8909	31.0157/0.8924
	10	27.7793/0.8954	28.8140/0.8973	30.0985/0.8974	30.1247/0.8978
	20	23.9799/0.9024	26.7437/0.9033	26.5292/0.9051	26.5692/0.9133
	30	22.1423/0.8937	22.9751/0.9102	22.6868/0.9188	22.9827/0.9124
Cameraman	5	33.1716/0.9021	29.8097/0.8991	31.8976/0.9000	32.2383/0.9006
	10	27.8417/0.9056	29.3846/0.9015	30.3292/0.9034	31.4454/0.9059
	20	23.9418/0.9140	27.0952/0.9084	27.6781/0.9082	27.7247/0.9144
	30	22.8149/0.8902	23.1143/0.8907	23.5803/0.8905	24.6137/0.8911
Pepper	5	34.6851/0.8698	32.6693/0.8620	34.5614/0.8713	34.7247/0.8729
	10	28.8834/0.8764	32.0109/0.8689	32.9939/0.8786	33.8284/0.8767
	20	24.2617/0.8853	28.6892/0.8808	29.1745/0.8902	29.3960/0.8913
	30	22.5738/0.8937	24.4121/0.8938	25.4702/0.8950	29.0906/0.9016
Foot	5	37.5861/0.8536	37.4102/0.8524	35.4676/0.8713	35.7180/0.8561
	10	33.8384/0.8340	33.9815/0.8312	34.0148/0.8311	34.6360/0.8350
	20	27.1051/0.8286	28.2964/0.8279	28.7614/0.8253	29.8339/0.8272
	30	25.4720/0.8050	25.6541/0.8144	25.9459/0.8193	27.5594/0.8074
Head	5	33.5699/0.8953	33.4973/0.8874	34.6343/0.9030	33.6142/0.8423
	10	32.3203/0.8470	31.0355/0.8309	32.1153/0.8576	32.3258/0.8310
	20	22.9786/0.7933	22.8761/0.78237	27.6259/0.8117	28.1583/0.7446
	30	19.4785/0.7403	19.4509/0.7371	23.3897/0.7649	24.5390/0.6468

. Figure 6 and Figure 7 demonstrate the denoising effects of the Pepper and Foot images with $\sigma =30$. Overall, the algorithm proposed can achieve better results.

Next, a second set of comparative experiments [51,52,53] was conducted by adding Gaussian noise with $\sigma =10,20,30$ to Cameraman, Boats, Barbara, and Pepper. The experimental results are shown in

Table 2. The PSNR and SSIM of each algorithm under the second set of comparative experiments.

Image	$\mathit{\sigma }$	PSNR/SSIM [51]	PSNR/SSIM [52]	PSNR/SSIM [53]	PSNR/SSIM Proposed
Cameraman	10	31.3925/0.8948	33.1347/0.9049	25.0623/0.9089	33.6375/0.9101
	20	28.2944/0.8272	28.5723/0.8309	24.9803/0.9218	28.5996/0.9241
	30	27.1039/0.7700	27.2705/0.7315	24.8771/0.9277	26.3916/0.9279
Boats	10	32.3014/0.9017	32.3438/0.9065	25.6609/0.9032	32.4058/0.8889
	20	28.6036/0.8802	28.7042/0.8840	25.5432/0.9019	28.3939/0.9011
	30	27.5687/0.8204	27.5816/0.8253	25.3779/0.9087	26.0102/0.9112
Barbara	10	30.0487/0.9186	31.8501/0.9340	22.0302/0.8974	30.1247/0.8978
	20	26.1845/0.9259	26.4952/0.9247	21.9840/0.9051	26.5692/0.9133
	30	25.7061/0.8790	25.7384/0.8790	21.9057/0.9188	22.9827/0.9124
Pepper	10	33.3949/0.8950	33.9108/0.9032	28.8618/0.9034	34.7247/0.8729
	20	29.2765/0.8593	29.5572/0.8897	28.6174/0.8596	29.3960/0.8913
	30	28.2187/0.8746	28.2854/0.8731	28.4424/0.8905	29.0906/0.9016

. Based on the experimental results, it can be observed that the results obtained in this study are slightly higher than those of other algorithms. In conclusion, the algorithm proposed is effective in a certain sense.

5. Conclusions

In summary, the achievements obtained in this article are as follows:

*

The proposed MOTV model plays a significant role in improving the staircase effect and eliminating residual noise.

*

The appropriate fractional-order $\alpha$ and the selection of $\lambda$ based on the deviation principle have a certain effect on the performance of the algorithm.

*

According to the results of two groups of comparative experiments, it can be concluded that the algorithm proposed outperforms other algorithms and is effective.

In fact, the algorithm proposed also has certain limitations. Specifically, when the noise variance increases, the noise suppression effect of the algorithm may decrease. In future research, it will be considered by the authors to incorporate the ideas of this paper into other methods mentioned in the introduction section. Additionally, the application of regularizers will be considered to achieve better results. Furthermore, efforts will be made to reduce the computation time.