Euler’s elastica and curvature based model for image restoration

Minimization functionals related to Euler’s elastica energy has a broad range of applications in computer vision and image processing. This paper proposes a novel Euler’s elastica and curvature-based variational model for image restoration corrupted with multiplicative noise. It combines Euler’s elastica curvature with a Weberized total variation (TV) regularization and gets a novel Euler’s elastica energy and TV-based minimization functional. The combined approach in this variational model can preserve edges while reducing the blocky effect in smooth regions. The implicit gradient descent scheme is applied to efficiently finding the minimizer of the proposed functional. Experimental results demonstrate the effectiveness of the proposed model in visual improvement, as well as an increase in the peak signal-to-noise ratio, compared to the PDE-based methods.


Introduction
During some phases of the manipulation of an image, some random noise is usually introduced. Therefore, image restoration is the fundamental problem in image processing. Among the restoration models, the variational model has been extremely successful in a wide verity of image restoration problems and remains one of the most active areas of research in mathematical image processing and computer vision. The problem includes additive noise removal and multiplicative noise removal. Multiplicative noise appears in various image processing applications, such as synthetic aperture radar (SAR), medical images, single-particle emissioncomputed technology, and positron emission tomography. The additive noise problem can be molded as where f is the observed image, g is the original image and η 1 is the additive noise. In literature, many effective numerical methods have been utilized to tackle such models connected with image de-noising having additive noise, for instance, in [1][2][3][4][5][6]. In the last two decades, variational techniques have been widely studied and investigated for image processing tasks, some of them are related image de-noising. For further details, see [7][8][9][10][11][12][13][14].
Multiplicative noise removal problem can be stated as follow; f ðx; yÞ ¼ gðx; yÞ Z 2 ðx; yÞ; ð2Þ where f is the observed image, g is the original image and η 2 is multiplicative noise. In this direction, researchers have used various total variation-based approaches to solve the multiplicative noise removal problem. The interested reader is referred to [15][16][17][18][19][20] for more details.
The main advantage of TV-based regularization is that it preserves the edges well, but the images resulting from the application of this technique in the presence of noise lead to piecewise constant function. Thus, the finer details in the original image may not be recovered satisfactorily, and the ramps will give stairs (piecewise constants) even though some efforts have been made to reduce the staircase effect. For further details, see [21][22][23][24][25][26].
Euler's elastic was first introduced as a prior curve model by Mumford [27] in computer vision. Then, the Euler's elastica curvature-based total variation-based regularization was applied to find image inpainting [28,29], and useful results were obtained in image inpainting and restoration.
In this paper, we develop a new model with an improved regularized term, i.e., Euler's elastica curvature-based TV regularizer with fourth-order, nonlinear Euler-Lagrange equation. This model restores the images well and substantially reduces the staircase effect while preserving the edges and textures.
The rest of the paper is organized as follows. In section 2, the literature has been reviewed regarding image processing by using TV-based methods and Euler's elastica and curvaturebased energy in image restoration and inpainting models. Section 3, discusses Huang et.al [19] model for the removal of multiplicative noise. The proposed model which is Euler's elastica and curvature-based total variation model for multiplicative noise removal has been discussed in section 4. Section 5, presents the discretization and numerical implementation of the proposed model. Section 6, describes experimental results to compare the two models regarding the visual quality of restoration (PSNR). Section 7, contains the comparison of our proposed model with two other variational based models for image restoration. The sensitivity of the parameters of the proposed model is explained in Section 8. Section 9, concludes the paper. Finally, the derivation of the proposed Euler-Lagrange PDE has been given in appendix A.

