Proximity algorithms for the ${\mathit{L}}^{1}{\mathit{L}}^{2}/{\mathit{T}\mathit{V}}^{\mathit{\alpha }}$ image denoising model

1.   Introduction

Throughout the evolution of digital image processing, a variety of processing technologies have been formed, including the wavelet transform, partial differential equation (PDE), and stochastic model. In image processing, the edge of an image is the most important visual feature. In 1992, Rudin et al. proposed the well-known total variation (TV) model ^[1], which has been named the ROF model. The ROF model can balance edge preservation and noise removal because it can take advantage of the inherent regularity of the image. The ROF model is as follows:

where $\mathrm{\Omega }\subset {R}^{n}$ is an open bounded set, $n\ge 2$ ^[2], ${u}_{0}(x, y)$ denotes the noisy image and $u(x, y)$ denotes the desired clean image. λ denotes a real positive number and ${‖u‖}_{TV}$ denotes the TV of $u(x, y)$ , which is defined as ${‖\nabla u‖}_{1}$ . The ROF model has played an important role in image denoising, deblurring and inpainting. However, the solution of the ROF model is a piecewise constant function, so it is easy to generate a blocky effect in the flat region. To reduce the block effect, scholars have proposed a fourth-order PDE ^[3] and LLT ^[4], which can effectively remove the noise and reduce the blocky effect. The LLT model is as follows:

where $\mathrm{\Delta }u = ({\partial }_{x}^{2}u, {\partial }_{y}^{2}u)$ and ${‖\mathrm{\Delta }u‖}_{1} = \left|{\partial }_{x}^{2}u\right|+\left|{\partial }_{y}^{2}u\right|$ . The disadvantage of the LLT model is that it produces excessive smoothing in the edge region. To solve this problem, an adaptive fourth-order PDE has been proposed ^[5]. Both the ROF model and LLT model have the ${L}^{2}$ -data fidelity term. The type of noise that corrupts the image typically affects the data fidelity term selection. In general, images are affected by different types of noise. If the image is only affected by a mixture of Gaussian noise and Poisson noise, the noises can be converted into additive Gaussian noise. This is probably why most of the literature is devoted to removing Gaussian noise. The ${L}^{2}$ - data fidelity term is suitable for removing Gaussian additive noise, but it is almost invalid for other noises. The ${L}^{1}$ -data fidelity term can effectively remove non-additive Gaussian noise, such as Laplacian noise and impulse noise ^[6,7]. The ${L}^{1}/TV$ model is as follows:

The ${L}^{1}/TV$ model has some unique features. It does not destroy the geometric structures or morphological invariance of the images under processing ^[8,9]. Therefore, the ${L}^{1}/TV$ denoising image model is widely used in practical applications, such as face recognition ^[10], shape denoising ^[11] and image texture decomposition ^[12]. In fact, images are generally not corrupted by only one type of noise. The mixture of Gaussian and salt-and-pepper noise is considered in this paper. In particular, salt-and-pepper noise is a simple type of impulse noise ^[13]. An ${L}^{1}$ - ${L}^{2}$ -data fidelity term was introduced and proved to be suitable for the removal of a mixture of Gaussian and impulse noise in ^[14]. The ${L}^{1}{L}^{2}/TV$ model is as follows:

where $\lambda$ , $\mu \ge 0$ . The ${L}^{1}{L}^{2}/TV$ model (4) is a generalization of (1) and (3). For example, if we set $\lambda = 0$ in (4) then we get the ${L}^{1}/TV$ model. If we set $\mu = 0$ then we get the ${L}^{2}/TV$ model. In particular, the choice of parameters critically affects the quality of image restoration. Small values of $\lambda$ and $\mu$ lead to an oversmoothed reconstruction, which eliminates both noise and detail in the image. In contrast, large values of $\lambda$ and $\mu$ retain noise ^[15]. An improvement of the ${L}^{1}{L}^{2}/TV$ model has been proposed in ^[16], where ${‖Wu‖}_{1}$ replaces the TV. In ^[17], the authors used second-order total generalized variation ^[18] as a regularization term and incorporated box constraints.

