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

In order to separate completely the objects with larger occluded boundaries in an image, we devise a new variational level set model for image segmentation combing the recently proposed Chan-Vese-Euler model with elastica and landmark constraints. For computational efficiency, we deign 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 solved via simple Gauss-Seidel iterative method, or generalized soft thresholding formulas with projection methods respectively. Numerical experiments show that the proposed model not only can recover larger broken boundaries, but also can improve segmentation efficiency, decrease the dependence of segmentation on tuning parameters and initialization.


Introduction
Variational level set methods [1] have been widely applied to image segmentation in image analysis based on image features, such as edge, region, texture and motion, etc. [2][3][4]. For the images with occluded objects, prior shape model can be augmented into a variational model to inpaint the missing boundaries based on pre-defined shapes [5][6][7][8]. But in numberous cases, obtaining a prior shape is not an easy work, so appending some landmark constraints in a variational model to get the desired complete boundaries maybe a good alternative.
Motivated by image registration with landmarks [9][10][11], Pan et al. [12] proposed a Chan-Vese model [13] with landmark constraints (CVL) under the variational level set framework. This model not only can enforce the curve for segmentation to pass through some obviously feature points accurately, but also can improve computational efficiency and weaken the dependences on initialization. However, since the Chan-Vese model uses the total variation (TV) [14] of Heaviside function of level set function to approximate the length of contours, the CVL prefers to recover the smaller occluded boundaries. Fortunately, the elastica regularizers proposed in earlier 1990s to depth segmentation [15] have been successively used to image inpainting with larger broken images [16], image restoration with smooth components [17,18], image segmentation with occlusions [19][20][21][22]. In [19], Zhu  The task of two-phase segmentation of a gray value image f (x) : Ω → R is to divide Ω into two regions Ω 1 , Ω 2 , such that Ω = Ω 1 Ω 2 and Ω 1 Ω 2 = ∅. The classical Chan-Vese model [2] as a reduced piecewise constant Mumford-Shah model [3] under variational level set framework leads the original image to be denoted as f (x) = c 1 χ 1 (φ (x)) + c 2 χ 2 (φ (x)), where, c 1 and c 2 are constant image intensities in Ω 1 , Ω 2 respectively, χ 1 (φ (x)) = H (φ (x)) ∈ [0, 1] and χ 2 (φ (x)) = H (φ (x)) ∈ [0, 1] are characteristic functions of Ω 1 , Ω 2 respectively. φ (x) is a level set function defined as a signed distance function, i. e. (2. its partial derivative with φ (x) is the Dirac function which is a generalized function. Usually, H (φ), δ (φ) are replaced by their mollified versions by introducing a small positive constant parameter ε , for instance [10] H ε (φ) = 1 2 1 + 2 π arctan φ ε , (2.6) After this we can estimate c 1 , c 2 by variational method, c 1 and c 2 can be estimated as (2.8) Using the above mentioned expressions, the Chan-Vese model can be written as given input image f . This fidelity term is the same as in the original Chan-Vese model [13]. The γ is a penalty parameter for the length term of the curve. The (2.10)

2.2
The Chan-Vese model with landmark constraints [12] Let x L = {x 1 , x 2 ...x l } be the landmark points, on this basis we can also select a threshold to mark the parts of the image that we are interested in,we can use a mask function η (x) ∈ {0, 1} to indicate their positions as Since the zero level sets are used to describe the boundary curve, the landmark constraints are So, the Chan-Vese model (2.9) is transformed into the following constrained optimization problem Introducing a penalty parameter µ > 0, (2.13) is approximated as a new form This is the so-called CVL model in [12].

The Chan-Vese model with elastica
In order to recover curves which are not determined by image features, [19,28] proposed the CVE model combining Chan-Vese model and the elastica term as where (∇ · ∇φ |∇φ| ) 2 is the elastica, i. e. the square of the curvature. Due to [19] studied the corresponding convex optimization counterpart of (2.15) as But the investigations of this paper are based on (14) for curve evolution control.

The CVE model with landmark constraints and its ADMM algorithm
Similar to (2.14), based on (2.15), we propose the Chan-Vese-Elastica model as with |∇φ| = 1.In order to simplify the implementation of (3.1), we introduce auxiliary variables u, p, m, n such that p = ∇φ, Due to |m ≤ 1|, (3.3) can be replaced by a new constraint |p| − p · m ≤ 1 and |m ≤ 1| for relaxification as [29]. Regarding p = ∇φ, the constraint |∇φ| = 1 is rewritten as |p| = 1. Additionally, we introduce a new variable n = m as [29] for splitting. The constraints (3.2)-(3.4) can summarized as (3.8) In order to design the ADMM algorithm of the problem, we introduce Lagrange multipliers λ 1 , λ 2 , λ 3 , λ 4 and penalty parameters γ 1 , γ 2 , γ 3 , γ 4 to rewrite (3.1) into the following Augement Lagrangian Function as where |p| = 1, R = m ∈ L 2 (Ω) : |m| ≤ 1 a.e. in Ω and δ R (m) is the characteristic function on the convex set R, which is given by Under the framework of ADMM, the Lagrangian multiplion are updated on (3.10). . (3.10) And the original problem is transformed into the following sub-problems in each loop for k = 0, 1, 2...K. 14) Using standard variational method, we get the solutions of (3.11) and (3.12) respectively as the Euler-Lagrange equations on φ are obtained from (3.13) as (3.20) p k+1 can be calculated using generalized soft thresholding formula and projection formula as The Euler-Lagrange equations on from (3.15) The generalized soft thresholding formula of with projection formula from (3.16) The generalized soft thresholding formula of with projection formula from (3.17) 4 Implementations of the relevant sub-problems of minimization To compute (3.18)-(3.24) and (3.10) numerically, we need to discretize them firstly and then design proper algorithms. For the sake of simplicity, we discretize the image domain pixel by pixel with indices for rows and columns. Then the gradients can be represented approximately by forward, backward and central finite differences, as where, .