Background work
In the last few years, many approaches have been used to solve the multiplicative noise removal problem [18,30,31]. One of these approaches is the local linear minimum mean square method [32][33][34][35], another approach is the anisotropic diffusion method [36][37][38]; these strategies are based on the statistical information of images and noise, so we have not discussed them in details. Our focus in this paper is on variational methods.
The variational approach is an important paradigm for solving the image de-noising problems when the image is defined in the continuous domain. The total variation is a powerful notion in such variational problems. In recent years, variational methods have got much attention in reducing the multiplicative noise owing to the use of total variation (TV) and nonlocal total variation (NLTV) regularization [10,39,40]. In literature, many total variation-based models have been proposed, see [41][42][43][44][45][46] for more details. The variational-based NLTV models have also been widely applied in image restoration [47,48]. Steidl and Teuber [48] employed the NLTV as a regularizer to recover multiplicative noise images. Since the NLTV-norm used the relevant image patches and hence gives good qualitative results.
As, in image model, TV(g) = R O |rg|dx does not take into account that our visual sensitivity to the regularity or local fluctuation |rg| depends on the ambient intensity level g [49]. Since all images are eventually perceived and interpreted by the Human Visual System (HVS); as a result, many researchers have found that the human vision psychology and psychophysics play a significant role in image processing. The classic example is the using of the Just Noticeable Difference Model (JND) in image coding and watermarking techniques [50,51]. In these fields, the NDJ model is used to control the visual perceptual distortion during the coding procedure and watermark embedding. Weber's law was first described by Weber [52]. The law reveals the universal influence of the background stimulus g on human's sensitivity to the intensity increment |rg|, or so called JND, in the perception of the both sound and light: Accoring to Weber's law [49,53], when the mean intensity of background is increasing with high value, then the intensity increment rg also has high value. In literature [19], the authors proposed a nonconvex variational functional of model (2) for multiplicative noise removal:ĝ where the first term is the regularization term and the second term is the nonconvex fitting/ fidelity term following the MAP estimator for multiplicative Gama noise. β 1 and β 2 are the two parameters, f > 0 in L 1 (O) is the given data. The first regularization term is TV functional, which preserve important structures such as edges, an important visual cue in human and computer vision. The second term E w ðgÞ≔TVðloggÞ ¼ R O jrgj g is the well known Weberized TV regularization term. To briefly explain the role of this term, we assume that g has a gradient rg 2 L 1 (O) 2 , then TV(logg) = R O |rg|dx and the Weberized local variation is which encodes the influence of the background intensity according to Weber's law (3). Recently many variational models involving higher order derivatives have been widely used in image processing because they reduce the staircase effect during the noise elimination. The use of Euler's elastic energy minimization model based fourth order derivatives damps out the high-frequency components of images faster than faster than the second order PDE based methods because the associated PDE to the Euler's elastica minimization is fourth order. So the Euler's elastica model can reduce the staircase effect, textures and can produce better approximations to the natural image. Indeed, it is also able to preserve the object edges while erasing noise. The Euler's elastica model is one of the most famous high orders models. It has been successfully applied in various problems image denoising, inpainting, and zooming. In [28] the authors derived are the Euler-Langrage equation, proposed some numerical schemes to solve the corresponding Euler-Language equation and also explained the effect of parameters α 1 and α 2 in Eq (3) in image inpainting.
The Euler' elastica can be described by using the curvature κ of the smooth curve Γ as follows where s is the arc length and α 1 and α 2 are the two positive parameters. In the above functional (6), the first term minimizes the total length while the send term minimizes the power of total curvature, where p can be set p = 1, 2 in [28,54]. In this work, we set p = 2. The Euler's elastica of all the level curves of an image g can be defined as where γl represents the level curve with g = l. The curvature κ can be expressed as Combining the Eqs (7) and (8) and using the co-area formula we get For image restoration applications, the Euler's elastica energy (9) can be used as a regularization term.
The Elastica model [28], minimizing the Euler's elastica energy for image inpainting problem is proposed in the following minimization functional where α 1 and α 2 are arbitrary positive constants, λ > 0 is a penalty parameter, p = 2 is usually chosen, g = g(x, y) is the true image to be restored and k ¼ r Á rg jrgj is the curvature. The virtue of Eq (10) is that the regularization using the Euler's elastica energy penalizes the integral of the square of the curvature along edges instead of only penalizing the length of the edges as the TV model does (if taking α 2 = 0) [10]. Motivated by the applications of Euler's elastica in image inpainting and restoration and Weberized TV-regularization in image restoration, we propose a new model based on the theory of Euler's elastica energy and Weberized TV-regularization for multiplicative noise removal problem.