In this paper, the fractional-order TV regularization term is the focus. We propose a combined model with a fractional-order TV regularization term, an ${L}^{1}$ -data fidelity term, and an ${L}^{2}$ - data fidelity term, which we term the ${L}^{1}{L}^{2}/{TV}^{\alpha }$ model. This model aims to remove the mixture of Gaussian noise and salt-and-pepper noise.

It is difficult to minimize the objective function because the fractional-order TV regularization term is non-differentiable. Numerous efforts have been devoted to addressing this issue. There are some methods to solve the fractional-order TV model, including the use of the primal-dual algorithm ^[19], fractional-order Euler-Lagrange equations ^[20], alternating projection algorithm for the fractional-order multi-scale variational model ^[21,22], and majorization-minimization algorithm ^[23]. The Split Bregman iterative algorithm ^[24] and alternating direction method of multipliers ^[25] can also effectively solve non-differentiable terms. Recently, proximity algorithms ^[26–30] for solving the ROF model or the ${L}^{1}/TV$ denoising image model have attracted widespread attention in digital image processing. The method mainly combines a convex function with a linear transformation to represent the non-differentiable term ${\Vert u\Vert }_{TV}$ . The issue of solving the proximity operator of the convex function can be reformulated into solving a fixed-point equation. Consequently, the proximity operator of the convex function can be obtained. The convergence of the fixed-point proximity algorithm has been proven ^[26]. The ${L}^{1}/TV$ model requires solving two fixed point equations due to the non-differentiability of the ${\mathrm{L}}^{1}$ -data fidelity term ^[28]. In this paper, the proximity algorithm is used to solve the ${L}^{1}{L}^{2}/{TV}^{\alpha }$ model.

The structure of the paper is as follows. Section 1 introduces the prior works and our motivation. Section 2 proposes the ${L}^{1}{L}^{2}/{TV}^{\alpha }$ model and proves the existence of its solution. The proximity algorithm is applied to solve the model and the convergence of the algorithm is proved. Section 3 presents several numerical experiments and shows the results. Finally, Section 4 concludes the paper.

2.   ${\mathit{L}}^{1}{\mathit{L}}^{2}/\mathit{T}{\mathit{V}}^{\mathit{a}}$ denoising model and proximity algorithm

2.1.   Proximity operator

This section first introduces two very important concepts of convex functions: the proximity operator and the subdifferential. The relationship between them will also be given.

Initially, we introduce some notations. We denote the $m$ -dimensional Euclidean space by ${R}^{m}$ . For $x, y\in {R}^{m}$ , we define the standard inner product of ${R}^{m}\;\mathrm{as} < x, y > ≔\sum _{i = 1}^{m}{x}_{i}{y}_{i}\;\mathrm{and}\; \mathrm{the}\;p-\mathrm{norm}\; \mathrm{of} \;\mathrm{a}\; \mathrm{vector}\;x\in {R}^{m}$ as ${\left|\left|x\right|\right|}_{p}≔{\left({\sum _{i = 1}^{m}{|x}_{i}|}^{p}\right)}^{\frac{1}{p}}$ . The proximity operator was introduced in ^[31]. We recall its definition as follows.

Definition 2.1. (Proximity operator): Let $f$ be a proper lower-semi-continuous convex function on ${R}^{m}$ , where ${R}^{m}$ is $m$ -dimensional Euclidean space. The proximity operator of $f$ is defined for any $x\in {R}^{m}$ by $pro{x}_{f}\left(x\right) = \underset{u}{\mathrm{arg}\mathit{min}}\{\frac{1}{2}{‖u-x‖}_{2}^{2}+f\left(u\right):u\in {R}^{m}\}$ .

Definition 2.2. (Subdifferential): Let $f$ be a proper lower-semi-continuous convex function on ${R}^{m}$ , where ${R}^{m}$ is $m$ -dimensional Euclidean space. The subdifferential of $f$ is defined for $y\in {R}^{m}$ by $\partial f\left(x\right)≔\{y\in {R}^{m}andf\left(z\right)\ge f\left(x\right)+\left\langle{y, z-x}\right\rangle, \forall z\in {R}^{m}\}$ .

