Poisson image restoration via an adaptive Euler’s elastica regularization

Abstract

Many recent studies have shown that Euler’s elastica regularization performs better than the famous total variation (TV) regularization on keeping image features in smooth regions during the process of denoising. In addition, an adaptive weighted matrix combined with total variation is a key technique for well restoring local features of image. Considering these two factors, in this paper, we propose an adaptive Euler’s elastica model for Poisson image restoration so as to well preserve both image features in smooth regions and local features of image. To solve this non-smooth and non-convex model efficiently, we design an alternating direction method of multipliers. Experiments on both natural and synthetic images illustrate the effectiveness and efficiency of the proposed method over the state-of-the-art methods in Poisson restoration and denoising, respectively. In particular, for Poisson restoration, our proposed method improves the TV method up to 2.46 about PSNR for dealing with the Peppers image with Gaussian blur and noise level σ = 1. In addition, the proposed method for Poisson denoising gets higher PSNR and SSIM values than the TAC method, while costing less CPU time.

Keywords

Poisson restoration Euler’s elastica regularization adaptive weighted matrix alternating direction method of multipliers (ADMM)

1 Introduction

Noise and blur are inevitable in the process of image formation, image transmission and image display. It is very necessary to remove them through image restoration technology. Image restoration has received extensive attention in many application fields, such as astronomical and medical imaging, computer science, electronics and engineering [17, 28]. Its goal is to obtain a clear image u from the degraded image f. It is viewed as an inverse problem. Due to the lack of some prior information, this inverse problem is ill-posed. Therefore, many regularization techniques [1] have been widely studied to overcome this disadvantage. The variational regularization problem for image restoration can be constructed as follows:

$min_{u} λ F (u, f) + R (u),$ (1.1) where $f \in ℝ^{M \times N}$ represents the degraded image, u stands for an image restored from f, λ > 0 is a regularization parameter to affect the balance between the fidelity term $F (\cdot, \cdot)$ and the regularization term $R (\cdot)$ . Typical examples of $R (u)$ include total variation (TV) regularization [27], high-order TV (HOTV) regularization [6 , 44], total generalized variation (TGV) regularization [4], curvature based regularization (including Euler’s elastica (EE) [6 , 39], Gaussian curvature [9] and mean curvature [10 , 48]) etc. Among them, TV regularization proposed by Rudin, Osher and Fatemi is one of the most famous methods, which has been broadly used in many other image applications, such as decomposition, segmentation, reconstruction and so on.

In numerous practical applications, the obtained images are more vulnerable to Poisson noise. Different from the Gaussian noise, Poisson noise depends on the image intensity, which makes Poisson image restoration very challenging. The restoration models proposed for Gaussian noise are not effective in Poisson noise removal [18].

To remove Poisson noise effectively, the TV-based regularization models were proposed in [23 , 45]. Owing to the desirable ability on edge-preserving of TV regularization, it has been widely used in image processing and computer vision [22 , 36]. However, the main disadvantage of TV regularization is that it often produces staircase effects in the flat regions of the restored images, especially for piece-wise smooth images. In order to overcome the staircase effects, many improved models were proposed for Poisson image restoration. Specifically, in [14 , 43], some HOTV-based regularizations were introduced. Compared with the TV regularization, these ones have better ability to maintain features in smooth regions. However, they usually result in blurring edges of restored images. To preserve edges while suppressing the staircase effects, a Poisson restoration model based on TGV regularization was proposed in [11]. Since the TGV regularization has the ability to construct the piece-wise polynomial function, the proposed model can describe intensity variations in smooth regions more precisely. On this basis, a directional TGV regularization [8] was proposed for Poisson image restoration. In [19], the authors integrated an anisotropic diffusion tensor into the TGV regularization and thus proposed an improved model for Poisson noise removal. This combination yielded a significant quality improvement visually and numerically. But the introduced diffusion tensor increases the computational burden, which needs more CPU time to optimize model. In terms of both preserving edges and overcoming the staircase effects, a non-convex EE model was presented for Gaussian noise removal [31]. This model can faster suppress the high frequency components of images. Thus, it can effectively alleviate the staircase effects and make the denoised images more natural. Due to its non-convexity, non-smoothness and nonlinearity, the solution of the EE model is a highly challenging task. To overcome these problems and reduce computational cost, several fast algorithms have been studied [13 , 46]. Inspired by the advantages of EE regularization, we consider it for Poisson image restoration.

Generally, an image contains many different features. For the traditional TV-based models, the two subvariables in the gradient operator use the same weight. In this case, these models cannot effectively couple with local features of the image [12 , 26]. Therefore, we hope that the proposed model can diffuse along the tangent direction of local features. To this end, Pang et al. [26] proposed an anisotropic TV (ATV) based on local information, which adds different weight to each component of the gradient operator. It can better describe the local features of image than other TV-based models. Inspired by this model, Wang et al. [33] constructed an anisotropic mean curvature model. This model combined the gradient operator with the adaptive weighted matrix in the mean curvature regularization term, and could describe the local features in image efficiently. In [20], the authors introduced a weighted anisotropic total variation (WATV) model for Gaussian image restoration, where a nonlinear and monotonically increasing function was used as the weight. In some special images, the texture details in images show obvious directionality. Thus, the designed model needs to be able to describe this geometric feature. The directional total variation (DTV) method which rotates the gradient operator by adding a rotation matrix can effectively couple with local structures of image [3 , 40]. However, while processing some special images with several dominant directions, the ideal restoration effect will not be achieved by using previous DTV-based models. The reason is that one fixed angle parameter in the rotation matrix cannot simultaneously describe multiple dominant directions [25]. Considering this fact, by combining a rotation matrix, an adaptive weighted matrix and the TV^p-quasinorm, the authors proposed an adaptive weighted TV^p (AWTV^p) regularization model in [25]. The l^p-quasinorm was used to promote the sparsity by setting p ∈ (0, 1). Specifically, the adaptive weighted matrix can enhance the diffusion along the tangential direction of the edge. Combining it with the rotation matrix, the proposed model can handle images with complex structures, which have several dominant directions. Motivated by [25], in our recent work, an adaptive weighted directional TV^p (AWDTV^p) regularization was proposed for Poisson image restoration [42].

Different from traditional TV-based models, the adaptive TV-based models can obtain better restoration effect owing to the use of some adaptive weighted matrices. Specifically, they can effectively preserve image details. However, it still exists some problems while restoring smooth images. Considering the superior performance of EE model on overcoming the staircase effects and maintaining the structural information [30, 41], we designed an adaptive Euler’s elastica regularization model for Gaussian denoising [37]. The experimental results demonstrate that the proposed method can obtain the highest PSNR and SSIM values in most cases. In addition, it can improve the image contrast and suppress the staircase effects. Encouraged by the successful application of this regularization, in this paper, we construct a new adaptive Euler’s elastica regularization for Poisson image restoration.

