New dual method for elastica regularization

The Euler’s elastica energy regularizer has been widely used in image processing and computer vision tasks. However, finding a fast and simple solver for the term remains challenging. In this paper, we propose a new dual method to simplify the solution. Classical fast solutions transform the complex optimization problem into simpler subproblems, but introduce many parameters and split operators in the process. Hence, we propose a new dual algorithm to maintain the constraint exactly, while using only one dual parameter to transform the problem into its alternate optimization form. The proposed dual method can be easily applied to level-set-based segmentation models that contain the Euler’s elastic term. Lastly, we demonstrate the performance of the proposed method on both synthetic and real images in tasks image processing tasks, i.e. denoising, inpainting, and segmentation, as well as compare to the Augmented Lagrangian method (ALM) on the aforementioned tasks.


Introduction
Traditional variational methods have been extensively applied to image restoration problems based on image features, such as texture, edge, and region, etc. [1][2][3][4]. In particular, the combination of the high-order TV term and Euler regularizers in variational models addresses certain problems that cannot be addressed by low-order models. However, the complexity of the terms makes the models more difficult to implement. It has become a challenge in recent years to design simpler and more effective solutions for the combined model. Before presenting our new work, we introduce the TV term and the Euler's elastica term below.
The problem of image restoration is finding the restored image u = u(x) given the damaged image f = f(x), x 2 < d , where < is a bounded domain with a Lipschitz boundary. For a grayscale image, the image repair model is f = Au + η, where η is the noise information and A is a blurring operator [5][6][7]. The image restoration problem can then be formulated as the minimization of the following energy functional, However, the minimization problem in (1) is ill-posed. To solve this, Tikhonov et. al [8] proposed a regularization technique. By adding a smoothing regularizer into the energy functional, the problem will obtain a unique solution. The side effect is that the model can no longer preserve edges in the image. The later proposed Rudin-Osher-Fatemi (ROF) model retains image edges by solving for a piecewise constant function in the space of bounded variation functions (BV). Nowadays, many methods based on TV regularization are used to deal with imaging problems such as image denoising [9][10][11][12] and image segmentation [13,14]. Another downside of the TV models is that results are often accompanied by blocky (staircase) effects and loss of image contrast [15][16][17][18]. Recently, scholars have proposed many solutions such as iterative regularization techniques [19] and the use of other high-order terms to mitigate the problem. Faster implementations have also been invented [20].
As for the Euler's elastica [21], it has attracted much attention due to its good properties in mathematical and physical systems. The Euler's elastica energy functional is defined as where, n ¼ rn jrnj . The term was first used in computer vision by Mumford [22], and has since proven to be effective in solving the problems present in the TV model. The Euler's elastica has also been widely applied in various fields of image processing such as image denoising [23][24][25], image segmentation [20,[26][27][28], inpainting [24,[29][30][31], illusory contour [32,33], and segmentation with depth [34][35][36]. Therefore, we believe it is important to design an efficient numerical solution for the combined model.
Due to the non-convexity, non-smoothness, and high-order of the derivatives of 2, it is a challenging task to design a fast and efficient solution. The ALM method has achieved good results in optimizing 2, Tai et al. [23] first proposed this ALM method to solve the image inpainting problems, then Zhu et al. [28] extend the ALM method to image segmentation field. So far, the primal-dual technique [20,37] performs better in optimizing 2.
In this paper, we propose a new primal-dual method for the solution of 2, that makes it easier to use. The key points of the proposed method can be summarised below: (i) We introduce the dual variable p to circumvent the curvature term. (ii) Using appropriate indicator functionals, we reformulate 2 as a minimization problem of u and a maximization of p. (iii) The subproblem u has an analytic solution, and the subproblem p can be solved by a gradient descent algorithm. Numerical experiments demonstrate the improved efficiency of the proposed method.
Compared to the ALM method, the advantages of this method can be summarized into three points, (i) The proposed method only introduces one dual variable p, while the ALM method introduces 8 intermediate variables. (ii) Due to the fewer variables, the proposed method has a weaker dependency on parameters. (iii) The CPU running time required for each iteration is greatly reduced.
The rest of this article is organized as follows. In the next section, we introduce the previous models and the associated numerical algorithms. In section 4, we propose our model for image denoising and image segmentation. The subproblems of energy minimization are solved in section 5. In section 6, we provide some numerical results to illustrate the effectiveness of the new algorithm. The last section presents the conclusions.