The following lemma describes the relationship between the proximity operator and the convex function subdifferential.

Lemma 2.1. (Proposition 2.6 in ^[27]): If $f$ is a convex function on ${R}^{m}$ and $x\in {R}^{m}$ , then

The proof of this lemma is given in ^[27]. Based on the Lemma 2.1, we can get that

2.2.   Model description and analysis

Recently in ^[13], it has been demonstrated that the ${L}^{1}{L}^{2}/TV$ model is effective at removing mixtures of Gaussian and impulse noise. In this approach, an image is restored by solving the following equation:

where ${p}_{0}\in {R}^{N\times N}$ denotes the noise image, $N$ is a positive integer, $p\in {R}^{N\times N}$ denotes the denoising image, and $\lambda, \mu$ are the parameters of ${L}^{2}$ -data and ${L}^{1}$ -data fidelity terms respectively. This model combines two kinds of data fidelity terms, ${L}^{1}$ and ${L}^{2}$ , which can combine the advantages of both norms. Therefore, it has a significant effect in the removal of mixtures noise of Gaussian noise and salt-and-pepper noise.

However, we observe that the numerical solution produced by the ${L}^{1}{L}^{2}/TV$ model yields a substantial block effect. Additionally, this model fails to completely remove salt-and-pepper noise. The fractional-order TV regularization term has been proved to effectively reduce the block effect. This section introduces a minimum optimization denoising model, termed the ${L}^{1}{L}^{2}/T{V}^{\alpha }$ model. The ${L}^{1}{L}^{2}/T{V}^{\alpha }$ model includes three terms: an ${L}^{2}$ -data fidelity term for Gaussian noise, an ${L}^{1}$ -data fidelity term for salt-and-pepper noise, and a fractional-order TV regularization term for a balance between detail preservation and noise reduction. The model is as follows:

where ${p}_{0}\in {R}^{N\times N}$ denotes the noise image and $p\in {R}^{N\times N}$ denotes the denoising image. ${‖p‖}_{T{V}^{\alpha }}$ is the $\alpha$ fractional-order TV of $p$ , and ${‖p‖}_{T{V}^{\alpha }}$ is defined as ${‖{\nabla }^{\alpha }p‖}_{1}$ , where ${\nabla }^{\alpha }p = ({\partial }_{x}^{\alpha }p, {\partial }_{y}^{\alpha }p)$ and ${‖{\nabla }^{\alpha }p‖}_{1} = \left|{\partial }_{x}^{\alpha }p\right|+\left|{\partial }_{y}^{\alpha }p\right|$ . In particular, note the following:

● When setting $\lambda = 0$ , the model (8) simplifies ${L}^{1}/T{V}^{\alpha }$ .

● When setting $\mu = 0$ , the model (8) simplifies ${L}^{2}/T{V}^{\alpha }$ .

The parameter settings of $\lambda$ and $\mu$ for these specialized models demonstrates the flexibility of the ${L}^{1}{L}^{2}/T{V}^{\alpha }$ model.

To prove the existence of a solution to the ${L}^{1}{L}^{2}/T{V}^{\alpha }$ model, it is critical to prove the boundedness of the potential solution ^[33].

Lemma 2.2. (Boundedness) Let ${p}_{0}\in {L}^{2}\left(\varOmega \right)$ , where $\mathrm{\Omega }\subset {R}^{n}(n\ge 2)$ is an open bounded set. Given $\underset{\varOmega }{inf}\;{p}_{0} > 0$ , if the model has a solution $\widehat{p}$ , then $\underset{\varOmega }{inf}\;{p}_{0} < \widehat{p} < \underset{\varOmega }{sup}\;{p}_{0}$ .

Proof of Lemma 2.2. Let $\omega = \underset{\varOmega }{inf}\;{p}_{0}$ and $\nu = \underset{\varOmega }{sup}\;{p}_{0}$ . When $p > {p}_{0}$ , functions $\left|p-{p}_{0}\right|$ and ${(p-{p}_{0})}^{2}$ increase monotonically. Then,

where $inf(p, \nu)$ is the lower bound of $p$ and $\nu$ . That is, $inf\left(p, \nu \right)$ is the minimum value of $p$ and $\nu$ .