As for tackling the optimization problem, several efficient numerical algorithms have been explored, such as the primal dual algorithm [35], the split Bregman method [29], the nonlinear multigrid method [44] and the alternating direction method of multipliers (ADMM) [7, 34]. In this paper, we adopt the ADMM which has been applied in numerous applications. The main contributions of this paper can be summarized as follows:

The models related to Euler’s elastica energy are rarely used in Poisson image restoration. By integrating a new adaptive weighted matrix into the Euler’s elastica regularization, we develop an adaptive Euler’s elastica model with better adaptability and more stronger restoration ability. This combination can effectively overcome the staircase effects and preserve image details, even for some smooth images.

Extensive experiments are carried out on natural and synthetic images to show the efficiency of our proposed model and algorithm. For Poisson image restoration, we compare with the TV model and the AWDTV^p model, while the comparison of the model and the TAC model [46] is conducted for Poisson denoising.

This paper is organized as follows. Section 2 mainly gives a brief introduction of some previous related models. In Section 3, an adaptive Euler’s elastica model for Poisson image restoration is proposed. Section 4 presents an efficient ADMM for solving this model. Section 5 shows some experimental results on natural images and synthetic images for Poisson image restoration and denoising respectively, which demonstrate the effectiveness and efficiency of the proposed method. Finally, we draw some conclusions in Section 6.

2 Related models

In this section, we mainly introduce two previous related models for Poisson image restoration. Experiments on these two models will be carried out and compared with the experimental results of our proposed model.

2.1 TV model

As is well-known, the TV regularization was first proposed by Rudin et al. in [27] for removing Gaussian noise, which has been used in many fields such as image restoration, image inpainting, medical image reconstruction and so on. Particularly, the TV-based Poisson restoration model [23 , 45] can be described as

$min_{u} ∥ \nabla u ∥_{1} + λ 〈 Ku - f log Ku, 1_{Ω} 〉,$ (2.2) where λ > 0 is a regularization parameter, and 1_Ω denotes the indicator function on a discrete grid of the image domain Ω. As one of the most popular regularizations, TV performs very well on dealing with piece-wise constant images. Namely, it can better preserve image edges in the process of restoration. However, it may cause staircase effects, especially for piece-wise smooth images. For all that, it has been widely used in various image processing problems owing to its advantages.

2.2 AWDTV^p model

In our recent work [42], an adaptive weighted directional TV^p regularization was designed for Poisson restoration. The AWDTV^p model is formulated as

$min_{u} ∥ T_{β} (f) R_{- θ} \nabla u ∥_{p}^{p} + λ 〈 Ku - f log Ku, 1_{Ω} 〉,$ (2.3) where T_β (f) and R_-θ stand for the adaptive weighted matrix and the rotation matrix respectively, which are defined by $\begin{matrix} T_{β} (f) & = [\begin{matrix} t_{1, β} (f) & 0 \\ 0 & t_{2} (f) \end{matrix}] \\ = [\begin{matrix} \frac{β (f)}{1 + ι | G_{\hat{δ}} * \nabla_{x} f |^{2}} & 0 \\ 0 & \frac{1}{1 + ι | G_{\hat{δ}} * \nabla_{y} f |^{2}} \end{matrix}], \end{matrix}$ and $R_{- θ} = [\begin{matrix} cos θ & sin θ \\ - sin θ & cos θ \end{matrix}] .$ Here $G_{\hat{δ}} (\cdot)$ represents the Gaussian convolution kernel, ι is a tuning parameter, β (f) is used to enhance the diffusion, and θ is the affine angle. Particularly, when $∥ G_{\hat{δ}} * \nabla f ∥_{2} > c$ , we set β (f) = β. On the contrary, we set β (f) =1.

For handling images which have obvious directionality, it is very useful to add a rotation matrix to rotate the gradient operator. When the image has several dominant directions, it needs to further combine the adaptive weighted matrix with the rotation matrix so as to get better recovery results. Besides, the l^p-quasinorm used in model (2.2) is to promote gradient sparsity.

3 The proposed model

Although the TV-based models can effectively preserve image edges, it fails to restore smooth images. It often causes staircase effects in flat regions. To overcome this problem, some high-order regularization models were proposed, including the popular curvature regularization. Inspired by the successful application of curvature regularization, in this section, we construct a new Poisson restoration model by adding an adaptive weighted matrix to the Euler’s elastica regularization.

In order to obtain restored images with high quality, the following model based on an adaptive Euler’s elastica regularization is proposed for Poisson restoration:

$min_{u} \sum_{i, j} g (κ (u_{i, j})) | T \nabla u_{i, j} | + λ 〈 Ku - f log Ku, 1_{Ω} 〉 .$ (3.1) Here, g (κ (u_i,j)) is a certain function of the curvature defined as follows:

$g (κ (u_{i, j})) = 1 + α | κ (u_{i, j}) |,$ (3.2) where α is a positive constant. For a 2D curve, the curvature κ can be expressed as a function of u. Namely,

$κ (u_{i, j}) = \nabla \cdot (\frac{\nabla u_{i, j}}{| \nabla u_{i, j} |}) .$ (3.3)T is an adaptive weighted matrix defined as

$\begin{matrix} T (i, j) = [\begin{matrix} t_{1} (i, j) & 0 \\ 0 & t_{2} (i, j) \end{matrix}] \\ = [\begin{matrix} \frac{β (f)}{1 + ι | G_{\hat{δ}} (i, j) * \nabla_{x} f (i, j) |^{2}} & 0 \\ 0 & \frac{1}{1 + ι | G_{\hat{δ}} (i, j) * \nabla_{y} f (i, j) |^{2}} \end{matrix}] . \end{matrix}$ (3.4) To effectively describe diffusion, β (f) is used to represent the smoothing and edge regions. If ι = 0 and β (f) =1, T will become the identity matrix $I$ . In this case, the proposed model is related to the Euler’s elastica model. In addition, if α = 0 and β (f) =1, the proposed model is also related to the ATV model proposed in [26].

4 The proposed algorithm

Since the minimization problem (3.1) is non-convex and non-smooth, we will face some numerical difficulties while solving it. To deal with this problem, we design an efficient alternating direction method of multipliers (ADMM), which can decompose a large global problem into a series of smaller local subproblems and then use local solutions to calculate the solution of the original problem [2, 5].

We first introduce three auxiliary variables v, w and q, and transform the original unconstrained optimization problem (3.1) into a constrained optimization problem as follows

$\begin{matrix} min_{u, q, v, w} \sum_{i, j} g (κ (u_{i, j})) | w_{i, j} | + λ 〈 q - f log q, 1_{Ω} 〉, \\ s . t . v = \nabla u, w = Tv, q = Ku . \end{matrix}$ (4.1)

