Augmented Lagrangian Method for A Mean Curvature Based Image Denoising Model

High order derivative information has been widely used in developing variational models in image processing to accomplish more advanced tasks. However, it is a nontrivial issue to construct eﬃcient numerical algorithms to deal with the minimization of these variational models due to the associated high order Euler-Lagrange equations. In this paper, we propose an eﬃcient numerical method for a mean curvature based image denoising model using the augmented Lagrangian method. We describe how to construct an augmented Lagrangian functional for the model and detail the procedures of ﬁnding the related saddle-points of the functional. We present numerical experiments to illustrate the eﬀectiveness and eﬃciency of the proposed numerical method, and show a few important features of the image denoising model such as keeping corners and image contrast. Moreover, a comparison with the gradient descent method further demonstrates the eﬃciency of the proposed augmented Lagrangian method.


Introduction
Image denoising provides a crucial step towards the capturing image signals from noisy images with lots of applications in medical imaging and video monitoring etc. During the last three decades, numerous methods have been applied to deal with this problem, such as variational methods, partial differential equations, statistical methods and so on [3,4,8,23,24,25,27,28,29,33]. Among these variational method is one of the most popular methods. A classical variational model on image denoising was developed by Rudin, Osher, and Fatemi (ROF) [29], in which achieving a clean image from a noisy one amounts to the minimization of the following functional: where f : Ω → R is a given noisy image defined on Ω (always a rectangle in R 2 ) and λ > 0 is a positive constant. The ROF model has proved to be very efficient in the noise removal 1 tasks, and it has also been widely used in other tasks of image processing such as image deblurring, image inpainting etc. Besides these, the idea of the total variation has broadly been utilized in many other disciplines.
Even though the ROF model is very powerful in removing noise while preserving edges and contours of objects, it also has some unfavorable properties. For instance, it yields the staircase effect, smears object corners and leads to the loss of image contrasts. To remedy these drawbacks, quite a few high order variational models have been proposed [1,2,21,22,10,19,38]. For instance, in [1], the authors used the Euler's elastica of level curves of a smooth function as a regularizer and proposed the following energy functional: In [19], Lysaker et al. directly incorporated second order derivative information into the image denoising process and their energy functional reads: In [38], two of the authors proposed a variational model that uses the mean curvature of the induced image surface (x, y, f (x, y)) to remove noise. Specifically, the model employs the total variation of mean curvature as a regularized term and the functional is written as follows: where λ is a tuning parameter. The term ∇ · ( ∇u √ 1+|∇u| 2 ) is the mean curvature of the surface φ(x, y, z) = u(x, y) − z = 0. The model tries to fit the given noisy image surface (x, y, f (x, y)) with a surface (x, y, u(x, y)) that bears small amount of mean curvature. The model can sweep noise while keeping object edges, and it also ameliorates the staircase effect. More importantly, as discussed in [38], the model is also capable of preserving image contrasts as well as object corners.
While all the models can effectively accomplish the noise removal task, the associated functionals are not easy to minimize. Since they involve high order derivatives, the related Euler-Lagrange equations are often fourth-order, which raises a nontrivial issue of developing effective and efficient numerical algorithms to solve them. Indeed, there are lots of high order model existing in the literature of image processing [5,8,9,12,13,14,15,26,32,36,37]. Therefore, there is a pressing need to develop efficient numerical methods for these models.
In this paper, we plan to develop a fast algorithm to minimize the functional (4). We employ the idea of augmented Lagrangian methods that have already been successfully utilized in the minimization of the ROF model [34] and Euler's elastica based functionals [31]. In these works, the augmented Lagrangian methods achieve much higher speed than other numerical methods, such as the gradient descent method, primal-dual methods. The reason of how the method is so efficient mainly lies in the following fact: the augmented Lagrangian methods proposed in [34,31] decomposed the original nontrivial minimization problem to be a few simple ones, some of which can be solved using the Fast Fourier transformation (FFT), while the other ones can be quickly dealt with due to their existing explicit solutions. Therefore, the construction of an efficient augmented Lagrangian method for the minimization of a given functional depends on whether one can break down the original functional into simple ones. In the functional (4), note that the mean curvature ∇ · ( ∇u √ 1+|∇u| 2 ) is not homogeneous in the variable u, how to introduce a new variable that can commute freely with ∇u is a key problem for developing the associated augmented Lagrangian method of the functional (4). To solve this problem, just as what we regard a 2D image function as a surface in R 3 , we introduce a new variable p =< ∇u, 1 >, instead of p = ∇u, to connect with ∇u. In this way, we can apply the augmented Lagrangian method to the minimization of the functional (4).
The rest of this paper is organized as follows. In Section 2, we will first review the mean curvature based image denoising model, recall the augmented Lagrangian method for Euler's elastica based variational models [31], and then will detail how to apply the augmented Lagrangian method to the mean curvature denoising model. In Section 3, we will discuss how to carry out the numerical implementation. In Section 4, we present numerical experiments to illustrate the effectiveness and efficiency of the proposed numerical method, which is followed by a conclusion in Section 5.

