Augmented Lagrangian method for an Euler’s elastica based segmentation model that promotes convex contours

In this paper, we propose an image segmentation model where an L1 variant of the Euler’s elastica energy is used as boundary regularization. An interesting feature of this model lies in its preference for convex segmentation contours. However, due to the high order and non-differentiability of Euler’s elastica energy, it is nontrivial to minimize the associated functional. As in recent work on the ordinary L2-Euler’s elastica model in imaging, we propose using an augmented Lagrangian method to tackle the minimization problem. Specifically, we design a novel augmented Lagrangian functional that deals with the mean curvature term differently than in previous works. The new treatment reduces the number of Lagrange multipliers employed, and more importantly, it helps represent the curvature more effectively and faithfully. Numerical experiments validate the efficiency of the proposed augmented Lagrangian method and also demonstrate new features of this particular segmentation model, such as shape driven and data driven properties.


Introduction
Image segmentation is one of the fundamental topics in the fields of image processing and computer vision.The aim is to partition an image region into sub-regions in order to capture meaningful objects.In the literature, numerous approaches have been proposed for this problem.For instance, in the classical snake and active contour model by Kass, Witkin, and Terzopoulos [22], segmentation contours are driven to object boundaries by internal forces, such as regularity, and external ones, like the intensity gradient.They proposed minimizing the following functional: where f : Ω → R 2 represents a given image, C(s) : [0, 1] → R 2 is a parameterized curve, and α, β are positive tuning parameters.The first two terms impose regularity assumptions on the contour while the third one imposes restrictions induced by the given image.In fact, as f takes on large gradients along object boundaries, the functional E(C) will take on a small value when the active contour C resides on these boundaries.Later on, in [7], Caselles, Kimmel, and Sapiro developed a variational segmentation model called geodesic active contours, where one needs to minimize the functional: over all C(s) ∈ S = {C : [0, 1] → R 2 : C(s) is piecevise C 1 , and C(0) = C(1)}.The function g o : [0, +∞) → [0, +∞) is an edge detection function satisfying g o (0) = 1, lim r→+∞ g o (r) = 0, and which is also strictly decreasing.The edge detection function takes on values close to 1 in regions with homogenous grey intensity while being close to zero at boundaries with sharp intensity gradients.Therefore, the above functional attains its minimum value when the active contour C(s) lies along edges in the image.Another classical segmentation model was proposed by Mumford and Shah [30].In this approach, the segmentation problem amounts to finding a piecewise smooth function which approximates the given image, where also the discontinuity set is restricted to be smooth.They proposed minimizing the following functional: with respect to both the clean image function u defined on Ω and the discontinuity set Γ ⊂ Ω. λ, µ > 0 are tuning parameters.The methodology has brought forth plenty of variational models in many topics of image processing, including segmentation, denoising, inpainting, etc.An interesting offspring of Mumford-Shah's model is Chan-Vese' model [8], where the discontinuity set Γ is restricted to be closed curves, leading to a segmentation of the image into two or several regions.Moreover, a nice feature of Chan-Vese's model lies in its treatment of the boundary Γ via level set functions [32] that can handle curves with topological changes easily.The above mentioned segmentation models are all grey intensity based models, that is, the segmentation result mainly relies on grey intensity values of the images.However, due to the complexity of real images, meaningful objects might be occluded by other ones, or some parts of them may not be distinguished from the background by the contrast.For instance, in medical image analysis, target organs may be blended with other ones, and some parts of them may be occluded by other organs or even missing due to the imaging conditions.Therefore, those grey intensity based segmentation models might not be well suited for segmenting the target objects from such images.To overcome this issue, one possible way is to incorporate prior knowledge about the shapes of the target objects into the segmentation process.A lot of work have also focused on this research topic -segmentation using prior shape knowledge, see e.g.[12,13,14,23].Another possible way is to impose constraints on the segmentation contours to restore the boundaries that are missing or not well defined by the grey intensity gradient of the images.Of course, those well-known models discussed above all contain regularity terms on the segmentation contours.However, the main purpose of these terms is to improve the smoothness of the contours.They add no other geometric preference on the contours intentionally.In this work, we consider a new variant of Chan-Vese's model by using the L 1 -Euler's elastica energy as regularization on the segmentation contours.As will be discussed later, the most interesting feature of this new model lies in its capability of producing segmentations with convex contours.Moreover, the new segmentation model also possesses other new features -data driven and shape driven properties, the later of which has never been discussed in the literature, to the best of our knowledge.
Before introducing the L 1 -Euler's elastica energy, recall first the standard Euler's elastica energy of the curve Γ: where κ represents the curvature of the curve and a, b > 0 are two parameters.In the field of mathematical imaging, Euler's elastica was first introduced by Nitzberg, Mumford, and Shiota in their famous work on segmentation with depth [31].Since then, it has been adapted to many fundamental problems in imaging, such as image inpainting [11,26], image restoration [1,2,39], and image segmentation [42,45].Notice that this elasticity energy contains a quadratic dependence on curvature.In fact, in [31] the authors considered another variant of Euler's elastica that linearly depends on the magnitude of curvature, that is, Γ (a + b|κ|)ds, based on the requirement of preserving corners in the segmentation contours.For the brevity of presentation, these two Euler's elastica terms will be called the L 2 -and L 1 -Euler's elastica terms respectively in the later sections.
Recently, a lot of research have focused on the development of fast and reliable numerical methods for minimizing curvature based functionals, such as the multigrid algorithm [5], the homotopy method [41], augmented Lagrangian method (ALM) based algorithms [39,40,45] and graph cut based algorithms [4,18].Another interesting line of research is the design of convex relaxation approaches [6,29,37,38].Their main advantage is the ability to avoid getting stuck in local minimums, but they also have several limitations: The computational cost and memory requirements are much higher due to the higher dimensional formulations of the problems; The curve must be represented discretely by a finite number of line segments; The convex relaxations only correspond approximately to the original problems, and in particular the tightest convex relaxation for the Euler's elastica model without additional boundary length regularization reduces to constantly zero [6].
In this work, inspired by our previous work [45], we focus on the exact non-convex formulation and intend to construct a novel ALM based algorithm with fewer Lagrange multipliers, which markedly reduces the effort of choosing appropriate parameters for obtaining minimizers of the proposed model.Moreover, by employing fewer auxiliary variables, more direct connections among those variables are ensured, which leads to a more accurate discrete representation of the curvature of the segmentation contours.The augmented Lagrangian method was first used to get fast numerical schemes for nonlinear image processing models in [40], including the total variation based Rudin-Osher-Fatemi (ROF) model.Extensions of the augmented Lagrangian method for the ROF model to a spatially continuous setting were proposed in [21], and an interesting discussion about convergence properties were given.The fast algorithms were extended to L 2 -Euler's elastica and surface curvature models in [39,44,45].For these non-linear higher order problems, there are many auxiliary variables to compute and constraints for the Lagrangian methods to handle.One of the contributions of this work is to reduce the number of auxiliary variables and constraints.We note that Myllykoski et al. [28] recently also considered an approach to reduce the number of Lagrange multipliers for an ALM based method for the surface curvature based image denoising model in [43].
As a summary, the new contribution of the current work lies in two aspects: • We propose a curvature based segmentation model that promotes segmentations with convex contours; • We develop an augmented Lagrangian algorithm for the new segmentation model, where curvature is treated in a novel and more accurate way, using fewer variables and parameters than in related work.
The rest of this paper is organized as follows.In the next section, we discuss the new segmentation model, where the L 1 -Euler's elastica is used as regularization of the segmentation contours, and also present the specific application to segmentation of objects with convex contours.In Section 3, we discuss a natural extension of our previous work [45] for this new model and address several problematic issues, such as choosing appropriate parameters.After than, a novel ALM algorithm is proposed to ameliorate those issues.We then discuss its numerical implementation in Section 4, which is followed by numerical experiments that validate the features of the new model and the efficiency of the proposed ALM algorithm in Section 5. Section 6 is devoted to conclusions.