Moreover, based on Lemma 2 in the literature ^[34], there exists ${TV}^{\alpha }\left(\mathit{inf}\left(p, \nu \right)\right)\le {TV}^{\alpha }\left(p\right)$ . Thus, we have

and the equation holds if and only if $p\le \nu$ .

Since $\widehat{p}$ is the minimum solution of optimization problem (8), the equation holds when $p = \widehat{p}$ and hence $\widehat{p}\le \nu$ . Similarly, $E\left(p\right)\le E\left(sup\right(p, \omega \left)\right)$ ; then, $\widehat{p}\ge \omega$ can be obtained. In summary, $\underset{\varOmega }{inf}\;f < \widehat{p} < \underset{\varOmega }{sup}\;f$ .

In what follows, we will give the existence of a solution for the optimization problem (8).

Lemma 2.3. (Existence): Let ${p}_{0}\in {L}^{2}\left(\varOmega \right)$ , where $\mathrm{\Omega }\subset {R}^{n}$ ( $n\ge 2$ ) is an open bounded set. Given $\underset{\varOmega }{inf}\;{p}_{0} > 0$ , the optimization problem (8) has at least one solution in the solution space ${BV}^{\alpha }\left(\varOmega \right)$ .

Proof of Lemma 2.3. The space of bounded variational functions ${BV}^{\alpha }\left(\varOmega \right)$ can be defined as follows: ${BV}^{\alpha }\left(\varOmega \right) = \{f:f\in {L}^{1}\left(\varOmega \right)\}$ , forming Banach spaces under the ${BV}^{\alpha }$ norm ${‖f‖}_{{BV}^{\alpha }} = {‖f‖}_{{L}^{1}}+{TV}^{\alpha }\left(f\right)$ .

Define $\omega = \underset{\varOmega }{inf}\;{p}_{0}$ and $\nu = \underset{\varOmega }{sup}\;{p}_{0}$ . Because $p = \nu \in {BV}^{\alpha }\left(\mathrm{\Omega }\right)$ , the solution space is not empty ^[35]. Consider that the optimization problem (8) has a minimization sequence $\left\{{p}_{n}\right\}\in {BV}^{\alpha }\left(\varOmega \right)$ with $\omega \le {p}_{n}\le \nu$ .

Because ${BV}^{\alpha }\left(\varOmega \right)$ is a Banach space and $\mathrm{\Omega }$ is bounded, it follows that