Augmented Lagrangian Method
In this section, we first review the mean curvature image denoising model and the augmented Lagrangian method to the Euler's elastica functional, then detail how to treat the mean curvature based image denoising model using the augmented Lagrangian method.

The Mean Curvature Based Image Denoising Model
In [38], two of the authors proposed a variational model that uses the mean curvature of image surface for image denoising. For simplicity, we will call it as MC model in the following context. The functional can be written as follows: where λ is a tuning parameter and the term κ u = ∇ · ( ∇u √ 1+|∇u| 2 ) is the mean curvature of the image surface φ(x, y, z) = u(x, y) − z = 0. This model aims for obtaining a surface that is close to the original induced surface (x, y, f (x, y)) and has a small amount of mean curvature. Therefore, noise part of an image can be removed to reduce the amount of mean curvature.
The motivation of the MC model is to remedy some unfavorable features of the classical ROF model [29] such as the staircase effect. To avoid this effect, as discussed above, numerous variational models that use high order derivative information have already developed [1,2,21,22,10,19]. The MC model regards a given 2D image as a surface in 3D space, and uses the total variation of mean curvature of the surface as the regularizer. Therefore, the model suppresses the oscillation of mean curvature and prefers small mean curvature that lead to piecewisely smooth image surface. This effectively ameliorates the staircase effect.
Besides removing noise, the MC model is also capable of preserving image contrasts. In fact, as discussed in [38], for a region R ⊂ Ω with a C 2 boundary, one gets Ω |κ hχ R | = H 1 (R) for the image hχ R , where χ R and H 1 (R) are the characteristic function and the perimeter of the region R, respectively, and h depicts the height. Note that the regularizer is independent of the height h, which results in the preservation of image contrasts. It demonstrates an important difference between this model and the ROF model whose regularizer Ω |∇(hχ R )| equals to hH 1 (R) that relies on the height h.
Moreover, the MC model is capable of preserving object corners. This suggests that the MC model can keep lots of visually meaningful clues. It is well-known that object corners will be smeared when the ROF model is applied.
Even though the MC model possesses the above merits, it is a challenging issue to develop numerical algorithms to minimize its functional as it involves high order derivatives and the non-differentiable regularization term. In fact, due to the non-differentiability of the function |x|, one cannot obtain the associated Euler-Lagrangian equation. Instead, as in [26], one may consider the following smoothed version of the functional (5) with and the associated Euler-Lagrangian equation can then be written as follows [38]: where I, P : , and the associated gradient descent equation of the functional (6) can be given as follows: with time t being an evolution parameter.