L 1 -Euler's elastica extension of Chan-Vese's model
In [45], we considered an extension of Chan-Vese's segmentation model by using the L 2 -Euler's elastica as a new regularizer on the segmentation contours instead of the standard curve length regularization.This simple modification introduces new features of the resulting model: 1) the model is able to connect the boundary around missing parts of objects; 2) it can also capture objects of relatively large size while omitting those of small sizes; 3) it is more suitable than the original Chan-Vese's model for keeping elongated structures.These features have particular applications in medical imaging.For instance, MRI images often contain the mixture of organs or tissues of different sizes.By applying this model, it is possible to locate the boundaries of organs of considerable sizes while ignoring tiny tissues.
In this work, we intend to study another modification of Chan-Vese's segmentation model, that is, instead of the L 2 -Euler's elastica, we utilize the L 1 -Euler's elastica as regularization on the segmentation contour.With the aid of a level set function [32], we may express this new variational model (called ECV -L 1 segmentation model) as follows: where f : Ω → R is a given image, φ is a level set function whose zero level set locates object boundaries, and a, b > 0 are parameters.The term ∇ • (∇φ/|∇φ|) denotes the mean curvature of the level curves of φ, and the second part of the above functional just represents the L 1 -Euler's elastica of the zero level set of φ.The most remarkable feature of this new modified segmentation model is that it promotes convex contours provided the parameter b is sufficiently large.It may therefore be used for the particular task of segmenting objects that are known to have a convex shape.This specific feature is supported by the following theorem in differential geometry [16]: Theorem 1.Let Γ be any closed piecewise smooth curve in R 2 and κ its curvature, then Γ |κ|ds ≥ 2π and equality holds when Γ is a convex curve.
We remind that a convex curve is a closed curve that is the boundary of a convex set.As this integral along convex curves is always 2π, one still needs the length term to prohibit excessively large contours in the segmentation.
A nice feature of using a level set function in the above model is that since multiple curves can be represented by a single level set function φ, one may impose convex shape priors on several objects simultaneously.Moreover, if the parameter b is chosen large, only those convex contours that encompass relatively large regions will be preserved in the final segmentation.
In contrast to the ECV -L 2 segmentation model introduced in [45], the ECV -L 1 segmentation model discussed in this paper assumes two new features: first, it prefers convex boundaries when the parameter b is chosen large; second, the model is capable of preserving sharp corners, since the L 1 -Euler's elastica linearly depends on the magnitude of curvature, while the ECV -L 2 segmentation model necessarily smears out corners.The corner preservation property of the L 1 -Euler's elastica energy was first pointed out in [31].
Just as in many other higher order variational models [1,2,31,43], the numerical treatment of this new model is also nontrivial.The difficulty is mainly due to the presence of the curvature term, which makes the model high order, nonlinear, and non-differentiable.Recently, quite a few research works have focused on the development of efficient algorithms for these models, such as the multigrid algorithms [5], the homotopy method [41], and augmented Lagrangian methods (ALMs) [40,39,45].In this work, we confine ourself to ALM in order to develop an efficient and reliable fast algorithm for the functional (Eq.5), which will be detailed in the following sections.