To solve the constrained problem (4.1), we introduce three Lagrangian multipliers η₁ = (η₁₁, η₁₂) ^T, η₂ = (η₂₁, η₂₂) ^T and η₃, and then reformulate it as a saddle-point problem. The corresponding augmented Lagrangian functional is

$\begin{matrix} L (u, q, v, w; η_{1}, η_{2}, η_{3}) \\ = \sum_{i, j} g (κ (u_{i, j}^{n})) | w_{i, j} | + λ 〈 q - f log q, 1_{Ω} 〉 \\ + 〈 η_{1}, v - \nabla u 〉 + \frac{γ_{1}}{2} ∥ v - \nabla u ∥_{2}^{2} \\ + 〈 η_{2}, w - Tv 〉 + \frac{γ_{2}}{2} ∥ w - Tv ∥_{2}^{2} \\ + 〈 η_{3}, q - Ku 〉 + \frac{γ_{3}}{2} ∥ q - Ku ∥_{2}^{2}, \end{matrix}$ (4.2) where γ₁, γ₂, γ₃ > 0 are three penalty parameters. In (4.1), it includes four subvariables, so we iteratively and alternatively solve each subvariable when using the ADMM until reaching the terminating condition; see Algorithm 1.

Algorithm 1

ADMM for solving (3.1)

1: Input: a blurred noisy image f

2: Parameters: λ, ι, δ, β, α, nMax, and ɛ ∈ mathdsR

3: Initialize: u⁰ = f, q⁰, v⁰, w⁰, $η_{1}^{0}$ , $η_{2}^{0}$ , $η_{3}^{0}$ , and n = 0