The Augmented Lagrangian Method to the Euler's Elastica functional
In [31], the augmented Lagrangian method was successfully applied to the Euler's elastica based image denoising model that was supposed to remove Gaussian white noise [1,2]. The functional reads As discussed in [31], one may convert the minimization of this functional as a constrained minimization problem by introducing two new variables p = ∇u and n = ∇u/|∇u|, and 4 the corresponding augmented Lagrangian functional can be written as follows: L(u, p, n; λ 1 , However, instead of dealing with this functional, in [31], the authors considered a more relaxed and easily handled functional that is given as follows: where R = {m ∈ L 2 (Ω) : |m| ≤ 1 a.e. in Ω} and δ R (·) is the characteristic function on R and can be expressed as follows: δ R (m) = 0, m ∈ R; +∞, otherwise.
One then needs to find its saddle points as these points correspond to the minimizers of the original functional (10). To this end, as in [31], one just needs to solve the associated subproblems of each variable by fixing other ones iteratively with new updated Lagrangian multipliers λ 1 , λ 2 , λ 3 , and λ 4 . There are totally four subproblems related to the functional, each of which involves only one variable of u, m, p, and n. Specifically, the equation for the variable u is a linear equation that can be efficiently solved using FFT, and the one for n can be locally solved. Moreover, the other two subproblems have explicit minimizers using some shrinkages. Therefore, the new augmented Lagrangian functional (12) will be solved very efficiently.
From the above discussion, whether an augmented Lagrangian method can be efficiently employed is highly determined by the original functional. Basically, a functional could be easily converted to be a constrained problem by introducing new variables. However, the searching for saddle points of the associated augmented Lagrangian functional may require quite a lot of effort. The success of the above augmented Lagrangian method is benefited by the simpleness of the associated subproblems.

Description of the Augmented Lagrangian Method to the Mean Curvature Based Image Denoising Model
In this section, we discuss how a similar augmented Lagrangian method will be applied to the MC model (5).
To develop an augmented Lagrangian functional for the mean curvature based image denoising model, we may naturally introduce three new variables q, w, and p and consider the following constrained minimization problem: We then write down the associated augmented Lagrangian functional: However, due to the term p/ 1 + |p| 2 , this functional is also not easy to handle. This differs from the Euler's elastica functional where p = ∇u and n = ∇u/|∇u| can be nicely and mutually converted. Therefore, a key issue of applying the augmented Lagrangian method to functionals involving the term 1 + |∇u| 2 or ∇u/ 1 + |∇u| 2 is to find a way that can easily solve p from 1 + |p| 2 . In fact, the term 1 + |∇u| 2 also appears in other occasions such as the surface area based image denoising model [17]: Before stating our idea of dealing with this issue, let's recall the idea of introducing the mean curvature denoising model. In this model, a given 2D image f (x) is regarded as a surface (x, y, f (x, y)) in R 3 . One thus considers the surface φ(x, y, z) = u(x, y) − z = 0 and the mean curvature κ = ∇ · (∇φ/|∇φ|) = ∇ · (< ∇u, −1 > /| < ∇u, −1 > |). Note that one introduces two variables p = ∇u and n = ∇u/|∇u| to tackle the Euler's elastica for its curvature term κ = ∇ · (∇u/|∇u|). This gives us a hint to treat the curvature term in our case, that is, we may introduce a variable p =< ∇u, −1 > instead of p = ∇u and n =< ∇u, −1 > /| < ∇u, −1 > | accordingly.
As there is no difference if p =< ∇u, 1 > is used, we then propose the following constrained minimization problem that is equivalent to the minimization of the functional (5): 6 and the associated augmented Lagrangian functional: where n, m, p ∈ R 3 , λ 1 , λ 3 ∈ R, λ 2 , λ 4 ∈ R 3 . All these terms are similar to the ones discussed in [31]. For the sake of the completeness of presentation, we want to emphasize important points of the above augmented Lagrangian functional as follows.
The introduction of the variable m is to relax the variable n that is supposed to connect with the variable p in terms of n = p/|p|, and the variable m is required to lie in the set R so that the term |p| − p · m is always non-negative. As discussed in [31], the benefit of this non-negativeness is that the L 2 penalization is unnecessary and we just use |p| − p · m as a penalization.
As the saddle points of the augmented Lagrangian functional (17) corresponds to the minimizers of the constrained minimization problem (16), one just needs to find the saddle points of (17). To this end, as in [31], we apply an iterative algorithm. Specifically, for each variable in (17), we fix all the other variables and seek a critical point of the induced functional to update this variable. Once all the variables are updated, the Lagrangian multipliers will also be advanced. Then we repeat all the process until the variables are convergent. The algorithm is summarized in Table 1 and Table 2.
Therefore, we need to consider the following subproblems and obtain critical points for each of them.
4. Measure the relative residuals and stop the iteration if they are small than a threshold ǫ r . Table 2: Alternating minimization method to solve the subproblems.
2. For l = 1, · · ·, L and fixed Lagrangian multiplier , solve the following subproblems : In what follows, we will discuss the minimizers of all the above functionals and the update of the Lagrangian Multipliers.
Similarly as what discussed in [31], the minimizers of the functionals ε 2 (q), ε 3 (p), and ε 5 (m) can be expressed explicitly, while the minimizers of the functionals ε 2 (u) and ε 3 (n) are determined by the associated Euler-Lagrange equations. For the sake of completeness of presentation, we provide the details here.
To obtain the minimizers of the functionals ε 2 (q) and ε 3 (p), we need the following lemmas [31]. Lemma 1. Let x 0 be a given scalar or vector, µ, r be two parameters with r > 0, then the minimizer of the following functional [18]) Let m 0 be a given vector, then the minimization problem