The development of a novel augmented Lagrangian method
In this section, we first consider a natural extension of our previous work [45] for the new segmentation model and point out its problematic issues, and then develop a novel ALM for the minimization of the original model (Eq.5).

Issues arising from the previously proposed ALM for the model
Recently, ALMs have proven to be very useful in designing fast algorithms for variational imaging models with higher order terms or/and non-differential terms [39,40,44].These models include the classical Rudin-Osher-Fatemi (ROF) model [36] and the Euler's elastica based models [2,11].The idea of ALMs is first to convert the minimization of those functionals to constrained problems, from which appropriate augmented Lagrangian functionals can be constructed.Based on the theory of optimization, the search of minimizers of the original functional amounts to finding saddle points of the resulting Lagrangian functionals.To get saddle points, one only needs to deal with several lower order subproblems, some of which can be solved using FFT while others have closed-form solutions.Therefore, with the utilization of ALMs, those issues relating to higher order terms and non-differentiable terms can be solved or avoided.
Following the idea in the work [10], we may introduce a new function u = H(φ), and as ∇ • (∇H(φ)/|∇H(φ)|) = ∇ • (∇u/|∇u|), one can rewrite the functional (5) as follows: where u is supposed to take on the value of either 0 or 1.Note also that the curvature makes pointwise sense only for C 2 functions.By following the ideas of [10,45], one could relax the value of u inside the interval [0, 1], introduce some auxiliary variables, and thus get the new constrained minimization problem: To deal with this minimization problem, just as in our previous work [45], one could construct an augmented Lagrangian functional as follows: where D = [0, 1] and R = {m ∈ L 2 (Ω) : |m| ≤ 1 a.e. in Ω}, and δ D (v) and δ R (•) are the characteristic functions on the sets D and R respectively: The r i s, i = 1, ..., 5 are positive parameters, while λ 1 , λ 2 , λ 3 , , λ 4 , λ 5 are Lagrange multipliers.As discussed in [39], a tricky technique used here is the introduction of the variable m, which helps simplify the associated subproblem for the variable p.In fact, as m ∈ R, that is, |m| ≤ 1, then |p| − p • m ≥ 0 and the equality holds if and only if m = p/|p|.One thus avoids the quadratic term Ω (|p| − p • m) 2 , which leads to a simpler subproblem for p.Moreover, the introduction of the variable q is to tackle the non-differentiable term of mean curvature.
Using the same technique as in [45], the saddle points of this Lagrangian functional can be found by solving the resulting subproblems using either Fast Fourier Transform (FFT) or their closed-form solutions.However, to do so, one needs to choose appropriate auxiliary parameters r i s, i = 1, ..., 5.In fact, the choice of these parameters for the current segmentation model is more involved than for denoising problems.This is because for denoising problems, the cleaned images should be close to the given noisy ones, which is an intrinsic constraint that is not possessed by segmentation problems.
Basically, the principle of choosing those auxiliary parameters naturally meets two fundamental requirements.First, segmentation results should be mainly determined by the grey intensity of the image under consideration, which give the regularization term a supportive role.Second, the quantity of mean curvature of segmentation contours needs to be expressed reliably, otherwise, the above functional (8) fails to faithfully reflect the essence of the ECV -L 1 segmentation model (Eq.5).To satisfy the first prerequisite, besides the model parameters a and b, those auxiliary parameters r i s need to be chosen relatively small so that the fitting terms remain dominant.On the other hand, to fulfill the second need, one has to set r i s relatively large in order to strengthen the connection among those auxiliary variables and thus express the curvature term effectively.These two seemingly opposite prerequisites considerably narrow the feasible range of the parameters r i s.Moveover, as discussed above, theoretically, the ECV -L 1 segmentation model promotes convex contours, which also requires the curvature to be expressed accurately so that this particular feature of the model can be embodied.Therefore, the choice of these auxiliary parameters becomes a crucial issue for the numerical implementation of the ECV -L 1 segmentation model.
To overcome this issue, one possible way is to utilize other methods to represent the curvature term directly, such as the well-known Γ-convergence [15].In the literature, the notion of Γ-convergence has been successfully employed to deal with the minimization of functionals involving the L 2 -Euler's elastica for image segmentation [20,25].With this technique, the original functional can be approximated by a sequence of more tractable functionals whose minimizers converge to that of the original functional.However, to the best our knowledge, there is no related work on functionals involving the L 1 -Euler's elastica.
In this paper, we still stick to augmented Lagrangian methods, but intend to develop a different augmented Lagrangian functional with fewer Lagrange multipliers for the ECV -L 1 segmentation model.Besides considerably ameliorating the issue of choosing auxiliary parameters, due to the fewer new variables involved, this novel ALM can capture the mean curvature more directly and reliably.