The TV model for image denoise and inpainting
The well-known TV model [2] for image denoising is an energy minimization problem on u, such that where, γ is a penalty parameter for the summed length of the curves. To use the model in image inpainting, we need to incorporate a mask function η, which is defined as where, x L = {x 1 , x 2 . . .x l } denote the damaged regions. The classical TV model combined with this mask function is It is evident that if η is the identity matrix, then the inpainting model above is the same as the denoising model. Therefore, we only focus on the image restoration model. The evolution equation of u can be derived via variational methods as @uðx; tÞ @t

The TV model reformulated via the dual method
To simplify the TV model (3), we introduce the dual variable p to circumvent the curvature term. Substituting p into the TV model, we get By using the dual method, we can successfully avoid the curvature term and significantly simplify calculations. The new evolution equations of u are @uðx; tÞ @t and p can be solved from @pðx; tÞ @t ¼ ru

The Euler's elastica model for image denoising
In order to recover edges and counter the staircase effect, Tai et al. [23] proposed the Euler's elastica model The boundary produced by this method is curved rather than straight. Its solution via the ALM proposed by Tai et al. [23] simplifies the calculation and increases the optimization efficiency. The Tai-Hahn-Chung (THC) formulation is Eðu; w; n; m; q; pÞ ¼ where, the functions d R ðvÞ and d R ðmÞ denote the constraints |m| � 1 and 0 � u � 1 respectively.

The Chan-Vese model with elastica for image segmentation
The task of two-phase segmentation The Chan-Vese model, a classical two-phase segmentation model, [38] is a reduced piecewise constant Mumford-Shah model [4] under the variational level set framework. Its form is In the above model, ϕ is the level set function and H(ϕ) is the Heaviside function of ϕ(x), stated as Replacing the TV regularizer with the Eular's elastica energy in the Zhu-Tai-Chan (ZTC) formulation [28] leads to the Chan-Vese model with elastica (CVE) below Adding another variable that relaxes the Heaviside function, u = H(x) and u 2 [1, 0], we can construct the following augmented Lagrangian functional Eðu; w; n; m; q; pÞ The functions d R ðvÞ and d R ðmÞ denote the constraints |m| � 1 and 0 � u � 1 respectively.

Image denoise and inpainting
Combining (7) and (10), we propose the dual formulation of the Chan-Vese model with elastica This minimization problem can be divided into two subproblems, and their solutions can be expressed as follows: u − subproblem. This subproblem is a minimization problem, and the objective function of optimization is We find that this function is almost identical to the TV model and get the same solution: In this formula, u has an analytic solution, so there is no need to iterate further and running time is greatly reduced. p − subproblem. The dual variable p can be solved by @pðx; tÞ @t ¼ gru where the value of g is directly determined by the u obtained in the previous step. This is a good way to avoid fourth-order terms and can solve the Euler's elastica term better. Since no additional parameters are introduced, and the iteration time required for each step is greatly reduced compared to the ALM solution.

Image segmentation
Combining (7) and (14), we propose the Chan-Vese model with elastica reformulated with the dual method shown below, The Chan-Vese model is a two-term segmentation model, c 1 and c 2 are mean values of the foreground and background, u − subproblem. Similar to (18), we can also obtain an exact solution for u, p − subproblem. Although the meaning of u has changed, p can still be solved the same way as (19).
In this section, we presented the dual formulation of the CVE model for image denoising and segmentation. Next, we will design the discretized numerical algorithms for the models.

Image denoising and inpainting
To compute the two subproblems numerically, we need to design discrete algorithms for each problem. For the sake of simplicity, we discretize the image domain pixel by pixel with the rows and column numbers as indices. Then, the gradients can be represented approximately by forward, backward, and central finite differences, where, The Euler's Elastica term of u can be stated as The other variables can be expressed in similar ways. Next, we will give a detailed explanation of the solutions to the subproblems obtained in section 4.
u − subproblem. The partial derivative with respect to E gives the analytic solution with respect to u as follows The u in the image inpainting model can be formulated as follows, In fact, this formulation is the same as the corresponding one in the solution of the TV model using the dual method.
p − subproblem. The dual variable p also can be solved by where no additional parameters have been introduced compared to the ALM, so the iteration time required for each step is greatly reduced.
In each iteration, the following error tolerance should be checked to determine convergence, i. e., where, Tol = 0.001. S k+1 are defined as In summary, the denoising al-gorithm is shown in Algorithm 1, and the inpainting algorithm is shown in Algorithm 2.

Image segmentation
In this subsection, we will apply discretization to solve the formulas obtained in the previous image segmentation section. c 1 , c 2 − mean value. (21) can be solved by u − subproblem. Same as TV model, u in (22) can be computed by p − subproblem. Here, the Euler's elastic term is solved directly by the differential equation. We can get Algorithm 3 is the summary of the dual elastica segmentation method.