Lemma 2. (Lions and Mercier
has an explicit solution Note that where c 1 is a constant independent of the variable q, and where c 2 is a constant independent of the variable p, based on Lemma 1, we obtain the minimizers of ε 2 (q) and ε 3 (p) as follows: 9 As where c 3 is a constant independent of the variable m, based on Lemma 2, we get As for the functionals ε 1 (u), ε 4 (n), we need to get the associated Euler-Lagrange equations. Standard procedures lead to the following equations of u and n: and where p =< p 1 , p 2 , p 3 >, m =< m 1 , m 2 , m 3 >, n =< n 1 , n 2 , n 3 > and the Lagrangian multiplier λ 2 =< λ 21 , λ 22 , λ 23 >, λ 4 =< λ 41 , λ 42 , λ 43 >. To update the variables u and n, one needs to solve these Euler-Lagrangian equations.

Discussion of the Smoothed Mean Curvature Based Image Denoising Model
As discussed previously, in [38], due to the non-differentiability of the MC model (5), the gradient descent equation of the smoothed version of the MC model (6) had to be employed. However, by using the augmented Lagrangian method, the difficulty can be easily overcome [34,31]. This is one appealing feature of the augmented Lagrangian methods. In this work, to fairly get a comparison of efficiency for the proposed augmented Lagrange method and the gradient descent method, we also consider the functional (6). In fact, to minimize this functional, we may use the same augmented Lagrangian functional (17) except that the term λ |q| that will be replaced by λ Φ ǫ (q) with Φ ǫ defined as follows: where ǫ > 0 is a parameter. This modification only leads to the change of the subproblem related to the variable q, and it becomes Therefore, the subsequent problem is how to determine its minimizer. To this end, we propose the following lemma.
Lemma 3. Let q 0 be a given scalar, µ, r be two parameters with r > 0, then the following functional with Φ ǫ being defined as (36) attains its minimizer at q 1 ≥ ǫ and q 2 ≥ ǫ, ǫ, q 1 ≥ ǫ and q 2 < ǫ, q 1 , q 1 < ǫ and q 2 ≤ ǫ, argmin{f (q 1 ), f (q 2 )}, q 1 < ǫ and q 2 > ǫ. (38) where f (q) = µΦ ǫ (q) + r 2 (q − q 0 ) 2 , and q 1 , q 2 are given as follows: Proof. Since no derivatives involved in the functional, one just needs to find the minimum point of the following scalar function: Note that f is piecewisely defined on the intervals (−∞, −ǫ), [−ǫ, ǫ], and (ǫ, +∞), to find the minimum, one just needs to calculate the critical points for each interval. As f is an even function, only the case with q 0 > 0 is considered below. The critical points of f defined in [−ǫ, ǫ] and (ǫ, +∞) are given as respectively. As f is upward convex in each of these intervals, f attains the minimum at the critical point if it is inside the interval where f is defined. Therefore, if q 1 ≥ ǫ, the minimum is either at ǫ or at q 2 , which can be determined by whether q 2 is inside (ǫ, +∞); if q 1 < ǫ, the minimum may be attained either at q 1 or q 2 . In summary, we get the minimum point of f as stated in this lemma.

Numerical Implementation
In this section, we present the details of how to solve the equations (33) and (34) and update the variables q, p and m for each iteration. Since the numerics are almost the same as what discussed in [31], we here only present the key points.
Before discussing the discretizations of the above equations, we need to note that the regularizer of the MC model, |κ u | = |∇ · ( ∇u √ 1+|∇u| 2 )| or Φ ǫ (κ u ), is not homogeneous in u. This is one of the most important features that differentiate it from other image denoising ones such as the ROF model [29] and the Euler's elastica based model [1,2]. Therefore, one shall introduce a spatial mesh size for the discretization of the terms involving derivatives.
Numerically, the spatial mesh size plays an important role in the performance of the model, as it determines the magnitude of gradient. Note that any given image is only defined on a finite number of grids, therefore, once the spatial mesh size is small enough, any tiny variation of u in a neighborhood of a grid will lead to large κ u at that point and may be preserved as a jump by the model. This surely fails to remove the noise but also suggests that a relative small spatial mesh size helps to preserve fine structures of images. On the other hand, if the spatial mesh size is large enough, and as the image intensity only ranges in [0, 255], then even visually salient jumps will still lead to a small κ u and be removed.
We then define the discrete backward and forward differential operators with periodic boundary condition and the spatial mesh size h as follows: and the central difference operators and the gradient operators are defined accordingly as Let's first discuss how to solve the equation (33). As periodic condition is imposed, we may use FFT to solve this linear equation. We employ the following discretization: We then apply the discrete Fourier transform F for both sides. Note that where z i = 2πi/M, i = 1, · · ·, M , z j = 2πj/N, j = 1, · · ·, N , we get where j). Therefore, once Fu is calculated, u can be obtained using the discrete inverse Fourier transform.
We also use FFT to solve the equations (34). To this end, we first discretize the first two equations of (34 ) as follows: then apply the discrete Fourier transform to both sides and get the following 2 × 2 system for each grid (i, j): As the determinant of the above 2×2 matrix equals to r 2 4 −2r 4 r 3 (cos 2πi/M +cos 2πj/N −2), which is always positive when r 4 > 0, we can easily solve the above system to get Fn 1 and Fn 2 , then apply the discrete inverse Fourier transform to them, and the real parts of the inverse Fourier give us the new n 1 and n 2 . The third component n 3 of n can be calculated directly.
In what follows, we discuss the update of the variables q, p, m as well as the Lagrangian multipliers.
As q is scalar defined on the grid points (i, j), based on the formulation (30), one gets

13
As for the variable p, we first calculate the three components of p, and then calculate the length of p and the updated p. Specifically, and based on the formulation (31) Similarly, we calculate and get the new m(i, j) using the formulation (32). Moreover, based on the formulations (35), we may update all the Lagrangian multipliers: j)).