A novel augmented Lagrangian method for the ECV -L 1 model
In this work, we adhere to the original form of the model (Eq.5), that is, we will use a level set function instead of a binary function as in (Eq.6), to represent the segmentation contour.This is because the mean curvature of the contour can be expressed more accurately in terms of a level set function than a binary function after discretization.
Besides keeping the level set function, we plan to develop a new augmented Lagrangian functional for the ECV -L 1 segmentation model with fewer Lagrange multipliers, but with a more tight connection among those auxiliary variables that describe curvature.In [39], as discussed above, in the Euler's elastica based image denoising model, the curvature term ∇•(∇u/|∇u|) was expressed as ∇•n with n = m, m = p/|p|, and p = ∇u.The introduction of m is to relax the direct connection of n = p/|p|.This new variable considerably simplifies the treatment of the associated subproblem for p. Due to the intrinsic constraint that denoised images should be close to the noisy ones, these relaxations often work very well for denoising problems [39].However, these relaxations become problematic when dealing with curvature terms in segmentation problems.As discussed above, the main reason lies in the fact that with too many relaxations, the curvature terms often cannot be expressed accurately.
To fix this issue, we first approximate H(φ) and ∇H(φ) as in [8] as: where > 0 is a small parameter.We can then develop a new augmented Lagrangian functional with a direct connection between p and n for the original model (Eq.5) as follows: In this new functional, the variables p and n are connected through |p|n − p = 0.This functional involves only three Lagrange multipliers, instead of five as the augmented Lagrangian functional (8), which considerably reduces the effort of choosing appropriate auxiliary parameters for the numerical implementation.
With this new Lagrangian functional (9), as discussed in [39,45], one just needs to find its saddle points in order to get local minimizers of the original functional (5).To do so, we employ the standard iterative strategy: we minimize the corresponding functional for each variable by fixing the other ones, and then update the Lagrange multipliers; this process will be repeated until all the variables are convergent, indicating that some saddle point has been obtained.In what follows, we discuss the subproblems associated with all the variables.Specifically, the corresponding functionals can be expressed as follows: As discussed in [39,45], the minimizers of the functionals ε 1 (q), ε 5 (c 1 ), and ε 6 (c 2 ) can be obtained explicitly, and the ones for ε 2 (φ) and ε 4 (n) are determined by the corresponding Euler-Lagrange equations.For the functional ε 3 (p), it possesses a non-differential term involving |p|p and therefore its Euler-Lagrange equation can only be obtained for its regularized version, which introduces approximation errors and makes the Euler-Lagrange equation intractable numerically.This explains the motivation of introducing the new variable m in [39].In this work, we try to find its minimizer directly.Let's first reformulate the functional ε 3 (p) as follows: where c is independent of p.Note that there is no spatial derivative of p in ε 3 (p).Therefore, one just needs to consider the minimizer of its integrand at each point of the region Ω.Specifically, we only need to find the minimizer of the function g : R 2 → R as follows: where λ, µ ∈ R, a, ν ∈ R 2 are given with µ > 0. To this end, we propose the following theorem.
Theorem 2. Assume that µ > 2|ν|.Let θ be the angle between the vector a and the minimum vector of g(x) listed in (17), and α the angle between a and ν.Then the following arguments hold: • if λ ≥ µ|a|, then g(x) attains its minimum at x = 0. • if λ < µ|a|, the minimum point of g(x) can be determined according to the following four cases: 1. if a = ν = 0, the minimum occurs at x = 0 if λ ≥ 0 and any vector of length 2. if a = 0, ν = 0, the minimum point can be expressed as )a; 3. if a = 0, ν = 0, the minimum occurs at λ µ − 2|ν| where [ν a] denotes the 2 × 2 matrix with the vectors ν and a being the first and second column respectively.
On the other hand, if λ < µ|a|, for any unit vector b, one has either µ(b • a) − λ ≤ 0 or µ(b • a) − λ > 0. For the first case, t * ≤ 0, and h(t) takes its minimum value µ 2 |a| 2 , while for the second case, the minimum value reads µ , which is strictly less than µ 2 |a| 2 .Therefore, to find the minimum of g(x), one needs to find the maximum of the function ψ In what follows, let's consider those cases under the assumption λ < µ|a| one by one.If a = ν = 0, then g(x) = λ|x| + µ 2 |x| 2 .As µ > 0 and λ < µ|a| = 0, simple calculation shows that the minimum value occurs at all those x with length −λ/µ.
If a = 0, ν = 0, ψ(b) degenerates to be ψ(b) = .Geometrically, as shown in Fig. 1, the smaller the angle θ is, the bigger the numerator of ψ(b) will be; while for its denominator, if the angle between b and ν is larger, the denominator will be smaller.These facts show that to attain the maximum value of ψ(b), b should sit on the further side of a with respect to ν.This requires to use the sign of the determinant of the 2 × 2 matrix [ν a] to determine the orientation of b with respect to a. Specifically, if det[ν a] ≥ 0, we may rotate the unit vector a/|a| counter-clockwisely to get b; while if det[ν a] < 0, b can be obtained by rotating the vector clock-wisely.In this manner, one finds the maximum b for ψ(b), and thus gets the minimum of the goal function g(x) as µ+2ν•b b, which proves the theorem.
Remark 1.Notice that if one expresses the integrand of ε 3 (p) as the function g(p) defined as (17), those parameters read λ and ν = −r 3 n.Then it is easy to see that the assumption µ > 2|ν| is surely satisfied, since r 3 (1 + |n| 2 ) ≥ 2r 3 |n| for any n.Therefore, the above theorem can be applied for the minimization of the functional ε 3 (p).
Remark 2. Based on this theorem, the closed-form minimizer of ε 3 (p) only exist for several extreme cases.In general, one has to solve the equation (18).In this paper, we solve it using Newton's method with an initial θ satisfying [cos θ, sin θ] = a/|a|, that is, choose the unit vector along with the vector a as the initial guess of b when maximizing the function ψ(b).
In fact, if the minimum θ is close to 0, by using the approximation that cos(θ + α) ≈ cos α and sin(θ + α) ≈ sin α, one gets an approximation and then we may also use this particular θ as the initialization when solving (18)  Remark 3. In fact, the constraint |p|n−p = 0 was also employed in [17].However, in this work, we propose a novel method to deal with the associated subproblem of p directly, while in [17], a fixed-point based technique was used.Specifically, during the iterative process, this constraint is relaxed to be |p k−1 |n − p = 0 in the k th iteration, and the minimization of the subproblem of p is circumvented and an approximation of its minimizer is obtained the same technique as in [39].
Similarly as was done in [39,45], as ε 1 (q) can be rewritten as where c 1 is independent of q, its minimizer reads with q * = ∇ • n − λ 2 /r 2 .And the minimizers of ε 5 (c 1 ), and ε 6 (c 2 ) are given as follows: As for the functionals ε 2 (φ) and ε 4 (n), standard procedures lead to the following Euler-Lagrange equations: In summary, as discussed above, to find saddle points of the new augmented Lagrangian functional (9), we minimize each of the above functionals ε i s, i = 1, • • •, 6 by fixing other variables, and once all these variables are updated, we then advance the Lagrange multipliers λ 1 , λ 2 , λ 3 as follows: (i, j) is called a pixel point.All the variable functions are defined on these pixel points.We first introduce the discrete backward and forward differential operators with periodic boundary condition as follows: and the central difference operators and the gradient operators are defined accordingly as We first discuss how to solve the equation ( 27) using FFT.As this equation contains nonlinear terms of φ, we employ the frozen coefficient technique to solve it, and consider the following equation where δ φ > 0 is a small parameter, and its discretization as follows: where p = p 1 , p 2 and λ 1 = λ 11 , λ 12 .Let's denote the right side by the function g and apply the discrete Fourier transform F for both sides.Notice that where Once Fφ is calculated, φ can be obtained using the discrete inverse Fourier transform.
As for the equation (28), notice that the coefficient r 3 |p| 2 varies over Ω and also might not always be positive.We hence consider the following equation where D = max x∈Ω (r 3 |p| 2 +δ n ) with δ n being a small positive number.As n = n 1 , n 2 , p = p 1 , p 2 , the discretization of (40) leads to a system, Denote the two right sides by the functions h 1 , h 2 respectively.One applies the discrete Fourier transform to both sides and gets the following 2 × 2 system for (Fn 1 , Fn 2 ) at each pixel point (i, j): with and Notice that the determinant of the above 2 × 2 matrix equals to D 2 − 2Dr 2 (cos z 1 i + cos z 2 j − 2) > 0. One can easily solve the above system to get Fn 1 and Fn 2 and then apply the discrete inverse Fourier transform to find the new updated n = (n 1 , n 2 ).
As for the update of the variable p for each iteration, in general, one needs to solve the equation (18).In this work, we employ the Newton's method with an initial θ that was discussed in Remark 3.2, that is, θ satisfies [cos θ, sin θ] = a/|a| with a being the vector defined in the function g(x) (Eq.17) or θ is the angle defined in Eq. (22).
the model ignores objects of relatively small sizes while only keeping bigger ones.This feature is also possessed by other models for image denoising, such as the T V -L1 image denoising model [9] and the mean-curvature denoising model [43].In fact, it is very natural for the proposed model to assume this feature, since as discussed in Theorem 1, each closed segmentation contour Γ leads to at least 2π for the term Γ |κ|ds, and if the region enclosed by Γ is relatively small, it will be omitted if this causes a change of the fitting terms that is less than 2π.
Secondly, from these plots, especially the middle in the second row, one observes that the staircase-like shape is taken out from the final segmentation contour.Notice that the three shapes in the first row of the image share the same area, that is, the fitting terms of the model will treat them equally.Therefore, it is the L 1 -Euler's elastica term that separates the two convex regions from the non-convex one.By the same argument as above, whether a region is kept in the final segmentation depends on the competition between the fitting terms and the L 1 -Euler's elastica term.This phenomenon also suggests a new shape driven property of the model, that is, it prefers objects of convex shapes.To the best of our knowledge, no discussion about this property has been given for other segmentation models, such as the Mumford-Shah model [30], the Caselles-Kimmel-Sapiro model [7] or the Chan-Vese model [8].
From the depicted results, one may also note that as the parameter b increases, the contour on the "pac-man" (the right on the second row) becomes more and more convex, since this reduces the L 1 −Euler's elastica energy.This observation again illustrates the effect of the curvature term in the model.