(4.2)
Then the discretized divergence operators and Laplacians can be stated as The other variables can be expressed in a similar way. The (3.17) and (3.18) can be calculated directly as For (3.19), we introduce intermediate variables to transform it into a consize form Based on its discretized form, the Gauss-Seidel iterative scheme can be easily designed as Fast Fourier transform(FFT) [26] can also be used here. . (4.9) Where F −1 denots as inverse Fourier transform, is the real part of a complex number, as for F i,j is given for (4.6). For (3.20), its discretized solution is For (3.21), it has the following iterative formula was same as φ: to transform it into a simple form Based on its discretized form, the Gauss-Seidel iterative scheme can be easily designed as , n k+1,0 2i,j = n k 2i,j U n k+1,l = n k+1,l i+1,j + n k+1,l i−1,j + n k+1,l i,j+1 + n k+1,l i,j−1 .

(4.13)
Here n can also be solved with FFT, the shape is very similar to (4.9), so I won't go into details here. For (3.22), the discretized analytical formula with projection formula is (4.14) For (3.23) is a very simple analytical solution where Tol means a given number,and it's value always set to 0.01, as for T k+1 s , Φ k+1 , Σ k+1 , there are the following methods to calculate The algorithm is summarized as

Numerical Experiments
In this section, we will enumerate the specific experimental results of the CVEL model. The design of these experiments is considered in four aspects: function, convenience, efficiency and practical application. Below we have carried out experiments on these four aspects.
Firstly, we mainly test the performance of the model by comparing different models. The models we selected as the comparison are CV model, CVL model, and CVE model. The main reason for this design is to compare the difference between the higher-order model and the lower-order model and the impact of the landmark on the result.