4: while n< nMax do ${\begin{matrix} u^{n + 1} = arg min_{u} L (u, q^{n}, v^{n}, w^{n}; η_{1}^{n}, η_{2}^{n}, η_{3}^{n}), (4.3 a) \\ q^{n + 1} = arg min_{q} L (u^{n + 1}, q, v^{n}, w^{n}; η_{1}^{n}, η_{2}^{n}, η_{3}^{n}), (4.3 b) \\ v^{n + 1} = arg min_{v} L (u^{n + 1}, q^{n + 1}, v, w^{n}; η_{1}^{n}, η_{2}^{n}, η_{3}^{n}), (4.3 c) \\ w^{n + 1} = arg min_{w} L (u^{n + 1}, q^{n + 1}, v^{n + 1}, w; η_{1}^{n}, η_{2}^{n}, η_{3}^{n}), (4.3 d) \\ η_{1}^{n + 1} = η_{1}^{n} + γ_{1} (v^{n + 1} - \nabla u^{n + 1}), (4.3 e) \\ η_{2}^{n + 1} = η_{2}^{n} + γ_{2} (w^{n + 1} - T v^{n + 1}), (4.3 f) \\ η_{3}^{n + 1} = η_{3}^{n} + γ_{3} (q^{n + 1} - {Ku}^{n + 1}) . (4.3 g) \end{matrix}$

5: Assign n by n + 1

6: If ∥uⁿ - u^n-1 ∥ ₂/∥ uⁿ ∥ ₂ ≤ ɛ, stop the iteration

7: end while

8: return u^∗ = uⁿ as the restored image

In Algorithm 1, the efficiency for solving the optimization problem relies on how to solve the subproblems with high accuracy. The above subproblems all have explicit solutions.

Subproblem (4.3a). This subproblem can be reformulated as $\begin{matrix} u^{n + 1} = arg min_{u} 〈 η_{1}^{n}, v^{n} - \nabla u 〉 \\ + \frac{γ_{1}}{2} ∥ v^{n} - \nabla u ∥_{2}^{2} + 〈 η_{3}^{n}, q^{n} - Ku 〉 + \frac{γ_{3}}{2} \\ ∥ q^{n} - Ku ∥_{2}^{2} . \end{matrix}$ (4.4) According to its optimality conditions, we easily derive the Euler-Lagrange equation of (4.3a): $\begin{matrix} (γ_{3} K^{*} K - γ_{1} div \nabla) u^{n + 1} \\ = K^{*} η_{3}^{n} + γ_{3} K^{*} q^{n} - div η_{1}^{n} - γ_{1} div v^{n} . \end{matrix}$ (4.5)

Under periodic boundary condition, the equation (4.5) can be solved efficiently by the fast Fourier transform (FFT) in the following way: $\begin{matrix} u^{n + 1} = F^{- 1} \\ (\frac{F (K^{*} η_{3}^{n} + γ_{3} K^{*} q^{n} - div η_{1}^{n} - γ_{1} div v^{n})}{γ_{3} F (K^{*} K) - γ_{1} F (Δ)}) . \end{matrix}$ (4.6)

Subproblem (4.3b). This subproblem can be expressed as $\begin{matrix} q^{n + 1} = arg min_{q} λ 〈 q - f log q, 1_{Ω} 〉 \\ + 〈 η_{3}^{n}, q - {Ku}^{n + 1} 〉 + \frac{γ_{3}}{2} ∥ q - {Ku}^{n + 1} ∥_{2}^{2}, \end{matrix}$ (4.7) which has a closed form solution. By solving a quadratic equation, the solution can be obtained by the following formula $\begin{matrix} q^{n + 1} = \frac{{Ku}^{n + 1} - \frac{λ + η_{3}^{n}}{γ_{3}}}{2} \\ + \sqrt{{(\frac{{Ku}^{n + 1} - \frac{λ + η_{3}^{n}}{γ_{3}}}{2})}^{2} + \frac{λ f}{γ_{3}}} . \end{matrix}$ (4.8)

Subproblem (4.3c). The v-subproblem can be written as $\begin{matrix} v^{n + 1} & = arg min_{v} 〈 η_{1}^{n}, v - \nabla u^{n + 1} 〉 \\ + \frac{γ_{1}}{2} ∥ v - \nabla u^{n + 1} ∥_{2}^{2} + 〈 η_{2}^{n}, w^{n} - T v 〉 \\ + \frac{γ_{2}}{2} ∥ w^{n} - Tv ∥_{2}^{2} . \end{matrix}$ (4.9) Considering its optimality conditions, we have $\begin{matrix} (γ_{1} I + γ_{2} T^{2}) v^{n + 1} \\ = - η_{1}^{n} + γ_{1} \nabla u^{n + 1} + T η_{2}^{n} + γ_{2} T w^{n}, \end{matrix}$ (4.10) where $I$ is a unit operator. Denote v = [v₁, v₂] ^T, w = [w₁, w₂] ^T, η_i = [η_i1, η_i2] ^T (i = 1, 2), then vⁿ⁺¹ satisfies the following linear equation system: $\begin{matrix} [\begin{matrix} γ_{1} + γ_{2} t_{1}^{2} & 0 \\ 0 & γ_{1} + γ_{2} t_{2}^{2} \end{matrix}] [\begin{matrix} v_{1}^{n + 1} \\ v_{2}^{n + 1} \end{matrix}] \\ = [\begin{matrix} - η_{11}^{n} + γ_{1} \nabla_{x} u^{n + 1} + t_{1} η_{21}^{n} + γ_{2} t_{1} w_{1}^{n} \\ - η_{12}^{n} + γ_{1} \nabla_{y} u^{n + 1} + t_{2} η_{22}^{n} + γ_{2} t_{2} w_{2}^{n} \end{matrix}] . \end{matrix}$ (4.11) Thus, the explicit solution of vⁿ⁺¹ can be obtained by simple calculation, ${\begin{matrix} v_{1}^{n + 1} = \frac{- η_{11}^{n} + γ_{1} \nabla_{x} u^{n + 1} + t_{1} η_{21}^{n} + γ_{2} t_{1} w_{1}^{n}}{γ_{1} + γ_{2} t_{1}^{2}}, \\ v_{2}^{n + 1} = \frac{- η_{12}^{n} + γ_{1} \nabla_{y} u^{n + 1} + t_{2} η_{22}^{n} + γ_{2} t_{2} w_{2}^{n}}{γ_{1} + γ_{2} t_{2}^{2}} . \end{matrix}$ (4.12)

Subproblem (4.3d). The w-subproblem can be written as $\begin{matrix} w^{n + 1} = arg min_{w} \sum_{i, j} g (κ (u_{i, j}^{n + 1})) | w_{i, j} | \\ + \frac{γ_{2}}{2} {∥ w - (T v^{n + 1} - \frac{η_{2}^{n}}{γ_{2}}) ∥}_{2}^{2} . \end{matrix}$ (4.13) Then we can get its closed form solution by the soft-thresholding operator

$\begin{matrix} w_{i, j}^{n + 1} & = soft ((T v^{n + 1} - \frac{η_{2}^{n}}{γ_{2}})_{i, j}, \frac{g (κ (u_{i, j}^{n + 1}))}{γ_{2}}) \\ = \max (| (T v^{n + 1} - \frac{η_{2}^{n}}{γ_{2}})_{i, j} | - \frac{g (κ (u_{i, j}^{n + 1}))}{γ_{2}}, 0) \\ \cdot sign ((T v^{n + 1} - \frac{η_{2}^{n}}{γ_{2}})_{i, j}) . \end{matrix}$ (4.14)

5 Numerical experiments

In this section, the proposed model for Poisson image restoration called AEEPR is compared with TV model [29] and AWDTV^p model [42]. Especially, two types of blur and two noise levels are tested to show the superiority of the proposed model. In addition, we apply the proposed method to Poisson image denoising and compare it with the previous curvature regularization model proposed recently in [46]. All the numerical experiments are conducted via MATLAB on a windows 10 (64bit) desktop computer with 2.90 GHz Intel(R) Core(TM) i7-10700 CPU and 32 GB RAM.

To assess the quality of restored images, two evaluation indexes are utilized, i.e. the peak signal-to-noise ratio (PSNR) and the structural similarity (SSIM). They are defined as follows:

$PSNR = 10 {log}_{10} \frac{{Max}^{2}}{\frac{1}{MN} \sum_{i = 1}^{M} \sum_{j = 1}^{N} (u_{i, j} - I_{i, j})^{2}},$ (5.1)

$SSIM = \frac{(2 μ_{I} μ_{u} + c_{1}) (2 σ_{Iu} + c_{2})}{({μ_{I}}^{2} + {μ_{u}}^{2} + c_{1}) ({σ_{I}}^{2} + {σ_{u}}^{2} + c_{2})},$ (5.2) where I_i,j and u_i,j denote the pixel values of the clean image and the recovery image, respectively. And Max stands for the maximum pixel value of the clean image I. μ_I and μ_u denote the local mean values of images I and u. σ_I and σ_u represent the respective standard deviations, and σ_Iu is the covariance value between the clean image I and the recovery image u. c₁ and c₂ are two constants used to avoid near-zero denominator values. Generally speaking, a larger PSNR value indicates the noise after denoising is less, while a larger SSIM value reflects higher similarity between the restored image u and the clean image I. The definition of the relative error in uⁿ is given as

$R (u^{n}) = \frac{∥ u^{n} - u^{n - 1} ∥_{2}}{∥ u^{n} ∥_{2}} .$ (5.3) In all experiments, the iteration will terminate when the relative error R (uⁿ) ≤10^-4 or the iteration reaches 500.

Six test images including natural images (“Peppers” (256 × 256), “Lena” (256 × 256), “Barbara” (512 × 512)) and synthetic images (“Dot” (256 × 256), “Shape” (128 × 128), “Triangle” (254 × 214)) are selected in experiments, as shown in Fig. 1.

Fig. 1

Test images.

5.1 Poisson restoration

For Poisson image restoration, there are some parameters in the proposed model and algorithm. For all parameters, we adopted a trial-and-error way to gradually adjust so as to determine the optimal parameter set that yielded the highest PSNR value. During the experiments, we first set several values of λ into a bigger range as [a, b]. Then we find the related λ value through the obtained maximum PSNR value to determine a suitable subset [a₁, b₁] ⊂ [a, b]. By repeating this strategy many times, we can find a more appropriate parameter value in subset [a_i, b_i]. When the difference between successive PSNR values is smaller than 0.01, we set this corresponding parameter value as the selected value of λ. And in turn, we use the same strategy for the choice of other parameters until the PSNR value almost unchanges during the cycle. The same strategy is employed to tune the parameters of the other two methods. During the parameters adjustment, the algorithms are terminated when the stopping criterion is satisfied. It is worth noting that in the process of parameters adjustment, we found that only a few parameters significantly affect the results. Generally, restoration results mainly depend on the regularization parameter λ, which is used to control the balance between the fidelity term and the regularization term. Besides, the main function of ι is to control the local adaptivity, δ is the standard deviation, and the weighted parameter β can enhance the diffusion. The positive parameter α can balance the effect of the curvature term on the length.

We compare the proposed AEEPR with TV and AWDTV^p for Poisson image restoration. The matlab functions $fspecial (^{'} {Gaussian}^{'}, [5, 5], \sqrt{3})$ and fspecial (′motion′, 10, 45) are chosen in the experiments to generate Gaussian blurring kernel and Motion blurring kernel, respectively. In addition, to test the effectiveness of the proposed method under different noise levels, the matlab routine poissrnd (B/σ) ∗ σ is applied to add Poisson noise. Here, B represents the blurry image degraded from the original clean image I. In the experiments, we test two noise levels σ = 1 and 5.

The optimal values of parameters for different restoration models under different blur kernels and noise levels are listed in Table 1. The number of parameters used in the AEEPR model (3.1) is less than that of AWDTV^p model (2.2). In these two models, the optimal values of the same parameters are taken in the same range. Table 2 records the PSNR and SSIM values of three methods for two natural and two synthetic images, and the corresponding CPU time and number of iterations are displayed in Table 3.

Table 1
Selected parameter sets for different restoration models under different blur kernels and noise levels

Image Blur Noise TV AWDTV^p AEEPR

σ λ (λ, ι, δ, β, p, θ) (λ, ι, δ, β, α)

Peppers Gaussian 1 65 (0.47,0.7,11,0,0.9,0) (0.5,0.7,7,1,0.3)

(256 × 256) 5 17 (0.1,0.6,10,1,0.9,0) (0.115,0.59,10,1,0.7)

Motion 1 70 (0.38,0.6,3,0,0.9,π/800) (0.8,0.67,5,1,0.0001)

5 17 (0.08,0.65,5,0,0.9,π/800) (0.19,0.71,5,1,0.01)

Lena Gaussian 1 53 (27,0.0007,0.3,2,0.99,π/32) (42,0.0006,0.7,2,0.001)

(256 × 256) 5 15 (5,0.0004,0.8,1,0.98,π/50) (10.2,0.0003,0.7,0,0.001)

Motion 1 54 (18,0.0004,0.5,1,0.99,0) (35,0.0004,0.4,0,0.008)

5 16 (6,0.0004,0.6,1,0.99,0) (11,0.0004,0.7,0,0.01)

Dot Gaussian 1 110 (12,0.005,0.1,1,0.85,0) (28,0.007,0.4,1,0.01)

(256 × 256) 5 25 (4.2,0.003,1.1,1,0.8,0) (12,0.004,1.1,0,0.02)

Motion 1 140 (42,0.0009,0.7,0,0.76,0) (130,0.0007,0.7,0,0.09)

5 35 (10,0.002,1.4,0,0.78,0) (15,0.002,1.6,0,0.1)

Shape Gaussian 1 170 (37,0.0006,1.2,2,0.53,π/512) (38,0.0007,1.3,2,0.002)

(128 × 128) 5 32 (6,0.002,2.6,10,0.53,π/512) (9,0.002,1.3,10,0.006)

Motion 1 160 (34,0.00005,4.2,0,0.53,0) (33,0.00006,0.6,1,0.02)

5 47 (9,0.00003,1.7,0,0.54,0) (10,0.00001,0.6,1,0.03)

Image	Blur	Noise	TV	AWDTV^p	AEEPR
		σ	λ	(λ, ι, δ, β, p, θ)	(λ, ι, δ, β, α)
Peppers	Gaussian	1	65	(0.47,0.7,11,0,0.9,0)	(0.5,0.7,7,1,0.3)
(256 × 256)		5	17	(0.1,0.6,10,1,0.9,0)	(0.115,0.59,10,1,0.7)
	Motion	1	70	(0.38,0.6,3,0,0.9,π/800)	(0.8,0.67,5,1,0.0001)
		5	17	(0.08,0.65,5,0,0.9,π/800)	(0.19,0.71,5,1,0.01)
Lena	Gaussian	1	53	(27,0.0007,0.3,2,0.99,π/32)	(42,0.0006,0.7,2,0.001)
(256 × 256)		5	15	(5,0.0004,0.8,1,0.98,π/50)	(10.2,0.0003,0.7,0,0.001)
	Motion	1	54	(18,0.0004,0.5,1,0.99,0)	(35,0.0004,0.4,0,0.008)
		5	16	(6,0.0004,0.6,1,0.99,0)	(11,0.0004,0.7,0,0.01)
Dot	Gaussian	1	110	(12,0.005,0.1,1,0.85,0)	(28,0.007,0.4,1,0.01)
(256 × 256)		5	25	(4.2,0.003,1.1,1,0.8,0)	(12,0.004,1.1,0,0.02)
	Motion	1	140	(42,0.0009,0.7,0,0.76,0)	(130,0.0007,0.7,0,0.09)
		5	35	(10,0.002,1.4,0,0.78,0)	(15,0.002,1.6,0,0.1)
Shape	Gaussian	1	170	(37,0.0006,1.2,2,0.53,π/512)	(38,0.0007,1.3,2,0.002)
(128 × 128)		5	32	(6,0.002,2.6,10,0.53,π/512)	(9,0.002,1.3,10,0.006)
	Motion	1	160	(34,0.00005,4.2,0,0.53,0)	(33,0.00006,0.6,1,0.02)
		5	47	(9,0.00003,1.7,0,0.54,0)	(10,0.00001,0.6,1,0.03)

Table 2

PSNR and SSIM obtained by different restoration models

Image			Peppers (256 × 256)		Lena (256 × 256)		Dot (256 × 256)		Shape (128 × 128)
Blur	Noise	Model	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
Gaussian	1	TV	25.568	0.899	26.812	0.880	22.289	0.871	17.270	0.745
		AWDTV^p	27.996	0.909	27.101	0.886	22.519	0.918	18.400	0.826
		AEEPR	28.023	0.910	27.063	0.885	22.535	0.920	18.554	0.831
	5	TV	24.413	0.862	25.163	0.833	20.691	0.852	16.677	0.732
		AWDTV^p	25.867	0.869	25.431	0.841	21.154	0.894	17.283	0.805
		AEEPR	25.895	0.869	25.412	0.839	21.177	0.893	17.334	0.806
Motion	1	TV	24.660	0.864	25.976	0.853	20.812	0.797	17.614	0.701
		AWDTV^p	26.388	0.870	26.389	0.862	20.870	0.864	18.803	0.776
		AEEPR	26.392	0.870	26.360	0.861	20.924	0.808	18.816	0.804
	5	TV	23.411	0.827	24.204	0.797	17.830	0.731	15.527	0.594
		AWDTV^p	24.703	0.828	24.646	0.807	17.928	0.793	16.417	0.704
		AEEPR	24.726	0.828	24.617	0.806	17.955	0.787	16.464	0.715

Table 3

CPU time and number of iterations for different restoration models

Image			Peppers (256 × 256)		Lena (256 × 256)		Dot (256 × 256)		Shape (128 × 128)
Blur	Noise	Model	CPU(s)	Iter	CPU(s)	Iter	CPU(s)	Iter	CPU(s)	Iter
Gaussian	1	TV	8.67	192	10.42	222	8.64	203	2.72	195
		AWDTV^p	24.06	352	12.31	194	18.39	274	6.81	348
		AEEPR	20.53	377	11.46	214	11.46	211	5.62	331
	5	TV	11.98	248	8.56	204	13.18	293	1.65	110
		AWDTV^p	35.13	500	22.03	311	20.25	281	5.38	254
		AEEPR	28.51	500	14.70	262	9.32	168	4.68	277
Motion	1	TV	14.28	322	9.56	217	22.19	500	7.00	500
		AWDTV^p	30.86	445	18.51	250	34.03	500	9.83	500
		AEEPR	26.64	500	13.78	249	28.60	500	8.31	500
	5	TV	20.39	436	9.26	179	15.32	357	7.30	500
		AWDTV^p	34.59	500	16.07	228	35.09	500	10.34	500
		AEEPR	29.62	500	11.42	213	27.59	500	8.35	500

From Table 2, we observe that AEEPR and AWDTV^p obtain almost the same PSNR and SSIM values, which are much higher than TV. For instance, AEEPR improves the TV up to 2.46 about PSNR for dealing with the Peppers image with Gaussian blur and noise level σ = 1. In addition, Table 3 shows that our AEEPR costs less CPU time than AWDTV^p in getting almost the same PSNR and SSIM values, especially in the case of higher noise level with the same blur.

To assess the visual quality of restored images obtained by different methods, we present the restored images and the locally enlarged images for two natural images in Figs. 3, while showing restored images and residual images f - u for two synthetic images in Figs. 5. From these figures, we can see that all methods remove most of the noise and obtain good recovery effects. However, the TV method produces obvious staircase effects in the flat regions, especially for the “Peppers" and “Lena" images. Obviously, the AEEPR and the AWDTV^p methods effectively overcome the staircase effects while retaining a lot of image details, such as neat edges and smoother regions. In particular, we observe this from a close-up of the nose and the head regions of “Lena" image. Figures 5 show that AEEPR and AWDTV^p remove more blur and noise than TV. On the whole, the AEEPR has better performance on restoring degraded images and costs less CPU time.

Fig. 2

Restored “Peppers” images and locally enlarged images obtained by different restoration models under Gaussian blur and noise level σ = 1.

Fig. 3

Restored “Lena” images and locally enlarged images obtained by different restoration models under motion blur and noise level σ = 5.

Fig. 4

Restored “Dot” images and their corresponding residual images generated by different restoration models under Gaussian blur and noise level σ = 5.

Fig. 5

Restored “Shape” images and their corresponding residual images generated by different restoration models under motion blur and noise level σ = 5.

At last, we empirically demonstrate the convergence of our proposed algorithm. We consider the case of Gaussian blur and noise level σ = 1, and the one of motion blur and noise level σ = 5. Furthermore, all four images are tested by Algorithm 4.1. From Figs. 7, all the plots of PSNR increase rapidly with respect to the CPU time while their plots of energy and relative error decrease accordingly. Thus, we can also draw a conclusion that our AEEPR is convergent.

Fig. 6

The plots of PSNR, relative error and energy versus CPU time for our AEEPR model under Gaussian blur and noise level σ = 1.

Fig. 7

The plots of PSNR, relative error and energy versus CPU time for our AEEPR model under motion blur and noise level σ = 5.

5.2 Poisson denoising

In this subsection, some Poisson denoising examples are given to further illustrate the advantages of our proposed regularization method. Especially, we compare with the TAC method proposed recently in [46]. The specific implementation is as follows.

5.2.1 Adaptive Euler’s elastica model

When K = I, the AEEPR model (3.1) is reduced to a Poisson denoising model called AEEPD, which can be written as follows:

$min_{u} \sum_{i, j} g (κ (u_{i, j})) | T \nabla u_{i, j} | + λ 〈 u - f log u, 1_{Ω} 〉 .$ (5.4) Here f is the noisy image polluted by Poisson noise. g (κ (u_i,j)) is a certain function of the curvature and T is an adaptive weighted matrix, which have been defined in (3.2) and (3.4) respectively.

To solve this minimization problem, we introduce two auxiliary variables v, w to rewrite it as the following constrained minimization problem:

$\begin{matrix} min_{u, v, w} \sum_{i, j} g (κ (u_{i, j})) | w_{i, j} | + λ 〈 u - f log u, 1_{Ω} 〉, \\ s . t . v = \nabla u, w = Tv . \end{matrix}$ (5.5) Then two Lagrangian multipliers η₁, η₂ are introduced to reformulate (5.5) as a saddle point problem, whose augmented Lagrangian functional is

$\begin{matrix} L (u, v, w; η_{1}, η_{2}) \\ = \sum_{i, j} g (κ (u_{i, j}^{n})) | w_{i, j} | + λ 〈 u - f log u, 1_{Ω} 〉 \\ + 〈 η_{1}, v - \nabla u 〉 + \frac{γ_{1}}{2} ∥ v - \nabla u ∥_{2}^{2} \\ + 〈 η_{2}, w - Tv 〉 + \frac{γ_{2}}{2} ∥ w - Tv ∥_{2}^{2} . \end{matrix}$ (5.6) While applying the ADMM, it needs to minimize three subproblems and update two Lagrangian multipliers with other variables fixed at each iteration. It should be pointed out that, we utilize a fixed-point strategy to solve the u-subproblem. More details of the ADMM for solving the above minimization problem (5.5) can be referred to [37, 46].

5.2.2 The previous curvature-based model

As a comparison, the variational model based on curvature regularization proposed in [46] is described as $min_{u} \sum_{i, j} g (κ (u_{i, j})) \sqrt{1 + | \nabla u_{i, j} |^{2}} + λ 〈 u - f log u, 1_{Ω} 〉 .$ (5.7)

Specifically, two auxiliary variables v and q are used to rewrite (5.7) as a constrained optimization problem,

$\begin{matrix} min_{v, q, u} \sum_{i, j} g (κ (u_{i, j})) \sqrt{1 + | v_{i, j} |^{2}} \\ + λ 〈 q - f log q, 1_{Ω} 〉, \\ s . t . v_{i, j} = \nabla u_{i, j}, q_{i, j} = u_{i, j}, \end{matrix}$ (5.8) where g (κ (u_i,j)) has the same definition as in (3.2). The corresponding augmented Lagrangian functional of (5.8) is given by

$\begin{matrix} L (v, q, u; η_{1}, η_{2}) \\ = \sum_{i, j} g (κ (u_{i, j}^{n})) \sqrt{1 + | v_{i, j} |^{2}} + λ 〈 q - f log q, 1_{Ω} 〉 \\ + 〈 η_{1}, v - \nabla u 〉 + \frac{μ_{1}}{2} ∥ v - \nabla u ∥_{2}^{2} + 〈 η_{2}, q - u 〉 \\ + \frac{μ_{2}}{2} ∥ q - u ∥_{2}^{2}, \end{matrix}$ (5.9) where η₁, η₂ are the Lagrangian multipliers and μ₁, μ₂ are the penalty parameters. It includes three subproblems v, q and u. Then, this requires to iteratively and alternatively solve these three subproblems and update two Lagrangian multipliers until reaching the terminating condition.

The v-subproblem is a typical nonlinear minimization problem, whose solution can be obtained via Newton method. In this way, one can get

$\begin{matrix} v^{n + 1, l + 1} = v^{n + 1, l} \\ - \frac{g (κ (u^{n})) v^{n + 1, l} (1 + | v^{n + 1, l} |^{2})^{- \frac{1}{2}} + μ_{1} (v^{n + 1, l} - \nabla u^{n}) + η_{1}^{n}}{g (κ (u^{n})) (1 + | v^{n + 1, l} |^{2})^{- \frac{3}{2}} + μ_{1}} . \end{matrix}$ (5.10)

The q-subproblem has a closed form solution, which can be given by

$q^{n + 1} = \frac{u^{n} - \frac{λ + η_{2}^{n}}{μ_{2}}}{2} + \sqrt{{(\frac{u^{n} - \frac{λ + η_{2}^{n}}{μ_{2}}}{2})}^{2} + \frac{λ f}{μ_{2}}} .$ (5.11) Here, one needs to solve a quadratic equation.

According to the optimality conditions, it is easy to derive the Euler-Lagrange equation of the u-subproblem. Assume periodic boundary condition and then use the fast Fourier transform (FFT) to solve this equation, one can get

$\begin{matrix} u^{n + 1} = F^{- 1} \\ (\frac{F (η_{2}^{n} + μ_{2} q^{n + 1} - div η_{1}^{n} - μ_{1} div v^{n + 1})}{μ_{2} - μ_{1} F (Δ)}) . \end{matrix}$ (5.12) More details of the numerical algorithm for solving (5.7) can be referred to [42, 46].

5.2.3 Comparison results

In experiments, the “Barbara" and “Triangle" images presented in Fig. 1 are chosen as the test images to show the competitive performance of AEEPD over the TAC method. The testing data are generated by adding Poisson noise with noise level σ = 1 or 5. The related numerical results for Poisson denoising are arranged in Table 4. Specifically, we show the restored “Barbara" images and the corresponding residual images at noise level σ = 1 in Fig. 8. At this moment, the optimal values of parameters used for the TAC method are λ = 38, α = 0.27, μ₁ = 0.39, μ₂ = 7, while the parameters of our AEEPD are set to be λ = 22, α = 0.24, ι = 0.0019, δ = 0.1, β = 0, γ₁ = 0.006, γ₂ = 2.7. In addition, the related visual results at noise level σ = 5 for the “Triangle" image are displayed in Fig. 9. Here, we choose λ = 2.5, α = 0.45, μ₁ = 0.36, μ₂ = 8 for the TAC method and λ = 5.14, α = 0.1, ι = 0.0002, δ = 0.5, β = 0, γ₁ = 0.0035, γ₂ = 0.1 for AEEPD, respectively.

Table 4
Numerical results for Poisson denoising

Image Barbara (512 × 512) Triangle (254 × 214)

Noise Model PSNR SSIM CPU(s) Iter PSNR SSIM CPU(s) Iter

1 TAC 30.147 0.915 52.88 101 38.781 0.987 17.50 139

AEEPD 30.598 0.926 29.68 209 40.509 0.991 14.29 500

5 TAC 25.444 0.802 108.17 206 34.519 0.981 25.62 273

AEEPD 25.790 0.820 43.25 322 37.747 0.984 17.25 500

Image		Barbara (512 × 512)	Triangle (254 × 214)
1	TAC	30.147	0.915	52.88	101	38.781	0.987	17.50	139
	AEEPD	30.598	0.926	29.68	209	40.509	0.991	14.29	500
5	TAC	25.444	0.802	108.17	206	34.519	0.981	25.62	273
	AEEPD	25.790	0.820	43.25	322	37.747	0.984	17.25	500

Fig. 8

(a), (c) Restored “Barbara” images by TAC and AEEPD; (b), (d) their corresponding residual images respectively; (e), (f), (g) the plots of PSNR, relative error and energy versus CPU time (under noise level σ = 1).

Fig. 9

(a), (c) Restored “Triangle” images by TAC and AEEPD; (b), (d) their corresponding residual images respectively; (e), (f), (g) the plots of PSNR, relative error and energy versus CPU time (under noise level σ = 5).

From Figs. 9, the restored images obtained by TAC are not as natural as AEEPD, and more noise are removed by AEEPD. This is owing to the use of the adaptive weighted matrix. We observe from Table 4 that our AEEPD method gets higher PSNR and SSIM values than the TAC method, while spending less CPU time. Also, this can be easily seen from the plots of PSNR versus CPU time in Figs. 9. Simultaneously, the plots of relative error and energy versus CPU time show evidences of algorithmic convergence.

6 Conclusion

In this study, we proposed an adaptive Euler’s elastica Poisson restoration model via an adaptive weighted matrix. The proposed model was solved by an efficient ADMM, where all the subproblems had explicit solutions. Numerical experiments on both natural and synthetic images demonstrated the effectiveness and superior performance of our proposed method by comparing with other state-of-the-art methods.

From the perspective of practical implication, the proposed regularization method can be applied to other problems in image processing and computer vision, such as image segmentation, hyperspectral unmixing and hyperspectral image fusion. A limitation is that a convergence theory for the proposed method is not given, although we conducted a convergence analysis. In the future, we will study the relevant convergence theory.

Data availability

The data used to support the findings of this study are available from the corresponding author upon request.

Competing interests

The authors declare that they have no competing interests.

Conflicts of interest

The authors declare that they have no conflicts of interest.

Funding

This work was partly supported by the Natural Science Foundation of Jiangxi Province (20192BAB211005).

References

Aubert

and Kornprobst

, Mathematical problems in image processing: partial differential equations and the calculus of variations, New York: Springer, 2006.

Aujol

J.F.

, Some first-order algorithms for total variation based image restoration, Journal of Mathematical Imaging and Vision 34(3) (2009), 307–327.

Bayram

and Kamasak

M.E.

, A directional total variation, IEEE Signal Processing Letters 19(12) (2012), 781–784.

Bredies

, Kunisch

and Pock

, Total generalized variation, SIAM Journal on Imaging Sciences 3(3) (2010), 492–526.

Chambolle

, An algorithm for total variation minimization and applications, Journal of Mathematical Imaging and Vision 20(1) (2004), 89–97.

De los Reyes

J.C.

, Schönlieb

C.B.

and Valkonen

, Bilevel parameter learning for higher-order total variation regularisation models, Journal of Mathematical Imaging and Vision 57(1) (2017), 1–25.

Deng

and Yin

, On the global and linear convergence of the generalized alternating direction method of multipliers, Journal of Scientific Computing 66(3) (2016), 889–916.

di Serafino

, Landi

and Viola

, Directional TGV-based image restoration under Poisson noise, Journal of Imaging 7(6) (2021), 99.

Elsey

and Esedoglu

, Analogue of the total variation denoising model in the context of geometry processing, Multiscale Modeling and Simulation 7(4) (2009), 1549–1573.

10.

Fairag

, Chen

and Ahmad

, An effective algorithm for mean curvature-based image deblurring problem[J], Computational and Applied Mathematics 41(4) (2022), 176.

11.

Gao

, Liu

and Yang

, Total generalized variation restoration with non-quadratic fidelity, Multidimensional Systems and Signal Processing 29(4) (2018), 1459–1484.

12.

Guo

and Chen

, Image denoising based on nonconvex anisotropic total-variation regularization, Signal Processing 186 (2021), 108124.

13.

, Wang

and Chen

, A penalty relaxation method for image processing using Euler’s elastica model, SIAM Journal on Imaging Sciences 14(1) (2021), 389–417.

14.

Jiang

, Huang

, Lv

X.G.

and Liu

, Restoring Poissonian images by a combined first-order and second-order variation approach, Journal of Mathematics 2013 (2013), 1–11.

15.

Jiang

, Huang

, Lv

X.G.

and Liu

, Alternating direction method for the high-order total variation-based Poisson noise removal problem, Numerical Algorithms 69(3) (2015), 495–516.

16.

Kongskov

R.D.

, Dong

and Knudsen

, Directional total generalized variation regularization, BIT Numerical Mathematics 59(4) (2019), 903–928.

17.

Lantéri

and Theys

, Restoration of astrophysical images-the case of Poisson data with additive Gaussian noise, EURASIP Journal on Advances in Signal Processing 2005(15) (2005), 1–14.

18.

, Chartrand

and Asaki

T.J.

, A variational approach to reconstructing images corrupted by Poisson noise, Journal of Mathematical Imaging and Vision 27(3) (2007), 257–263.

19.

and Liu

, Anisotropic total generalized variation model for Poisson noise removal[J]. Multimedia Tools and Applications 2023:1–14.

20.

M.M.

and Li

B.Z.

, A novel weighted anisotropic total variational model for image applications, Signal, Image and Video Processing 16(1) (2022), 211–218.

21.

Liu

, Huang

T.Z.

, Lv

X.G.

and Wang

, High-order total variation-based Poissonian image deconvolution with spatially adapted regularization parameter, Applied Mathematical Modelling 45 (2017), 516–529.

22.

Liu

, Yan

and Zeng

, Surface-aware blind image deblurring, IEEE Transactions on Pattern Analysis and Machine Intelligence 43(3) (2019), 1041–1055.

23.

X.G.

, Jiang

and Liu

, Deblurring Poisson noisy images by total variation with overlapping group sparsity, Applied Mathematics and Computation 289 (2016), 132–148.

24.

, Yu

and Zeng

, Sparse Representation Prior and Total Variation–Based Image Deblurring under Impulse Noise, SIAM Journal on Imaging Sciences 6(4) (2013), 2258–2284.

25.

Pang

Z.F.

, Zhang

H.L.

, Luo

and Zeng

, Image denoising based on the adaptive weighted TV^p regularization, Signal Processing 167 (2020), 107325.

26.

Pang

Z.F.

, Zhou

Y.M.

, Wu

and Li

D.J.

, Image denoising via a new anisotropic total-variation-based model, Signal Processing: Image Communication 74 (2019), 140–152.

27.

Rudin

L.I.

, Osher

and Fatemi

, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena 60(1-4) (1992), 259–268.

28.

Scherzer

, Handbook of mathematical methods in imaging, Springer Science and Business Media, 2010.

29.

Setzer

, Steidl

and Teuber

, Deblurring Poissonian images by split Bregman techniques, Journal of Visual Communication and Image Representation 21(3) (2010), 193–199.

30.

Shen

, Kang

S.H.

and Chan

T.F.

, Euler’s elastica and curvature-based inpainting, SIAM journal on Applied Mathematics 63(2) (2003), 564–592.

31.

Tai

X.C.

, Hahn

and Chung

G.J.

, A fast algorithm for Euler’s elastica model using augmented Lagrangian method, SIAM Journal on Imaging Sciences 4(1) (2011), 313–344.

32.

Thanh

D.N.H.

, Prasath

V.B.S.

, Hieu

L.M.

and Dvoenko

, An adaptive method for image restoration based on high-order total variation and inverse gradient[J], Signal, Image and Video Processing 14(6) (2020), 1189–1197.

33.

Wang

and Pang

Z.F.

, Image denoising based on a new anisotropic mean curvature model[J]. Inverse Problems and Imaging 2023.

34.

Wang

, Yang

, Yin

and Zhang

, A new alternating minimization algorithm for total variation image reconstruction, SIAM Journal on Imaging Sciences 1(3) (2008), 248–272.

35.

Wen

Y.W.

, Chan

R.H.

and Zeng

T.Y.

, Primal-dual algorithms for total variation based image restoration under Poisson noise, Science China Mathematics 59(1) (2016), 141–160.

36.

and Tai

X.C.

, Augmented Lagrangian method, dual methods, and split bregman iteration for ROF, vectorial TV, and high order models, SIAM Journal on Imaging Sciences 3(3) (2010), 300–339.

37.

Yang

, Ma

, Zhang

and Wang

, Noise removal using an adaptive Euler’s elastica-based model, The Visual Computer, 2022.

38.

Yashtini

and Kang

S.H.

, A fast relaxed normal two split method and an effective weighted TV approach for Euler’s elastica image inpainting, SIAM Journal on Imaging Sciences 9(4) (2016), 1552–1581.

39.

Zhang

, Wang

, Duan

, Li

, Hu

and Yan

, Euler’s elastica strategy for limited-angle computed tomography image reconstruction, IEEE Transactions on Nuclear Science 64(8) (2017), 2395–2405.

40.

Zhang

and Wang

, Edge adaptive directional total variation, The Journal of Engineering 2013(11) (2013), 61–62.

41.

Zhang

, Chen

, Deng

and Wang

, Fast linearized augmented Lagrangian method for Euler’s elastica model, Numerical Mathematics: Theory, Methods and Applications 10(1) (2017), 98–115.

42.

Zhang

, Li

, Yang

, Ma

and Deng

, Poisson image restoration using a novel directional TV^p regularization, Signal Processing 193 (2022), 108407.

43.

Zhang

, Ma

, Wu

and Deng

, High-order total bounded variation model and its fast algorithm for Poissonian image restoration, Mathematical Problems in Engineering 2019 (2019), 1–11.

44.

Zhang

and Yang

Y.F.

, Nonlinear multigrid method for solving the LLT model, Applied Mathematics and Computation 219(10) (2013), 4964–4976.

45.

Zhang

, Ng

M.K.

and Bai

, A fast algorithm for deconvolution and Poisson noise removal, Journal of Scientific Computing 75(3) (2018), 1535–1554.

46.

Zhong

, Yin

and Duan

, Image reconstruction by minimizing curvatures on image surface, Journal of Mathematical Imaging and Vision 63(1) (2021), 30–55.

47.

Zhu

and Chan

, Image denoising using mean curvature of image surface, SIAM Journal on Imaging Sciences 5(1) (2012), 1–32.

48.

Zhu

, Tai

X.C.

and Chan

, Augmented Lagrangian method for a mean curvature based image denoising model, Inverse Problems and Imaging 7(4) (2013), 1409–1432.