Experiments for real images and the effect of parameters
In Fig. 3, we apply the model to a real image with a boy wearing a big black hat.This hat is partially occluded by the boy's face and neck.Image intensity based models, such as the Chan-Vese model, generally find the visible part of the hat with a big notch near the boy's neck.However, when the proposed segmentation model is applied to this image and with suitable model and Lagrange parameters, this indentation on the segmentation contour can be replaced by an almost straight segment in order to form a convex contour as shown in Fig. 3.
In Fig. 4, we consider another real image inside which there is a mushroom.The mushroom consists of two parts -a stem and a cap, which form a nonconvex shape.By applying the proposed model with a relatively large curvature parameter b, one notices that only the larger cap part is included in the final segmentation, which cannot be achieved by the conventional intensity based segmentation models.
The above two experiments on real images demonstrate that the proposed segmentation model promotes convex segmentation contours once the curvature parameter b is chosen large.Moreover, to get those convex contours, the model may automatically restore missing or occluded parts or subtract nonconvex parts.
To see how the curvature parameter b affects the result of the segmentation, in Fig. 5 and Fig. 6, we present the results for the same two experiments with different values of b while keeping all the other parameters.From these results, one sees that once b is small, just as expected, the performance of the new model degenerates to that of the original Chan-Vese's model.These results also demonstrate that the curvature is captured accurately with the proposed algorithm.To check whether the iteration used for the proposed ALM converges to a saddle point of the augmented Lagrangian functional (Eq.9), we present the plots of the relative residuals (Eq.43), relative errors of the Lagrange multipliers (Eq.44), and relative error of the iterative u k (Eq.45) in Fig. 7.All these plots indicate that a saddle point of the functional and thus a minimizer of the proposed model (Eq.5) is approached.

