The Chan-Vese Model with Elastica and Landmark Constraints for Image Segmentation

TL;DR

This paper introduces a new variational level set model for image segmentation that combines the Chan-Vese model with elastica and landmark constraints, improving boundary recovery and segmentation efficiency.

Contribution

It proposes a novel model integrating elastica and landmark constraints with the Chan-Vese framework, along with an efficient ALM/ADMM optimization algorithm.

Findings

Enhanced boundary recovery in occluded regions

Improved segmentation efficiency and robustness

Reduced dependence on parameter tuning and initialization

Abstract

In order to completely separate objects with large sections of occluded boundaries in an image, we devise a new variational level set model for image segmentation combining the Chan-Vese model with elastica and landmark constraints. For computational efficiency, we design its Augmented Lagrangian Method (ALM) or Alternating Direction Method of Multiplier (ADMM) method by introducing some auxiliary variables, Lagrange multipliers, and penalty parameters. In each loop of alternating iterative optimization, the sub-problems of minimization can be easily solved via the Gauss-Seidel iterative method and generalized soft thresholding formulas with projection, respectively. Numerical experiments show that the proposed model can not only recover larger broken boundaries but can also improve segmentation efficiency, as well as decrease the dependence of segmentation on parameter tuning and initialization.

Tables

Algorithm 1: ADMM
1: Initialization: Set ${\alpha }_1,{\alpha }_2,\mu ,a,b$.
2: while any stopping criterion is not satisfied do
Calculate $c_1^{k+1},c_2^{k+1}$ from (3.18) and (3.19)
Calculate ${\varphi }^{k+1}$ from (3.20)
Calculate ${𝒑}^{k+1}$ from (3.21)
Calculate ${𝒏}^{k+1}$ from (3.22)
Calculate ${𝒎}^{k+1}$ from (3.23)
Calculate $q^{k+1}$ from (3.24)
Calculate ${\lambda }_1^{k+1},{𝝀}_2^{k+1},{𝝀}_3^{k+1},{\lambda }_4^{k+1}$ from $\left(3.10\right)$
3: end while

Equations