Comparisons with previous models
We first consider the broken letter UCLA. The result is shown in Fig.1. By comparing the CV model with the CVE model, we can find that the high-order model itself has a certain repair ability. The comparison between CVL model and CV model and the comparison between CVE model and CVEL model can be found that Landmark method has the effect of increasing repair capacity and increasing model stability.
Then we test the similarities and differences between the CVL model and the CVEL model through two experiments, just as shown in Fig.2 and Fig.3. Through those experiment, we found that both models can repair external damage, but the CVL model has a much higher dependence on the landmark point than the CVEL model. The CVEL model only needs a small number of landmark points to restore a relatively good circle, and the CVL model needs more dense points to be repaired. When the landmark point is fixed, the repair power of the CVEL model is stronger than that of the CVL model. As for the dependence on the landmark point, we will discuss it in the next study of parameter dependence.

The dependences on tuning parameters
This section mainly considers the influence of human factors on experimental results. The first factor is the change of parameters, and the other is the number of landmark points.
For the study of parameters, we should konw how the landmark dependence on the parameters affects. We compare the CVE model with the CVEL model by changing only one parameter with the same amplitude. The experiment that reduces γ 3 by 1 is as follows. As shown in Fig.3, the results of the CVE model change when the parameter is changed, and the CVEL model results unchanged. Then we also carried out other parameters research, found that the CVE model has lower requirements on the rule item constraint parameter α, but the γ requirements are higher, each time changing the value of γ will basically affect the experimental results, and the CVEL model does not have this limitation.  Those two ewperiment shows that in the case of insufficient landmark points, the demand for α in the CVL model becomes very demanding, and the CVEL model has a much looser demand for parameters. But how does the landmark point affect our model? We used a triangle with a broken corner as the experimental object, let the landmark point gradually observe the experimental results from scratch, and the result is shown in Fig.6.
From the experimental results, we can easily find that with the increase of landmark points, the experimental results are getting better and better, and there are many test results that are not repeated. From this large number of experiments, we have summarized the approximate location of the landmark points. The law is more demanding when it is close to the edge and the fixed point, while the demand in the middle is lower. Figure 6: Different experimental results caused by different landmark points.

Efficiency
We believe that by artificially adding landmark points, the convergence speed of the model can be accelerated, so the experiment in this section is used to test how the actual effect is. The first thing we consider is whether the landmark point can speed up the speed of the model. Since the point selection is human, we can mark the parts we are interested in through various techniques, such as the next experiment. Not for a single point, we mark the entire palm by binarization to check whether it will speed up the results. The result is shown in Fig.7. From the experimental results, we can see that the marking of the region can speed up the results, so we will analyze the convergence of the function. As we all know Iteration of ALM often leads to augmentation of the saddle point of the Lagrangian function, so we should to check the evolution of relative error in Lagrange multipliers , relative error of H (φ), and the evolution of energy E(φ) , In Fig.8, we present the plots of these quantities versus iterations for the above two examples.
The results shows that the function converges very fast. As for why the energy in the UCLA picture will rise, according to our research, it may be because the first iteration is fast to convergence, at this time c 1 , c 2 plays a major role, then the two items are divided into external landmark points. The effect is that we force the level set function to pass the lanmark point, and the energy may pick up.

Applications in real life
Here we use a simple example to introduce the apply of this model. In the image segmentation, there is a need for the segmentation of the image to be incomplete. There are also reasons for the image itself. For example, the CT image, because the information of the image itself does not change much, the model often leads to incomplete information when segmentation is performed. But our model can better segment the target image while ignoring the noise. We have done experiments on classic brain CT images. The results are shown in Fig.9. The results show that by adding landmark points, we can segment the results more accurately and faster.
Then the main part of our model is flexible, and the direction of segmentation can be determined by skillful artificial choices, more interesting that this model even can delete the c 1 and c 2 , such as Fig10 and Fig.11 is the result of this interesting experiments, Fig.10 shows that artificially calibrated landmark points can lead to two different results. Fig.11 indicates that there is less image information, as long as a large number of landmark points are used to mark the damaged area, you can get good results.