Numerical Experiments
In this section, we apply the proposed augmented Lagrangian method (ALM) for both the original MC model (5) and the smoothed MC model (6), and then compare the algorithm with the gradient descent method (GDM) in their efficiency for the smoothed MC model. To check whether the ALM converges to some saddle point, we monitor the following relative residuals as in [31]: with In all the following numerical experiments, we use the relative residuals (45) as the stopping criterion, that is, given a small threshold ǫ r , once R k i < ǫ r for i = 1, ..., 4 and for some k, the iteration process will be terminated. To check the convergence of the iteration process, as in [31], we also check the relative errors of Lagrange multipliers: and the relative error of the solution u k Besides all these quantities, we also consider how the energy (5) is evolving during the iteration by tracking the amount E(u k ). For the presentation purpose, all the above quantities are shown in log-scale in the following results. Moreover, to illustrate what signals are removed as noise, we also present the associated residual image f − u, besides the given noisy image f and the cleaned one u.
We first consider some real images. In Fig 1, we present the results for the images "Cameraman" and "Peppers", including the original noisy images, cleaned ones and the residual images. The cleaned images demonstrate that noise part and small scale signals have been removed, which can be observed in the related residual images. One can easily find that the cleaned images are composed of piecewisely smooth surface patches, that is, they are free from the stair case effect that is a typical unfavorable feature of the well-known ROF model. In Fig 2, the plots of the relative residuals (45), relative errors in Lagrange multipliers (46), relative error in u k (47), and the energy versus iteration are listed for these two examples. These series of plots illustrate the convergence of the iteration process using the proposed Augmented Lagrange method, and also show that the process leads to some saddle points of the constrained minimization problems.
We then apply the proposed Augmented Lagrange method for some synthetic images. In Fig 3, the results of two synthetic images are presented. These two images have some special visual clues. For instance, the image "lattice" has corners and sharp edges, and the image "rings" consists of a few surface patches with both flat patches and a curved one inside the image. From the two cleaned images, one can check that the noise has been effectively removed and the corresponding residual images nearly contain no meaningful signals for both examples. A careful check of the cleaned lattice image shows that the corners are preserved very well, which is an important feature of the mean curvature model [38]. Moreover, in the plots (a4) and (b4), we present the slices of the original noise free images, noisy images and the cleaned ones. These slices provide another perspective to check the denoising results. In the plot (a4), we can see that the cleaned slice curve (blue) almost overlaps with the original noise free slice curve (red), which demonstrates the preservation of image contrast, which is another feature of the mean curvature model. The plot (b4) shows that the original noise free "rings" has a smooth arch in its center, which is also almost fully restored by the model.
We have applied the proposed ALM for real and synthetic images. To explicit how efficiency the proposed method is, in Table 3, we present the cpu time for all the above Table 3: The presentation of the sizes, the SNRs, the numbers of total outer iteration, and the computational time for the images in Figures 1 and 3.
image size SNR Number of iterations time (sec) Fig. 1-(a2)  As mentioned in the introduction, a few other higher order image denoising models have been existed in the literature such as Euler's elastica image denoising model [1] and the Lysaker et al.'s model [19]. In Figure 4, we apply the MC model, Euler's elastica model and the ROF model for a synthetic image with sharp edge and corners. The results show that the MC model keeps corners better than the other two models. This can be observed in the related residual images where some corner signals exist explicitly for the two models. Moreover, from the cleaned images, one can see that the MC model also preserves image contrast better than the ROF model. In fact, the result obtained by the ROF shows less darkness of the object than that of the MC model.
The above three examples demonstrate that the proposed ALM can efficiently accomplish the task of sweeping noise and keeping corners and image contrasts that are assumed by the MC model.
In the next experiment, we demonstrate an important feature of the MC model, that is, the clean image function u depends discontinuously on the regularization parameter λ. To see this point clearly, we compare the results using the ROF model and the MC model respectively in Figure 5. In this experiment, we consider a noise-free synthetic image with a bunch of squares of different sizes. We first list the results of the MC model for two different regularization parameters. One can see that when λ = 2.5 × 10 4 , the squares with the smallest size are swept as noise and show up in the residual image, while the other squares are well preserved in the cleaned image. When we increase the parameter to be λ = 4.0×10 4 , then the squares of larger size will be removed as shown in the residual image. This phenomenon demonstrates that the results of the MC model do not continuously depend on the regularization parameter λ. In contrast, the ROF model depends on its regularization parameter continuously. In fact, for either a relatively small parameter λ = 10.0 and a relatively larger one λ = 100.0, all the squares, regardless of the sizes, are partially swept as noise, and the removed part depends on the square size.
This specific feature of the MC model suggests a data-driven scale selection approach, that is, as the regularization parameter λ gradually increases, objects of small sizes will be removed from the given image, while large size objects will be almost completely preserved. This particular property is also possessed by the TV-L1 model [6] with lots of applications [7,35]. However, the MC model is able to keep corners, while the TV-L1 model cannot, which suggests that the MC model might be more suited than the TV-L1 model to applications for which the preservation of geometry is crucial.
In the last experiment, we compare the efficiency of the proposed ALM and the GDM by applying them to the smoothed MC model (6). To make the comparison fairly, for each method, we stop the iteration process once the following L 2 error in u k with satisfies: where ǫ e = 5.5 × 10 −3 , u exact denotes the convergent state of u k . Note that the MC model is not convex, the two algorithms may converge to two different local minimizers. Therefore, for each algorithm, we calculate its own u exact . For instance, for ALM, we get the u exact using a far large number of iterations; while for GDM, besides the huge iteration number, we choose a much small time step. By using the quantity 1 |Ω| u k −u exact L 2 , one can check which algorithm can converge to u exact more quickly. In Figure 6, we list the results for both algorithms. The two cleaned images show almost no difference, however, the energy associated with ALM drops dramatically within one hundred iterations. From the efficiency comparison Table 4, we also see that ALM takes much fewer iterations and less cost time than GDM, which demonstrate that ALM is much more efficient than GDM.