$$\phi\left(x\right)=\left\{\begin{array}[]{rcl}&d\left(x,\Gamma\right),&\text{if}\;x\in\Omega_{1}\\ &0,&\text{if}\;x\in\Gamma\\ &-d\left(x,\Gamma\right),&\text{if}\;x\in\Omega_{2}\\ \end{array}\right.,$$

$$|\nabla \varphi \left(x\right)|=1.$$

$$H\left(\phi\left(x\right)\right)=\left\{\begin{array}[]{rcl}&1,&\text{if}\;\phi\left(x\right)\geq 0\\ &0,&\text{otherwise}\end{array}\right.,$$

$$\delta \left(\varphi \right)=\frac{\partial H\left(\varphi \right)}{\partial \varphi },$$

$$H_{\epsilon }\left(\varphi \right)=\frac{1}{2}\left(1+\frac{2}{\pi }arctan\left(\frac{\varphi }{\epsilon }\right)\right),$$

$${\delta }_{\epsilon }\left(\varphi \right)=\frac{\partial H_{\epsilon }\left(\varphi \right)}{\partial \varphi }=\frac{1}{\pi }\frac{{\epsilon }^2}{{\epsilon }^2+{\varphi }^2}.$$

$$\text{min}E(c_1,c_2,\varphi )={\int }_{\Omega }{\left(f-c_1\right)}^2{H}_{\epsilon }\left(\varphi \right)dx+{\int }_{\Omega }{\left(f-c_2\right)}^2\left(1-H_{\epsilon }\left(\varphi \right)\right)dx+\gamma {\int }_{\Omega }|\nabla H_{\epsilon }\left(\varphi \right)|dx,s.t.|\nabla \varphi |=1.$$

$$c_1=\frac{{\int }_{\Omega }f\left(x\right)H_{\epsilon }\left(\varphi \left(x\right)\right)dx}{{\int }_{\Omega }{H}_{\epsilon }\left(\varphi \left(x\right)\right)dx},$$

$$c_2=\frac{{\int }_{\Omega }f\left(x\right)\left(1-H_{\epsilon }\left(\varphi \left(x\right)\right)\right)dx}{{\int }_{\Omega }\left(1-H_{\epsilon }\left(\varphi \left(x\right)\right)\right)dx}.$$

$$\text{min}E(c_1,c_2,\varphi )={\int }_{\Omega }Q\left(c_1,c_2\right)H_{\epsilon }\left(\varphi \right)dx+\gamma {\int }_{\Omega }|\nabla H_{\epsilon }\left(\varphi \right)|dx,s.t.|\nabla \varphi |=1.$$

$$\left\{\begin{split}&\frac{\partial\phi\left(x,\;t\right)}{\partial t}=\left(\nabla\cdot\left(\frac{\nabla\phi\left(x,\;t\right)}{|\nabla\phi\left(x,\;t\right)|}\right)-Q\left(c_{1},c_{2}\right)\right)\delta_{\varepsilon}\left(\phi\left(x,\;t\right)\right)&t>0,\;x\in\Omega\\ &\frac{\partial\phi\left(x,\;t\right)}{\partial N}=0&t>0,\;x\in\partial\Omega\\ &\phi\left(x,\;0\right)=\phi^{0}\left(x\right)&t=0,\;x\in\Omega\end{split}\right.,$$

$$\eta\left(x\right)=\left\{\begin{array}[]{rcl}&1,&\text{if}\;x\in x_{L}\\ &0,&\text{otherwise}\end{array}\right..$$

$$\varphi \left(x\right)=0,\text{if}\eta \left(x\right)=1.$$

$$\text{min}E(c_1,c_2,\varphi )={\int }_{\Omega }Q\left(c_1,c_2\right)H_{\epsilon }(\varphi )dx+\gamma {\int }_{\Omega }|\nabla H_{\epsilon }\left(\varphi \right)|dx+\frac{\mu }{2}{\int }_{\Omega }\eta {\varphi }^2,s.t.|\nabla \varphi |=1,$$

$$\text{min}E(c_1,c_2,\varphi )={\int }_{\Omega }Q\left(c_1,c_2\right)H_{\epsilon }(\varphi )dx+\gamma {\int }_{\Omega }\left(a+b{\left(\nabla \cdot \left(\frac{\nabla \varphi }{|\nabla \varphi |}\right)\right)}^2\right)|\nabla H_{\epsilon }\left(\varphi \right)|dx,s.t.|\nabla \varphi |=1,$$

$$\frac{\nabla H_{\epsilon }\left(\varphi \right)}{|\nabla H_{\epsilon }\left(\varphi \right)|}=\frac{\nabla \varphi \delta \left(\varphi \right)}{|\nabla \varphi |\delta \left(\varphi \right)}=\frac{\nabla \varphi }{|\nabla \varphi |}.$$

$$\text{min}E(c_1,c_2,\varphi )={\int }_{\Omega }Q\left(c_1,c_2\right)\varphi dx+\gamma {\int }_{\Omega }\left(a+b{\left(\nabla \cdot \left(\frac{\nabla \varphi }{|\nabla \varphi |}\right)\right)}^2\right)|\nabla \varphi |dx,s.t.|\nabla \varphi |=1.$$

$$\begin{array}{rl}\text{min}E(c_1,c_2,\varphi )& ={\int }_{\Omega }Q\left(c_1,c_2\right)H_{\epsilon }(\varphi )dx+\frac{\mu }{2}{\int }_{\Omega }\eta {\varphi }^{2}dx\\ & +\gamma {\int }_{\Omega }\left(a+b{\left(\nabla \cdot \left(\frac{\nabla \varphi }{|\nabla \varphi |}\right)\right)}^2\right)|\nabla H_{\epsilon }\left(\varphi \right)|dx,\\ & s.t.|\nabla \varphi |=1.\end{array}$$

$$\mathbf{p}=\nabla \varphi ,$$

$$\mathbf{m}=\frac{\mathbf{p}}{|\mathbf{p}|},$$

$$q=\nabla \cdot \mathbf{n}.$$

$$\mathbf{p}=\nabla \varphi ,$$

$$|\mathbf{p}|-\mathbf{p}\mathbf{m}=0,|\mathbf{p}|=1,$$

$$\mathbf{n}=\mathbf{m},|\mathbf{m}|\le 1,$$

$$q=\nabla \cdot \mathbf{n}.$$

$$\begin{array}{rl}E\left(c_1,c_2,u,\varphi ,\mathbf{v},\mathbf{p},\mathbf{n},\mathbf{m}\right)& ={\int }_{\Omega }Q\left(c_1,c_2\right)H\left(\varphi \right)dx\\ & +\gamma {\int }_{\Omega }\left(a+b{q}^2\right)|\mathbf{p}|{\delta }_{\epsilon }\left(\varphi \right)dx+\frac{\mu }{2}{\int }_{\Omega }\eta {\varphi }^{2}dx\\ & +{\int }_{\Omega }{\lambda }_1\left(|\mathbf{p}|-\mathbf{p}\cdot \mathbf{m}\right)dx+{\gamma }_1{\int }_{\Omega }\left(|\mathbf{p}|-\mathbf{p}\cdot \mathbf{m}\right)dx\\ & +{\int }_{\Omega }{\mathbf{\lambda }}_2\left(\mathbf{p}-\nabla \varphi \right)dx+\frac{{\gamma }_2}{2}{\int }_{\Omega }|\mathbf{p}-\nabla \varphi {|}^{2}dx\\ & +{\int }_{\Omega }{\mathbf{\lambda }}_3\left(\mathbf{n}-\mathbf{m}\right)dx+\frac{{\gamma }_3}{2}{\int }_{\Omega }{\left(\mathbf{n}-\mathbf{m}\right)}^{2}dx+{\delta }_{\mathcal{R}}(\mathbf{m})\\ & +{\int }_{\Omega }{\lambda }_4\left(q-\nabla \cdot \mathbf{n}\right)dx+{\gamma }_4{\int }_{\Omega }{\left(q-\nabla \cdot \mathbf{n}\right)}^{2}dx\end{array},$$

$$\delta_{\cal{R}}(\bm{m})=\left\{\begin{array}[]{ll}0&\text{if}\;\bm{m}\in{\cal{R}}\\ +\infty&\text{otherwise}\end{array}\right..$$

$$\{\begin{array}{rl}& {\lambda }_1^{k+1}={\lambda }_1^k+{\gamma }_1\left(|{\mathbf{p}}^{k+1}|-{\mathbf{p}}^{k+1}\cdot {\mathbf{m}}^{k+1}\right)\\ & {{\mathbf{\lambda }}_2}^{k+1}={{\mathbf{\lambda }}_2}^k+{\gamma }_2\left({\mathbf{p}}^{k+1}-\nabla {\varphi }^{k+1}\right)\\ & {{\mathbf{\lambda }}_3}^{k+1}={{\mathbf{\lambda }}_3}^k+{\gamma }_3\left(|{\mathbf{n}}^{k+1}|-{\mathbf{m}}^{k+1}\right)\\ & {\lambda }_4^{k+1}={\lambda }_4^k+{\gamma }_4\left(q^{k+1}-\nabla \cdot {\mathbf{n}}^{k+1}\right)\end{array},$$

$$c_1^{k+1}=\text{arg}\underset{c_1}{\min }E\left(c_1,c_2^k,{\varphi }^k,{\mathbf{p}}^k,{\mathbf{n}}^k,{\mathbf{m}}^k,q^k\right),$$

$$c_2^{k+1}=\text{arg}\underset{c_2}{\min }E\left(c_1^{k+1},c_2,{\varphi }^k,{\mathbf{p}}^k,{\mathbf{n}}^k,{\mathbf{m}}^k,q^k\right),$$

Full text

††journal: Journal of LaTeX Templates

\addtokomafont

labelinglabel

url]2017021234@qdu.edu.cn myfootnotemyfootnotefootnotetext: Since 1880.

url]huizhu.pan@postgrad.curtin.edu.au

url]w.liu@curtin.edu.au

url]zisen_xu@126.com

1 Introduction

Variational level set methods[1] have been widely applied to image segmentation problems based on image features such as edge, region, texture and motion, etc.[2, 3, 4, 5]. For images containing occluded objects, variational models using shape priors can inpaint missing boundaries based on pre-defined shapes[6, 7, 8, 9]. However, obtaining shape priors is often not easy. In certain scenarios, using landmarks is a good alternative to using shape priors in the segmentation of occluded objects. Motivated by image registration with landmarks[10, 11, 12], Pan et al.[13] proposed a Chan-Vese model[14] with landmark constraints (CVL) under the variational level set framework. The model not only enforces the segmentation contour to pass through some pre-selected feature points but also improves computational efficiency and weakens the dependence of the segmentation result on initialization. However, since the Chan-Vese model uses the total variation (TV) [15] of the Heaviside function of the level set function to approximate the length of contours, sometimes the details are not handled well enough. So the CVL model performs better for recovering smaller boundaries. On the other hand, the elastica regularizers proposed in the early 1990s in depth segmentation[16] have been successively applied to image inpainting with larger broken images[17], image restoration with smooth components[18, 19], and image segmentation with larger damaged areas or occlusions[20, 21, 22, 23]. In [20], Zhu et al. propose a modified Chan-Vese model with elastica (CVE), combining the classic Chan-Vese model (CV) and the elastica regularizer to inpaint, or interpolate, segmentation curves. Due to the ill-posed nature of this model, the segmentation results rely heavily on the involved penalty parameters, which makes it hard to for curves to pass through desired points in the occluded regions. In this paper, we propose a CVE model with landmark constraints (CVEL) that combines the CVL and CVE to more accurately and robustly complete missing curves. Different from the CVE proposed in [20] which uses piecewise constant level set functions or binary label functions, we use the Lipschitz smooth level set function defined as a signed distance function to describe curve evolution. In order to solve the proposed model with the signed distance property and landmark constraints, we devise its Augmented Lagrangian Method (ALM), i.e., Alternating Direction Method of Multipliers (ADMM)[24, 25, 26, 27] solution by dividing the original problem into several simple sub-problems and optimizing them alternatively. The sub-problems can be solved respectively by the Gauss-Seidel iterative method and generalized soft thresholding formula with projection [28].

The paper is organized as follows. In section 2, we present the classical Chan-Vese model, the Chan-Vese model with elastica, and the Chan-Vese model with landmark constraints for comparisons consequently. In section 3, we give the CVE model with landmark constraints under the variational level set framework and design its ADMM method. Section 4 covers the solutions to all sub-problems derived in section 3. Numerical examples are presented in section 5 to show the performance of the proposed model and algorithm. Finally, concluding remarks are drawn in section 6.

2 The Previous Works

2.1 The Chan-Vese model for image segmentation

The task of two-phase segmentation of a gray value image $f\left(x\right):\Omega\rightarrow R$ is to divide $\Omega$ into two regions $\Omega_{1},\Omega_{2}$, such that $\Omega=\Omega_{1}\bigcap\Omega_{2}$ and $\Omega_{1}\bigcap\Omega_{2}\neq\varnothing$. The classical Chan-Vese model[14] is a reduced piecewise constant Mumford-Shah model[3] under the variational level set framework. The original image is denoted as $f\left(x\right)=c_{1}\chi_{1}\left(\phi\left(x\right)\right)+c_{2}\chi_{2}\left(\phi\left(x\right)\right)$, where $c_{1}$ and $c_{2}$ are the average image intensities in $\Omega_{1},\Omega_{2}$, and $\chi_{1}\left(\phi\left(x\right)\right)=H\left(\phi\left(x\right)\right)\in\left[0,1\right]$ and $\chi_{2}\left(\phi\left(x\right)\right)=1-H\left(\phi\left(x\right)\right)\in\left[0,1\right]$ are characteristic functions of $\Omega_{1},\Omega_{2}$ respectively. $\phi\left(x\right)$ is a level set function defined as a signed distance function form point $x$ to curve $\Gamma$, i. e.

[TABLE]

with the property

[TABLE]

(2.2) is the Eikonal equation, i. e., a kind of Hamilton-Jacobi equation. $H\left(\phi\left(x\right)\right)$ is the Heaviside function of $\phi\left(x\right)$, stated as

[TABLE]

its partial derivative with $\phi\left(x\right)$ is the Dirac function

[TABLE]

which is a generalized function. Usually, $H\left(\phi\right)$, $\delta\left(\phi\right)$ are replaced by their mollified versions by introducing a small positive constant parameter $\varepsilon$ , for instance [14]

[TABLE]

[TABLE]

The well-known Chan-Vese model[14] for two-phase image segmentation is an energy minimization problem on $c_{1}$, $c_{2}$ and $\phi$, such that

[TABLE]

Where $\gamma$ is a penalty parameter for the length term of the curve.

$c_{1}$ and $c_{2}$ are estimated as

[TABLE]

[TABLE]

By introducing $Q(c_{1},\;c_{2})=\alpha_{1}(c_{1}-f)^{2}-\alpha_{2}(c_{2}-f)^{2}$, (2.7) can be rewritten as

[TABLE]

Therefore, the evolution equation of $\phi\left(x\right)$ can be derived via variational methods and gradient descent as

[TABLE]

2.2 The Chan-Vese model with landmark constraints[13]

Let $x_{L}=\{x_{1},x_{2}...x_{l}\}$ be the given landmark points, represented through a mask function

[TABLE]

Since the zero level set describes the boundary curve and the landmarks are positioned on the boundary, the landmark constraint is

[TABLE]

Thus, the Chan-Vese model (2.9) can be transformed into the following constrained optimization problem

[TABLE]

where $\mu>0$ is the penalty parameter.

2.3 The Chan-Vese model with elastica

In order to recover curves which are not determined by image features, for instance the boundary of an occluded object, [20, 29] proposed the CVE model by combining Chan-Vese model and the elastica term

[TABLE]

where $(\nabla\cdot\left(\frac{\nabla\phi}{|\nabla\phi|}\right))^{2}$ is the elastica, i. e. the square of the curvature. The contours obtained by this method tend to be curved rather than straight. Later, Zhu et al. [20] considered the relation

[TABLE]

and studied the following convex optimization problem instead of (2.15)

[TABLE]

However, the investigations of this paper are based solely on (2.15). Because the reason for the simplification is to facilitate the calculation, and will not affect the experimental results.

3 The CVE model with landmark constraints and its ADMM algorithm

Combining (2.14) and (2.15), we propose the Chan-Vese model with elastica and landmark as

[TABLE]

By adding landmark points, you can force the contour to pass through some feature points to get good results. This is the design idea of The CVE model with landmark (CVEL). In order to simplify the implementation of (3.1), we introduce auxiliary variables $\bm{p},\bm{m},\bm{n},q$, The main reason for adding so many intermediate quantities is to avoid the curvature term appearing in the calculation and to simplify the calculation. Those are such that

[TABLE]

[TABLE]

[TABLE]

Considering that $|\bm{m}|\leq 1$, (3.3) can be substituted by a more relaxed set of constraints,  $|\bm{p}|-\bm{p}\cdot\bm{m}\leq 1$ and $|\bm{m}|\leq 1$ [3]. Since $\bm{p}=\nabla\phi$, the constraint $|\nabla\phi|=1$ can be rewritten as $|\bm{p}|=1$. Additionally, we introduce a new variable $\bm{n}=\bm{m}$ [3] for splitting. Thus, the constraints (3.2)-(3.4) can be summarized as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Next, to design the ADMM algorithm for the problem, we introduce the Lagrange multipliers $\lambda_{1},\bm{\lambda_{2}},\bm{\lambda_{3}},\lambda_{4}$ and penalty parameters $\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4}$ and rewrite the energy function in (3.1) as the following Augmented Lagrangian Function

[TABLE]

where $\left|\bm{p}\right|=1,{\cal{R}}=\left\{{m\in{L^{2}}\left(\Omega\right):|m|\leq 1\;\text{a.e.}\;\text{in}\;\Omega}\right\}$ and $\delta_{\cal{R}}(\bm{m})$ is the characteristic function on the convex set ${\cal{R}}$, given by

[TABLE]

Under the framework of ADMM, the Lagrangian multipliers are updated for iteration $k=0,1,2...K$ as

[TABLE]

and the original minimization problem is split into the following sub-problems

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The solutions to the sub-problems are presented below.

Using standard variational methods, we solve (3.11) and (3.12) respectively and get

[TABLE]

[TABLE]

For the sub-problem (3.13), the Euler-Lagrange equations on $\phi$ are

[TABLE]

where, $F^{k+1}={{Q^{k+1}}{\delta_{\varepsilon}}\left(\phi^{k}\right)+\left({a+b{q^{2}}}\right)\left|{\bm{p}^{k}}\right|\nabla{\delta_{\varepsilon}}\left(\phi^{k}\right)+\nabla\cdot\bm{\lambda}_{2}^{k}+{\gamma_{2}}\nabla\cdot{{\bm{p}}^{k}}}$.

To solve (3.14), we can derive $\bm{p}^{k+1}$ via a generalized soft thresholding formula and projection formula as below

[TABLE]

For (3.15), the Euler-Lagrange equations on $\bm{n}$ is

[TABLE]

$\bm{m}$ in (3.16) can be obtained as an exact solution. Considering the constraint in (3.7), the projection formula should be augmented as follows

[TABLE]

Lastly, $q$ in (3.17) also has an analytical solution

[TABLE]

In this part we introduce the solution of each variable in the CVEL model and the iterative process of the algorithm. In the next section we will combine the experiments to introduce the role of the model in detail.

4 Implementations of the relevant sub-problems of minimization

To compute (3.18)-(3.24) and (3.10) numerically, we need to design discrete algorithms for the sub-problems. 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

[TABLE]

where,

[TABLE]

The discretized Laplacian of $\phi$ can be stated as

[TABLE]

The other variables can be expressed in similar ways.

(3.18) and (3.19) can be calculated directly as

[TABLE]

[TABLE]

where M and N are the numbers of rows and columns of the image $f$.

Next, to discretize the formula of $\phi$ obtained in (3.20), we introduce the following intermediate variables

[TABLE]

and write the original Eular-Lagrange equations in the more concise form below

[TABLE]

Based on (4.3) and (4.7), we can easily design the Gauss-Seidel iterative scheme of $\phi$ as

[TABLE]

Alternatively, $\phi$ can be solved by Fast Fourier transform (FFT)[27].

The discretized solution of $\bm{p}$ as obtained from (3.21) is

[TABLE]

Since the form of $\bm{n}$ in (3.22) is similar to that of $\phi$, the solution of $\bm{n}$ can also be written similarly. Again, to simplify the equation, we introduce

[TABLE]

and (3.22) becomes

[TABLE]

Introducing the discretized form of $\bm{n}$, its Gauss-Seidel iterative scheme can be easily designed as

[TABLE]

Here, $\bm{n}$ can be solved with FFT as well.

For $\bm{m}$ in (3.23), its discretized analytical solution with projection formula is

[TABLE]

The $q$ obtained in (3.24) can also be drawn into a simple analytical solution

[TABLE]

After calculating (3.18) - (3.24), the Lagrange multipliers are updated as (3.10).

In each iteration, the following error tolerances should be checked to determine convergence, i. e.,

[TABLE]

where $\text{Tol}=0.01$. $T_{s}^{k+1},\;\Phi^{k+1},\;\Sigma^{k+1}$ are defined as

[TABLE]

[TABLE]

The complete algorithm is summarized in Algorithm 1.

5 Numerical Experiments

In this section, we devise four sets of numerical experiments with distinct goals in mind. The first set compares the performance of the CV, the CVL, the CVE, and the CVEL model in contour inpainting or interpolation, in cases where smaller and larger regions are missing from the original images. The second set of experiments examines the dependence of the segmentation result on the number of landmark points. The third one demonstrates how landmark points improve the segmentation efficiency of the CVEL model, and the fourth one is an application in medical image segmentation.

5.1 Comparisons with previous models

Since the CVE, CVL, and CVEL model are different extensions of the CV model for the purpose of missing contour recovery, we design some experiments to compare their performance. First, we segment an image of the letters ’UCLA’ with small damaged regions [20]. We give the original broken image in Fig.1a, and the initialization of the zero level set for all of the models in Fig.1b. Segmentation results obtained via the CV, CVL, CVE, and CVEL are shown in Fig.1c, Fig.1d, Fig.1e, and Fig.1f respectively. The parameters used in the CVEL model are $\gamma_{1}=1,\gamma_{2}=3,\gamma_{3}=5,\gamma_{4}=10,\alpha_{1}=0.5,\alpha_{2}=0.5$. One landmark was placed in the middle of each piece of missing contour. Results show that the CVL, CVE, and CVEL can all recover small sections of the missing contours, though the CVEL produces smoother curves due to the elastica regularizer.

Next, to compare the performance of the CVL and CVEL in recovering larger missing contours, we conducted the second experiment shown in Fig.2. Fig.2a shows a triangle with a missing corner, Fig.2b shows the initial contour and landmark points, Fig.2c and Fig.2d are segmentation results via the CVL and CVEL model respectively, and the parameters are $\gamma_{1}=1,\gamma_{2}=3,\gamma_{3}=5,\gamma_{4}=10,\alpha_{1}=1.1,\alpha_{2}=0.9$ for both models. Although they obtained similar results, the CVL needed 26 landmark points whereas the CVEL needed only 20.

Furthermore, we set up an experiment to compare the performance of the CVL and CVEL, as shown in Fig.3. Fig.3a shows the original broken image, Fig.3b presents the initial zero level set and landmark points, and Fig.3c and Fig.3d give the segmented results via CVL and CVEL model respectively. The parameters of CVEL model are $\gamma_{1}=7,\gamma_{2}=20,\gamma_{3}=5,\gamma_{4}=2,\alpha_{1}=1.1,\alpha_{2}=0.9$. The conclusion is the same as [20], i. e., it is hard to inpaint the external missing curve via the CVL, but the CVEL works well. This is mainly because the curvature term has good curve repair properties.

5.2 The dependence on tuning parameters and landmark points

In this part of experiments, we study the effect of the number of landmark points and their positions on the segmentation result. Setting different amounts of landmark points lead to different results. The more landmarks we set within a certain limit, the more accurate the result tends to be. However, increasing the number of landmarks beyond this limit will not increase segmentation accuracy, as shown in Fig.4.a to Fig.4d. In Fig.4e to Fig.4h, the number of landmark points is 2, 10, 18, 24, respectively. As we can see, the corresponding results became increasingly better with the more and more additional landmarks. However, setting over 24 landmarks did not improve the result further.

The placement of landmark points is also essential, especially when the total number of landmarks is small. Using Fig. 4d as an example of a well-segmented image, we proceed to take away landmarks from different locations. In Fig.5a, Fig.5b and Fig.5c, we remove two landmarks from the bottom, top, and middle of the broken tip of the triangle, respectively. As a result, the recovered contour in Fig.5d has distortions around the base section, the sharpness of the tip is not well maintained in Fig.5e, and the result in Fig.5f does not change significantly. Therefore, we observe that it is more effective to place landmarks at the vertices or corners of an object. The better the landmark positions, the fewer landmarks we need.

5.3 Efficiency

In this section, we examine how the CVEL performs in terms of efficiency, specifically the convergence time. We first consider whether the CVEL speeds up segmentation by constructing the experiment in Fig 6, where we mark the entire contour of the palm in the CVEL and compare performance with the CVE. Fig.6a is the original image, and Fig.6b is the initial zero level set for both models. We obtain the segmentation results shown in Fig.6c and Fig.6d via CVE by five steps and CVEL by two steps, respectively. This shows that using landmarks in the CVEL model can increase efficiency.

We then check the convergence time of the Lagrangian multiplier, the level set function, and the total energy (4.16-4.17) using the CVEL model, where Fig.7a, Fig.7b, and Fig.7c map their respective values in the experiment in Fig.6, and Fig.7d, Fig.7e, and Fig.7f map the same values in the experiment in Fig.1. On the one hand, we see that convergence is reached quickly in both experiments. On the other hand, we observe that the total energy increased towards the end in Fig.7f. This phenomenon is due to the model setting the landmark points onto an illusory contour. In the initial stage, the CV and elastica terms play a major role in moving the curve towards the broken boundary. As the curve approaches the landmarks, it moves away from the natural boundary as indicated by the separation of regions according to image information. This process inevitably raises the total energy.

5.4 Applications in real life

In this sub-section, we present some applications of the CVEL in segmenting CT images which often encounter difficulties due to the presence of fine details.

Fig.8a gives an original CT image of brain, Fig8.b presents the initial level set function for segmentation, Fig8.c shows the result via CV model which fails to separate the brain tissue completely. On the contrary, the CVEL model with proper landmark points and using the same level set initialization can produce good result as presented in Fig8.e. Since the missing parts in Fig8.c are areas with low image intensities, We attempt to improve quality through adding landmarks in these areas. The results in Fig8.d-Fig.8h obtained via marking points over intensity of 20,40,60,80,100, but Fig8.d presents over segmentation result due to redundant information, and Fig8.f-Fig.8h present deficient segmentation results due to less mark points.

Fig.9 gives another example of CT image segmentation with noises. The main difficulty in this experiment is to segment the adjacent areas in the image. For the original image Fig.9a, we initialize the level set function as Fig.9b. The Fig.9c and Fig.9d show the results obtained via CV model and CVEL model respectively. While using CVEL, we add landmark points over 100 of intensity for a satisfactory result.

6 Concluding remarks

In this paper, we present a Chan-Vese model with elastica and landmarks (CVEL) under the variational level set framework. The new model combines the classical Chan-Vese model (CV), the Chan-Vese model with landmarks (CVL), and the Chan-Vese model with elastica (CVE). We then design an ADMM algorithm for the solution. A variety of numerical experiments show that the CVEL performs better than the CVE in segmentation accuracy, and can recover larger broken boundaries than the CVL.For future work, we can integrate automatic landmarks detection methods as well as other techniques involving priors and ultimately aim to achieve automatic segmentation.

7 Acknowledgement

The authors thank the editor and anonymous reviewers for their helpful comments and valuable suggestions.