Reinitialization of the level set function φ
In the proposed model (Eq.5), a level set function φ is used for locating the segmentation contours.Differently from the alternative model (Eq.6) with u ∈ [0, 1], no constraint on the value of φ is imposed.Hence, when a stationary state of the segmentation contour is attained, that is, the region {φ > 0} stops evolving, the function φ may still keep changing, even though these updates become smaller and smaller.This is because of the choice of the smooth regularized version of the Heaviside function (Eq.9).In fact, in the region {φ < 0}, when φ → −∞, the term 1  2 + 1 π arctan( φ ) → 0, which decreases the first part of the fitting term.Similarly, in the region {φ > 0}, the value of φ will keep increasing during the iterative process.This phenomenon also explains the slow change of the plots of the relative error of the iterative u k in Fig. 7, even though the set {φ > 0} is invariant just after several hundred of iterations.
As discussed previously, the reason of sticking to a level set function, instead of a binary function, in the model (Eq.5) is to more accurately capture the mean curvature of segmentation contours.Indeed, since the level set function φ presents a few layers of its value around its zero level curve, the curvature of segmentation contour can be represented accurately.Therefore, the above issue that φ keeps changing won't affect the performance of the model as long as this changing is not abrupt.
However, if the curvature parameter b has to be chosen very large for some applications, the evolution of φ might become more intractable.The well-known procedure called reini- Figure 7: The plots of relative residuals (Eq.43), relative errors in Lagrange multipliers (Eq.44), and relative error in u k (Eq.45) for the two examples "hat" (Left column) and "mushroom" (Right column).These plots demonstrate the convergence of the iterative process.
tialization [34] for φ is recommended during the iterative process.This will help rationalize the values of φ around its zero level set so that the mean curvature can be effectively captured in order to realize the essence of the proposed model.