Li-Li Huang model (M1)
Li-Li Huang et al. proposed the non-convex multiplicative noise removal model using the total variation (TV) filter in [19]. This model achieved some useful restoration results.
The minimization functional by this approach for model (2) is given in (4) which is provided as follow; Here, f > 0 in L 1 (O) is the given data in the model. In (11) the first term is the TV-functional which preserve the important structures such as edges, an important cue in the human and computer vision and its second term is called Weberized TV regularization term [53]. We suppose that g has the gradient rg 2 L 1 (O) 2 , then the above Eq (12) can be re-written as with is the Weberized local variation, which encodes the influence of backward intensity according to Weber's law [53]. The corresponding Euler-Lagrange of the minimization functional (11) can be defined as Since g > 0, so the above Eq (15) can be re-defined as Definel ¼lðgÞ ¼ 1 gðb 1 gþb 2 Þ , then Eq (17) implies with Neumann boundary conditions. Additive operator splitting (AOS) method [19,55] has been utilized to solve (18).

The proposed model (M2)
In this section, we aim to introduce a new model using both Euler's elastica curvature energy and Weberized total variation (TV)-norm as a regularizer. This model apparently uses the advantages of both Euler's elastica curvature energy and Weberized TV-norm, which leads to good restoration results. In this model, the Euler's elastica energy, denoised images have smooth connections in the level curves of images which lead to good performance regarding image restoration (PSNR), elimination of staircase effect, and preservation of textures, while the Weberized total variation preserves jump discontinuities and sharp the edges in images. This is believed to be better estimation to a natural image than a piecewise constant approximation in the smooth regions. Hence using this model, consistent improvement in PSNR values is obtained.
As the minimization approach for model (2) by Huang et al. [19] is defined aŝ where R(g) = |rg|. But Euler's elastica and curvature-based inpainting model [28] is given by the equation. having where α 1 and α 2 are the parameters and κ is the curvature and is defined as follows.
Hence, (23) is a new high-order curvature-based total variational minimization functional for multiplicative noise removal problem. The first term is called regularization term and the second term is called non-convex data fitting term, where β 1 , β 2 , α 1 , and α 2 are the four regularization parameters which usually depends on the noise level and type of the image. The corresponding Euler-Lagrange equation of (23) is given as and with boundary conditions @g @n ¼ 0 and @ða 1 þ a 2 k 2 Þ @n ¼ 0. Where rg in (28) is the normal vector and r ? g is the corresponding tangent vector i.ẽ From Eq (26) U = (U 1 , U 2 ), where and

Numerical implementation
To discretize the Eq (26), which is a fourth-order PDE, we use the finite difference method discussed in [56][57][58]. Here, we include the details for the sake of the completeness of the present discussion.
In this discussion, we use the notation u (i, j) = u(iΔx, jΔy), where i, j are integers, Δx, Δy are the space step sizes along the x and y directions, respectively; and (iΔx, jΔy) represents a discrete point. Here At(i, j) we can write Eq (32) Using the Eq (33), we approximate U 1 iþ 1 2 ;j and U 1 iÀ 1 2 ;j as follow; Curvature term. These are approximated by min-mode of two adjacent whole pixels.
where min À mod m; Partial derivatives in x. By the central difference of two adjacent whole @x g iþ 1 Here Partial derivatives in y. By min-mod of @y, s at two adjacent whole points @y g iÀ 1 and Also As the steady state form of Eq (26) is At (i, j), we approximate Eq (50) by using (33) Since, implicit gradient descent scheme has achieved good results in [14], so we utilize this scheme to solve the PDE (51), which can be expressed as under; where g n ði;jÞ shows the value of g at ndt times, dt is time step size, and Δx = Δy.

Numerical experiments
In this section, we provide some numerical results from applying our proposed model M2. We also compare them with the results obtained by using the model M1. The test images used for numerical experiments are "Moon", "Rose", "Synthetic1", "Synthetic2", "Synthetic3", "Syn-thetic4" and "Lena", "Boat", "House", "Peppers", "Baboon", and "Texture", respectively, which are shown in Fig 1. We test these images on speckle noise (uniform distribution) with mean value 1 and variance σ 2 .