Moreover, because $\left\{{p}_{n}\right\}$ is a minimization sequence, there exists a constant $C > 0$ such that $E\left({p}_{n}\right)\le C$ . Because ${\int }_{\varOmega }\left|\right|p-{p}_{0}|{|}_{2}^{2}+||p-{p}_{0}|{|}_{1}d\varOmega$ is nonnegative, there is a constant $C\mathrm{\text{'}} > 0$ and

Equations (12) and (13) yield that $\left\{{p}_{n}\right\}$ is consistently bounded. Due to the compactness of ${BV}^{\alpha }\left(\varOmega \right)$ , there exists a subsequence $\left\{{p}_{{n}_{j}}\right\}$ of $\left\{{p}_{n}\right\}$ and a function ${pϵBV}^{\alpha }\left(\varOmega \right)$ such that

Using the Lebesgue control convergence theorem, we obtain

According to the lower semi-continuity of the function, the following inequality holds:

Since $\left\{{p}_{n}\right\}$ is a minimization sequence, $p$ is the smallest solution to the optimization problem (8).

2.3.   Numerical algorithm and convergence analysis

Consider an image represented by a grid of $N\times N$ pixels. The discretization of the data term is given by

where $(i, j)$ denotes the coordinates at the points. For the fractional-order TV term, we obtain the following discretization:

where ${C}_{k}^{\left(\alpha \right)} = {\left(-1\right)}^{\alpha }\frac{\mathrm{\Gamma }(\alpha +1)}{\mathrm{\Gamma }\left(k+1\right)\mathrm{\Gamma }(\alpha -k+1)}$ and $\mathrm{\Gamma }\left(x\right)$ is the gamma function.

Considering that the proximity algorithm is suitable for vectors, we respectively transform the image matrices $p$ and ${p}_{0}$ into vectors $u$ and ${u}_{0}$ by using the formulas ${p}_{i, j} = {u}_{\left(j-1\right)n+i}$ and ${p}_{0i, j} = {{u}_{0}}_{\left(j-1\right)n+i}$ , $i, j = \mathrm{1, 2}, \dots., N$ . We describe the minimization problem (8) as follows:

where $u\in {R}^{m}$ and ${u}_{0}\in {R}^{m}, m = {N}^{2}$ .

The proximity operator of ${‖{\nabla }^{\alpha }u‖}_{1}$ is not easy to compute. To overcome this difficulty, we treat ${‖{\nabla }^{\alpha }u‖}_{1}$ as the composition of a convex function with a fractional-order difference operator by using the formula ${‖{\nabla }^{\alpha }u‖}_{1} = (\phi \circ {B}^{\alpha })\left(u\right)$ . In the formula, $\phi :{R}^{2m}\to R$ is defined as the norm ${‖\cdot ‖}_{1}, {B}^{\alpha }$ is a $2m\times m$ matrix, and ${\nabla }^{\alpha }u$ can be represented as ${B}^{\alpha }u$ . The $(i, j)$ component of ${\nabla }^{\alpha }u$ can thus be represented as a multiplication of the vector $u\in {R}^{m}$ by a matrix ${B}_{n}^{\alpha }\in {R}^{2\times m}$ for $n = 1, 2, ..., m$ :

where the matrix ${B}_{n}^{\alpha } = {\left[{B}_{1}^{\alpha }, {B}_{2}^{\alpha }, \dots, {B}_{N}^{\alpha }\right]}^{T}\in {R}^{2\times 2m}$ ^[29]. Therefore, we describe the minimization problem as follows:

Consider $\varphi$ to be a convex function on ${R}^{m}$ at $u\in {R}^{m}$ , as follows:

Therefore, we can describe the above minimization problem as follows:

Proposition 2.1. Let $\phi$ be a proper convex function on ${R}^{m}$ ; ${B}^{\alpha }$ is a $2m\times m$ matrix. If $u\in {R}^{m}$ is a solution of model (21), then for any positive numbers ${\beta }_{1},$ ${\beta }_{2} > 0$ , there exists a vector $b\in {R}^{2m}$ such that

On the contrary, if $b\in {R}^{2m}$ and $u\in {R}^{m}$ satisfies (22) and (23) for some positive ${\beta }_{1}$ , ${\beta }_{2} > 0$ , then $u$ is a solution of (21).

Proof. If $u\in {R}^{m}$ is a solution of (21), then, by Fermat's theorem on convex analysis, it follows that

By the chain rule

then

For any ${\beta }_{1}$ , ${\beta }_{2} > 0$ , we choose two vectors $a\in \frac{1}{{\beta }_{1}}\partial \varphi \left(u\right)$ and $b\in \frac{1}{{\beta }_{2}}\partial \phi \left({B}^{a}u\right)$ such that

By (5) and $a\in \frac{1}{{\beta }_{1}}\partial \varphi \left(u\right)$ , we have that

Using (25), we conclude that $a = -\frac{{\beta }_{2}}{{\beta }_{1}}{\left({B}^{\alpha }\right)}^{T}b$ ; by substituting $a$ into (26), we obtain (22). By applying the definition of the proximity operator and $b\in \frac{1}{{\beta }_{2}}\partial \phi \left({B}^{\alpha }u\right)$ , we obtain (23). Conversely, if there exist ${\beta }_{1}, {\beta }_{2} > 0$ , $b\in {R}^{2m}$ , and $u\in {R}^{m}$ satisfying (22) and (23), then by Lemma 2.1, we obtain that $b\in \frac{1}{{\beta }_{2}}\partial \phi \left({B}^{\alpha }u\right)$ and $-\frac{{\beta }_{2}}{{\beta }_{1}}{\left({B}^{\alpha }\right)}^{T}b\in \frac{1}{\alpha }\partial \varphi \left(u\right)$ . We can yield that

This implies that $u\in {R}^{m}$ is a solution of (21).

According to Proposition 2.1, we can conclude the following corollary.

Corollary 2.1. Suppose that ${u}_{0}\in {R}^{m}$ is given, $\lambda$ , $\mu$ are two positive numbers, ${B}^{\alpha }$ is a $2m\times m$ matrix, $\varphi$ is the function defined by (12), and $\phi$ is a differentiable convex function on ${R}^{2m}$ . If $u\in {R}^{m}$ is a solution of (21), then for any ${\beta }_{1}$ > 0,

Conversely, if for some ${\beta }_{1}$ > 0 there exists $u\in {R}^{m}$ satisfying (27), then $u\in {R}^{m}$ is a solution to (21).

Proof. By Proposition 2.1, a solution $u\in {R}^{m}$ of (21) satisfies (22) and (23). If the function $\phi$ is differentiable, then $\partial \phi \left(u\right) = \{\nabla \phi (u\left)\right\}$ , where $\partial \phi \left(u\right)$ is the gradient of $\phi$ at $u$ . Therefore, (6) and (23) imply that $b = \frac{1}{{\beta }_{2}}\partial \phi \left({B}^{\alpha }u\right)$ . Hence, (22) yields the fixed-point equation (27).

The fixed-point equation (27) can be viewed as an instance of the split forward-backward formula ^[31]. Suppose that $\partial \mathrm{\phi }$ is Lipschitz continuous with a Lipschitz constant $L$ , that is

and that ${\beta }_{1}$ is chosen to satisfy

It was proved in ^[33], that for any initial point ${u}^{0}\in {R}^{m}$ , the Picard iteration

converges to a fixed point of (27), which is a minimum of (21).

Let $Hu: = u-\frac{1}{{\beta }_{1}}({B}^{\alpha }{)}^{T}\partial \phi ({B}^{\alpha }{u}^{k})$ and $Qu: = (pro{x}_{\frac{1}{{\beta }_{1}}\varphi }\circ H)u$ . To prove that (30) is convergent, we only need to prove that $H$ and $Q$ are non-expansive averaged operators. We recall the definitions of non-expansive operators ^[31].

Definition 2.3. (Non-expansive operator): An operator $T$ on ${R}^{m}$ is non-expansive if it satisfies the following condition $\forall x, y\in {R}^{m}:{‖Tx-Ty‖}_{2}\le {‖x-y‖}_{2}$ .

Both $pro{x}_{f}\left(x\right)$ and $(I-pro{x}_{f})\left(x\right)$ are operators; see ^[31].

Definition 2.4. (Firmly non-expansive operator): An operator $T$ on ${R}^{m}$ is firmly non-expansive if it satisfies the following condition $\forall x, y\in {R}^{m}:{‖Tx-Ty‖}^{2}\le < x-y, Tx-Ty >$ .

Definition 2.5. (Non-expansive averaged operators): A non-expansive operator $Q$ on ${R}^{m}$ is a non-expansive averaged operator if there exists $k\in \left(\mathrm{0, 1}\right)$ and it satisfies the following condition $\forall x, y\in {R}^{m}:Q = kI+(1-k)P$ , where $P$ is a non-expansive operator. If $k = \frac{1}{2}$ , then $Q$ is a firmly non-expansive operator.

Both $pro{x}_{f}\left(x\right)$ and $(I-pro{x}_{f})\left(x\right)$ are firmly non-expansive operators; see ^[32].

Proposition 2.2. If $\phi$ is a convex function and ${B}^{\alpha }$ is a $2m\times m$ matrix, then $H$ is firmly non-expansive.

Proof. First, by the definition of the operator $H$ , $\forall x, y\in {R}^{m}$ , we have

We have

According to the sub-gradient inequalities of convex functions, we have

Substituting (35) into (33), we have

Combining (36) with (34), we have

We have

This completes the proof.

If $H:{R}^{m}\to {R}^{m}$ is firmly non-expansive, then $H$ is a non-expansive $\frac{1}{2}$ -averaged operator (see Lemma 3.8 in ^[26]). Thus $Q$ is a non-expansive averaged operator (see Lemma 3.7 in ^[26]).

We prove the convergence of (30). To simplify (30) and find an iterative format that is equivalent to (30), we make the following substitution

Let $v\in {R}^{m}$ be a given vector and $x\in {R}^{m}$ ; we denote the proximity operator of $\frac{1}{{\beta }_{1}}\varphi$ for the given $v\in {R}^{m}$ as follows

We have

Let $g$ and $f$ be two functions on ${R}^{m}$ ; then, we have

Because the function $g$ is differentiable, it can be expanded by applying the Taylor formula to $\left(v-{u}_{0}\right)\in {R}^{m}$ :

where $r$ denotes a constant greater than $1$ .

We can use (43) to find the following minimum value problem:

By (41), we can get

Using (45), we obtain

Therefore, substituting (42), (44) and (46) into (40), we conclude that

We can combine (38) and ${u}^{k+1} = pro{x}_{\frac{1}{{\beta }_{1}}\varphi }\left(v\right)$ with (47) to obtain

Substituting ${b}^{k} = \frac{1}{{\beta }_{2}}\partial \phi \left({B}^{\alpha }{u}^{k}\right)$ into (48) shows that (49) and (50) are equivalent iterations of (30).

Hence, according to the iterative equations (48) and (49), We can propose the following algorithm.

3.   Numerical results

This section describes several image denoising experiments that were conducted to demonstrate the behavior of the proposed algorithm. The peak signal to noise ratio (PSNR) is currently the most widely used tool for objectively evaluating image quality, and it is consistent with human subjective perception. A larger value of PSNR indicates better quality of the recovered image. It is defined as follows:

where ${u}^{\mathrm{*}}$ is the original image and $u$ is the denoised image. All experiments' iterations were ceased when the following criterion was satisfied:

In this study, images of size 256×256 pixels were used to conduct numerical experiments with $r = \frac{{\beta }_{1}}{{\beta }_{1}+2\lambda }$ . We used the ${L}^{1}{L}^{2}/T{V}^{\alpha }$ model to remove Gaussian noise, salt-and-pepper noise, and mixed noise. Original images of the experiment are shown in Figure 1. In particular, the different noise regimes yielded different results, as shown in Figure 2. Salt-and-pepper noise involves setting a value of a pixel to the minimal or maximal value of the image intensity range. Gaussian noise may extend this intensity range. We considered adding salt-and-pepper noise to the original image after Gaussian noise.

This study included a total of four groups of experiments. The first experiment was to restore images affected with $\sigma = 20$ , which is the level of Gaussian noise. The second experiment was to restore images affected with $s = 0.03$ , which is the level of salt-and-pepper noise. The third experiment was to restore images affected by the mixed noise. The fourth experiment was to explore the convergence of our proposed fractional-order TV denoising algorithm.

We began by investigating the effects of different parameters on the experimental results. Inspired by ^[27], we consistently chose $\alpha = 6, \beta = 128$ . We determined the most suitable values $\lambda$ and $\mathrm{\mu }$ through trial and error. When Gaussian noise with $\sigma = 20$ was added to the image 'Lena', we found that $\lambda = 0.07, \mu = 0$ performed better. When salt-and-pepper noise with $s = 0.03$ was added to the image 'Square', we found that $\lambda = 0, \mu = 3$ performed better. We verified that these selected parameters were effective for other images with the same noise levels. Additionally, we increased $\alpha$ from 0.8 to 1.9 in increments of 0.1.

In the first experiment, the Gaussian noise was added to the 'Lena' image at different levels. We chose $\lambda = 0.07, \mu = 0$ to deal with noisy images. Table 1 shows the values of PSNR, while Figure 3 shows the experimental results. In addition, Gaussian noise was added to the other images at $\sigma = 20$ . Table 2 shows the values of PSNR. The first experimental results demonstrated that $\alpha$ has an impact on the denoising results. The best denoising result often did not appear when $\alpha = 1$ . Therefore, the fractional-order TV model can be applied to improve the denoising performance of the TV model.

In the second experiment, we chose $\lambda = 0,$ $\mu = 4.8$ to deal with the 'Square' image corrupted by the salt-and-pepper noise at noise levels of 0.01, 0.02, 0.03, 0.05. Table 3 shows the values of PSNR, while Figure 4 shows the experimental results. In addition, salt-and-pepper noise was applied to the other images at $s = 0.03$ . Table 4 shows the values of PSNR. The experimental results indicate that when $\alpha$ is larger, the effect is better.

In the third experiment, we added Gaussian noise at $\sigma = 20$ and salt-and-pepper noise at $s = 0.03$ to four images and explore the performance of the algorithm. We chose $\lambda = 0.009$ , $\mu = 2.3$ . Table 5 shows the values of PSNR, while Figure 5 shows the experimental results Figure 6 shows the original image, the noisy image, and the denoised image for different values of $\alpha$ (from $0.8$ to $1.9$ ). The third and fourth rows represent their corresponding contour map. The data from Table 5 indicate that a larger $\alpha$ yields better denoising performance. Consequently, the fractional-order TV model outperformed the traditional TV model under mixed noise.

Based on the PSNR values and denoised images from the first three experiments, we can see that the fractional-order TV model can effectively reduce the block effect and perform better than the TV model.

The fourth experiment focused on the convergence of the algorithm. We applied $\alpha = 1.8$ as an example in the $\alpha$ range of 0.8 to 1.9. The PSNR value was recorded at each iteration. Figure 6 shows the experimental results. The blue line represents the noisy image results for $\sigma = 15, s = 0.01$ , the red line represents the noisy image results for $\sigma = 20, s = 0.03$ , and the yellow line represents the noisy image results for $\sigma = 20, s = 0.05$ . From Figure 6, it is obvious that our proposed fractional-order TV denoising algorithm is convergent.

Furthermore, we will show that our proposed model demonstrated good performance on the task of removing mixed noise. For this purpose, we added Gaussian noise at $\sigma = 20$ and salt-and-pepper noise at $s = 0.03$ to image 'Lena'. We chose $\alpha = 1.9$ . Figure 7 shows the 60th and 100th rows of the 'Lena' image from a one-dimensional perspective. The original, noisy and denoised images are represented by black, pink and blue lines, respectively. The blue solid line and the black solid line nearly coincide, which indicates that our proposed model exhibited good denoising performance. Figure 8 shows that the histogram for the noisy image was completely different from that of the original image, while the histogram for the denoised image was similar to the histogram for the original image. We took a small part of the 'Lena' image and marked it with a red rectangle; the experimental results can be seen in Figure 9.

4.   Conclusions

In this paper, we developed a fractional-order TV ( ${L}^{1}{L}^{2}/{TV}^{\alpha }$ ) model to remove mixtures of Gaussian noise and salt-and-pepper noise, by incorporating an ${L}^{1}$ -data fidelity term and ${L}^{2}$ -data fidelity term into the model. The existence of the solution of this model has been proved. We solved the proposed model by using the proximity algorithm, which prevents non-differentiability of the fractional order TV regularization terms. The convergence of the algorithm has been proved. The numerical experiments revealed the following: (1) The ${L}^{1}{L}^{2}/{TV}^{\alpha }$ model can effectively reduce the block effect and achieve better denoising performance than the ${L}^{1}{L}^{2}/TV$ model. (2) The ${L}^{1}{L}^{2}/{TV}^{\alpha }$ model effectively removes the mixture of Gaussian noise and salt-and-pepper noise owing to the proximity algorithm. (3) In the ${L}^{1}{L}^{2}/{TV}^{\alpha }$ model, $\alpha$ should range from 0.8 to 1.9. Different images will have different optimal values of $\alpha$ .

Author contributions

Donghong Zhao: Conceptualization, funding acquisition, supervision, methodology; Ruiying Huang: Writing-original draft, writing-review & editing, software, formal analysis; Li Feng: Writing-original draft, software.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this article.

Acknowledgments

This research was funded by Graduate Online Course Research Project of USTB (2024ZXB002), the National Natural Science Foundation of China (grant number 12371481) and Youth Teaching Talents Training Program of USTB (2016JXGGRC-002).

Conflict of interest

All authors declare no conflict of interest.