Numerical experiments
In this section, we show the results from numerical simulations to illustrate the effectiveness of our algorithms in image denoising, inpainting, and segmentation. All experiments are running in MATLAB R2020a.

Testing of Algorithm 1 and Algorithm 2
First of all, we tested our algorithm on three synthetic images. The size of the synthetic images is 256 x 256.
The image denoising result is shown in Fig 1, where (a) is the picture of a clipped triangle with Gaussian noise, (b) is the result of Algorithm 1, and (c) is the relative error, i.e. noise removed through denoising. Fig 2 shows  Besides the qualitative performance of our algorithm in image deniosing, inpainting, and segmentation, computational efficiency is another major point of interest. In Fig 1, the algorithm took 9 iterations to satisfy (28) and the CPU running time was 0.0230 seconds to achieve a PSNR score of 27.9992. In Fig 2, our algorithm ran for 500 iterations and each iteration took 0.0034 seconds, so the total CPU time was 1.7048 seconds. The last segmentation experiment

Comparisons to ALM in real image
In this section, we compare our proposed method against the ALM solution. The images used in the denoising and inpainting experiments are taken from the Set12 dataset [39], Set14 dataset [40]and the BSD68 dataset [41], and the images in segmentation experiments are from the COVID-CT dataset [42] and PASCAL-VOC2012 dataset [43].
In Fig 6, both the THC algorithm and the proposed algorithm successfully removed noise and avoided the staircase effect. However, the THC model is more difficult to tune due to having additional parameters. As seen in the qualitative results, the dual method maintained edges better.
In order to further compare with the ALM method, we present some quantitative results and convergence times. First, we compare the convergence speed of the two methods in Fig 7. The convergence rate of the proposed method is slightly slower than that of the THC method. The main reason is that the p k in the calculation of u k+1 in 24 uses u k , but the ALM method does not. In Table 1, we compare the similarity of the denoising results to the ground truths (via PSNR) and the efficiency (via the number of iterations to reach convergence) of the two algorithms over denoising the three images in Fig 6. Four resolutions of each image were used to more thoroughly investigate the two methods. Results show that proposed method requires more iterations to achieve the same quality of denoising, but requires less time each iteration compared to the THC algorithm. For example, for the 256 � 256 image of 'Baboon', the THC method required 0.081 seconds per iteration whereas the Algorithm 1 only needed 0.013 seconds. The reason why our algorithm requires less time per iteration is the reduction in variables, as our algorithm only uses two variables u and p, and neither u nor p needs to be solved iteratively. On the whole, the dual elastica algorithm is far better in efficiency with the same qualitative results.    Table 2, less time is required of our algorithm to produce similar results as the ZTC.
In the next part, we will segment some COVID-CT images to show the effectiveness of our two-phase segmentation algorithm. We used the ZTC algorithm [28] for comparison with our  algorithm. Fig 9 shows two examples of two-phase segmentation of real images. The column (a) is the initialization of the level set function and the original picture, (b) shows the results of the ZTC method, and (c) shows the results of the proposed method. Visually, those results appear similar. However, as shown in Table 3, our algorithm is more efficient in both the execution time per iteration and the total number of iterations.
To further compare the segmentation results numerically, we use the Dice metric to measure the segmentation quality as,

PLOS ONE
where, N gs is the number of pixels in the object that are correctly segmented, N g is the number of pixels in the ground truth object, N s is the number of pixels in the segmented object.
In Table 4, we list the Dice metric numbers obtained in PASCAL-VOC2012 dataset to evaluate the quality of our segmentation results. The data shows that while the proposed method can achieve similar results as the ALM method, the runtime had been reduced to approximately one-sixth.  Table 3. The number of iterations to reach convergence (iters) and the total CPU time in seconds (CPU(s)) for different images of different sizes by using our algorithm and the ZTC algorithm in the segmentation problem.

Concluding remarks
In this paper, we used a dual method to solve the Euler's elastica regularizer for image denoising and segmentation. Our method can efficiently reduce the number of parameters, and formulate a more concise algorithm. There are two main contributions. Firstly, by introducing the dual operators, the optimization problem can be divided into simpler sub-problems, and the Chan-Vese model with elastica can be solved more easily. Secondly, our proposed algorithm can effectively reduce the number of parameters, to reduce the dependence on of parameter tuning. Numerical experiments show that compared with the ALM method our algorithm can obtain similar experimental results with less running time.