Quality of restoration
Peak signal to noise ratio (PSNR) is used to analyze the quality of the image that has been restored. So here we check the restoration quality of the two models by PSNR as follow: whereĝ is the given image, g is the restore image and M × N the size of the image. Example 1: In this first example, the two methods M1 and M2 are applied and tested on the natural images "Moon" and "Rose," having speckle noise with noise variance σ 2 = 0.1, and σ 2 = 0.1, respectively, which are shown in Figs 2 and 3. In both Figs, (a) and (b) are the original and noisy images, while Figs (c) and (d) depict the restored images by methods M1 and M2, respectively. In each case, we can notice that the visual quality of restoration by proposed method M2 is much better than that of method M1. Moreover, the PSNR values for the two images "Moon" and "Rose" for methods M1 and M2 are also listed in Table 1. The bigger the PSNR value, the better the de-noising performance. It can be seen from Table 1 that the PSNR values of procedure M2 are greater than that of method M1 for the two images, which shows the better restoration performance of M2 over M1. Hence, the results in Figs 2 and 3, and Table 1 can show that our proposed method M2 can improve the visual quality of the restored images and PSNR significantly better than M1. The optimal experimental values of the parameters of our technique M2 (β 1 , β 2 , α 1 , α 2 ), for the two images "Moon" and "Rose" are (0.003, 0.05, 0.40, 1.38) and (0.007, 0.01, 0.37, 1.41), respectively. In this example, we also set 4x = 4y = 10, and dt = 0.0002.

Example 3:
Here, the image size of "Lena", 256 2 , is taken as the test image. The main logic of our proposed model is to use Euler's elastic curvature and Weberized TV-norm to recover both jumps and smooth signals accurately. In Figs 8, 9 and 10, and Table 1, we can observe that M2 yields better image restoration results than M1 since it preserves edges and minimizes the staircase effect.

Example 5:
In this example, the homogeneity is checked, and loss (or preservation) is examined for the two techniques M1 and M2 while being applied to "Lena". For this purpose, different lines of the original image are compared to noisy and restored images that are shown in Fig 12. It is clear that the line restored by proposed method M2 (shown in Fig 12(c)) is far better than what is acquired utilizing method M1 that is presented in Fig 12(b).

Example 6:
In this example, we test our proposed model M2 for different values of time step dt for a real image for the same values of parameters (β 1 , β 2 , α 1 , α 2 ). The values of these parameters for this example are selected as (0.007, 0.093, 0.23, 1.29). So, from the Fig 13 and Table 2, we can Image restoration using Euler's elastica energy and total variation regularization  Therefore, it is reasonable to conclude that Euler's elastica curvature-based model M2 is best in the sense that it has piecewise smooth intensities, has sharp edges, and minimizes the staircase effect better than model M1.
where ξ = log(f) and v is a vector field of the image ξ. Also λ > 0, α > 0, and β > 0 are the regularizer parameters. Here, f 0 is the observed image and f is the original image. The alternating direction method has been used to solve Eq (54), which can be written into the following two minimization subproblems: x kþ1 ¼ arg min For the solution of the subproblem (55), the authors used the Chambolle's dual method which is given as follows; where q ¼ q 11 q 12 , p 1 = (q 11 , q 12 ), p 2 = (q 21 , q 22 ), x k x ; x k y represent first order forward and backward differences, respectively. Then, the authors solved p 1 and p 2 by fixed point iteration, which are given as under.
where, p k;0 1 ¼ p k;0 2 ¼ 0; and l = 1. For the solution of the subproblem (56), the authors used Bregman method by letting η = rξ, which leads to the following unconstrained problem.
where μ > 0 is a penalty parameter. The split Bregman algorithm is defend as under; x kþ1 ¼ arg min and T represents the thresholding operator defined by The solution of the first term of Eq (62) is given as under: The closest solution is obtained as follow.
In this subsection, we have compared the two models, i.e., method M3 and proposed method M2 for image restoration for the same images with the same size and noise variance along with same parameters values as selected in [59]. We can see that the results obtained by our proposed method M2 are outstanding in the visual quality of restoration (PSNR), eliminating the staircase effect and preserving the textures. These obtained results are shown in Figs 14 and 15, and Table 3, respectively. The optimal values of parameters for our proposed method M2 (β 1 , β 2 , α 1 , α 2 ) for the two images "Boat" and "House" are (0.01, 0.008, 0.23, 1.05), and (0.01, 0.0062, 0.17, 1.04), respectively. In this case, we choose 4x = 4y = 10, and dt = 0.0002.  [60] and hence some good restoration results have been obtained. The minimization functional for this model is given as follow: where the second term is called the weighted total variation term, function b : O ! ½0; ffiffi ffi a p ; is the positive parameter such that 0 < < 1. Also the second term  and (1 + ) R O |Du|. The authors used an alternating iteration process directly of (67) which includes two minimization problems as follows.
First by fixing b, they solve Second, with u fixed, they solve The given algorithm shows the alternating iteration process for the task of removing multiplicative noise. The only challenging issue is to minimize the functional (70) in Algorithm 2, or to minimize the general functional in the minimization problem (68). For this purpose, the authors convert the minimization functional (68) equivalently for constructing more efficient solver, which is given as under.
The authors then used the alternating direction of multipliers(ADMM) method to solve the constrained minimization problem (73) which is shown in the following iteration process in Algorithm 3.