Conclusion
Recently, the augmented Lagrangian method has been successfully used to minimize the classical ROF functional [34] and Euler's elastica based functionals [31]. In this paper, we apply this method to deal with a mean curvature based image denoising model, which is nontrivial to develop an efficient algorithm due to higher order derivative information and non-differentiability of the model. The numerical experiments demonstrate that the proposed augmented Lagrangian method is much faster than the gradient descent method in converging to minimizers of the mean curvature based functional. (a3) (b3) Figure 1: The denoising results for the images "Cameraman" and "peppers". The noisy, cleaned, and residual images are listed from the first row to the third row respectively. For both examples, we set the spatial step size h = 5.0, choose r 1 = 40, r 2 = 40, r 3 = 1 × 10 5 , and the remaining parameters r 4 = 1.5 × 10 5 and ǫ r = 6 × 10 −3 for "Cameraman" while r 4 = 10 5 and ǫ r = 6 × 10 −3 for "peppers".  Figure 2: The plots of relative residuals, relative errors in Lagrange multipliers, relative error in u k and energy for the previous examples "Cameraman" and "peppers". The left column lists the plots for "Cameraman" while the right one for "peppers". The plots demonstrate the converges of the iteration process, showing that some saddle points will be achieved. (a4) (b4) Figure 3: The denosing results for two synthetic images "lattices" and "rings". The noisy, cleaned, residual images are listed from the first row to the third row respectively. In the fourth row, a slice of the original noise free image, noisy and the cleaned images for each example is presented. For the "lattices" example, we choose the parameters as follows: h = 1.0, r 1 = 20, r 2 = 20, r 3 = 1 × 10 4 , r 4 = 1 × 10 5 and ǫ r = 5 × 10 −3 , while for the "rings" example, h = 5.0, r 1 = 20, r 2 = 20, r 3 = 1 × 10 5 , r 4 = 2 × 10 5 and ǫ r = 8 × 10 −3 . 22 by MC model by Euler's elastica model by ROF model Figure 4: The comparison of three models: the mean curvature (MC) model, the Euler's elastica model and the ROF model. The first row lists the noise free image and the noisy one. The second row presents the results by applying the above three models respectively, while the associated residual images are shown in the third row. The results illustrate that the MC model keeps corners better than the other two models, which can be observed in the related residual images where some corner signals exist explicitly for the two models. Moreover, the cleaned image by the ROF shows less darkness of the object than that of the MC model, which suggests that the MC model preserves image contrast better than the ROF model.  Figure 6: The efficiency comparison of ALM and GDM for the smoothed MC model. The first row lists the two results by both algorithms; the second row shows how the energy is evolving for the two algorithms; and the third row presents the L 2 -norm of the difference between the calculated solution and the exact solution, which is calculated with high accuracy. The common parameters used for these two methods are: λ = 5 × 10 3 , h = 5.0, and ǫ = 0.2. The parameters for ALM are r 1 = 50, r 2 = 100, r 3 = 5 × 10 4 , and r 4 = 5 × 10 4 . The time step size for GDM is dt = 2 × 10 −4 .