Conclusion
This work described a variational image segmentation model, where the L 1 -Euler's elastica energy was used as a regularizer on the segmentation contour.It was shown that an interesting feature of this new model is the ability to promote segmented regions with convex shapes.We also developed a novel augmented Lagrangian method for minimizing the associated functional.Compared with previous work on the L 2 -Euler's elastica energy [39,44], the proposed algorithm employed fewer Lagrange multipliers, which led to more direct connections among the variables that described curvature and thus helped to express the curvature effectively and accurately.Experiments demonstrated that the new segmentation model could be useful for segmenting objects that were known to have a convex shape.
In case occlusions prevented one from observing the entire object, or the transition with the background was not clear by the contrast, meaningful results were still obtained due to the additional convexity information.The experiments also indicated data driven and shape driven properties of the segmentation model, the latter of which has never been discussed in the literature for conventional variational segmentation models to our knowledge.

Figure 1 :
Figure 1: Two vectors ν, a and a possible unit vector b that maximizes the function ψ(b).

Figure 2 :
Figure 2: The original image (the first row) and the segmentation results (the second row)with different curvature parameter b.Two features can be observed from these results: 1) data driven property: as the parameter b increases, objects of relatively small size will be omitted in the final segmentation; 2) shape driven property: with the same parameter b, among those objects with equal areas, the one without a convex shape is taken out from the final segmentation (see the staircase-like shape).In this experiment, the other parameters are: a = 10 −4 , r 1 = 50, r 2 = 20, r 3 = 5.

Figure 3 :
Figure 3: The original image and the segmentation result.The result demonstrates that the invisible part of the hat can be restored by the proposed segmentation model with a relatively large curvature parameter b.In this experiment, the parameters are a = 10 −3 , b = 80, r 1 = 60, r 2 = 40, r 3 = 10.

Figure 4 :Figure 5 : 1 b = 5 . 0 Figure 6 :
Figure 4: The original image and the segmentation result.The result shows that only the "cap" part of the mushroom is captured while the "stem" is omitted in order to form a convex segmentation contour by the proposed segmentation model with a relatively large curvature parameter b.In this experiment, the parameters are a = 10 −4 , b = 18, r 1 = 20, r 2 = 10, r 3 = 5.