Algorithm 3: Algorithm for model Minimization functional (73)
Initialization process: Given parameters λ and μ, set the iteration index n = 0, and let u 0 = v 0 , d n = 0; Iteration process: Determination process: Once the sequence {u n } and {v n } converge when n ! 1, stop the algorithm, otherwise, set n = n + 1 and return to the iteration process.
The Euler Lagrange equation of the functional in the minimization problem (74) is given as under; The solution of above Eq (77) is obtained by using the few times Newton iterations. The solution to the (75) is obtained by the following Chambolle's projection algorithm; where q: O ! R × R shows the vector function, div notation shows the divergence operator. The function q is the limit of the sequence q m (the index m differs from the outer loop index n) generated from the following fixed point iteration: given an initial q 0 and a time-step %, the following iterative scheme is used: According to the authors, due to the nonconvex nature of minimization problem (69) it is hard to prove the global convergence about the whole algorithm in Algorithm 2. For, further details, see [60].
Here, the model M4 that is proposed in [60] is compared with our proposed method M2 for the same images having the same size and noise variances and same parameters values that have been selected in [60]. Again, from Figs 16 and 17, and Table 4 we can observe that our proposed technique M2 has better performance in the visual quality of restoration (SNR), reducing the staircase effect and preserving the textures compared to the model M4. The values of the parameters selected for our proposed model M2 (β 1 , β 2 , α 1 , α 2 ) for the two images "Peppers" and "Lena" are (0.01, 0.0085, 0.24, 1.05), and (0.02, 0.005, 0.20, 1.02), respectively. In this case, we select 4x = 4y = 10, and dt = 0.0002.

Sensitivity analysis of parameters
To briefly comment on the choice of the various parameters used in the algorithm of our model M2, it must be noted that β 1 , β 2 , α 1 and α 2 are more complicated to choose according to our experience. The difficulty of tuning these parameters is that they not only depend on the noise level, but also on the type of images. However, their optimal values are adjusted and tuned according to the noise variance and the image. It was observed that the ranges of values allowed are: β 1 2 [0.0025, 0.0097], β 2 2 [0.0082, 0.074], α 1 2 [0.29, 0.50] and α 2 2 [1.22, 1.61] for natural and synthetic images according to the noise variance σ 2 = 0.1. Using, these ranges, better restoration results could be achieved with improved PSNR results. Results are presented in Tables 5 and 6. For brevity, the following notations are utilized in these tables.

Conclusion
In this paper, a new high-order model was introduced using Euler's elastica curvature combined with the total variation regularization for image restoration with speckle noise. The implicit gradient descent scheme was exploited for solving nonlinear PDE arisen from the minimization of the proposed functional. The experimental results demonstrated that the Table 5. The PSNR value of the restored image "Moon" with optimal values of β 1 , β 2 , α 1 and α 2 is 26.23. Parameter sensitivity analysis of our proposed model M2 by percentage increased in values of the parameters β 1 , β 2 , α 1 and α 2 , with the resultant percentage increase or decrease in PSNR of the de-noised image of size (300 2 ). proposed model improves PSNR and can preserve edges, textures and minimize the staircase effect compared with existing PDE-based models. The sensitivity analysis of parameters was also discussed in details. Our model has also been compared with other variational PDE-based models for image restoration when the noise variance is large. In our model the resulting Euler-Lagrange equation has fourth order derivatives and is also anisotropic and highly nonlinear, and thus the conventional algorithms struggle to solve it efficiently due to the stability restriction. This problem is under intense study and results will be reported in the subsequent paper.