

ORIGINAL ARTICLE 

Year : 2018  Volume
: 4
 Issue : 1  Page : 515 

A variational level set method image segmentation model with application to intensity inhomogene magnetic resonance imaging
Chun Li^{1}, Jinhe Su^{1}, Longlong Yu^{1}, Le Wang^{1}, Luo Ze^{2}
^{1} eScience Technology and Application Laboratory, Computer Network Information Centre, Chinese Academy of Sciences; University of Chinese Academy of Sciences, Chinese Academy of Sciences, Beijing, China ^{2} eScience Technology and Application Laboratory, Computer Network Information Centre, Chinese Academy of Sciences, Beijing, China
Date of Web Publication  18May2018 
Correspondence Address: Chun Li eScience Technology and Application Laboratory, Computer Network Information Centre, Chinese Academy of Sciences; University of Chinese Academy of Sciences, Chinese Academy of Sciences, Beijing China
Source of Support: None, Conflict of Interest: None  Check 
DOI: 10.4103/digm.digm_44_17
Background and Objectives: In this article, we propose an image segmentation model based on ChanVese (CV) for image segmentation. By taking into account the local features of the image, the new proposed model can successfully segment images with intensity nonuniformity. Materials and Methods: We quantitatively compare our method with other two stateoftheart algorithms, namely, CV model and local binary fitting (LBF) model in segmenting synthetic MR images with the ground truth from BrainWeb; the data can be available at: https://www.mni/mcgill.ca/brainweb/. For segmenting the missing and weak boundaries, to deal with the intensity inhomogeneity, based on the LBF model, we introduced the convex total variation regularization term, for explicit smoothing of the level set function ø. The evolution equation will be solved through the level set method of calculus of variations. Results: In the experimental processing, we use some real images and magnetic resonance imaging brain images as the experimental images, to validate the stabilization of algorithm. The experimental results on comprehensive and sincerity images show the outstanding of our proposed model with reference to stabilization and availability. Conclusions: We propose a new segmentation local information of an image and introduce a new regularization functional is to keep the level set function smooth. Finally, various experimental results on real and lowcontrast image, showing which is a powerful type of images, including some that would be difficult to segment with gradientbased methods. In addition, the advantages of the proposed model are better than CV model and the LBF model. Our new model can effectively segment a real image.
Keywords: ChanVese model, image processing, image segmentation, level set method
How to cite this article: Li C, Su J, Yu L, Wang L, Ze L. A variational level set method image segmentation model with application to intensity inhomogene magnetic resonance imaging. Digit Med 2018;4:515 
How to cite this URL: Li C, Su J, Yu L, Wang L, Ze L. A variational level set method image segmentation model with application to intensity inhomogene magnetic resonance imaging. Digit Med [serial online] 2018 [cited 2020 Jun 1];4:515. Available from: http://www.digitmedicine.com/text.asp?2018/4/1/5/232708 
Introduction   
Segmentation is not only one of the extremely important assignments in image processing but also is a complex task and fundamental problem in the field of computer vision. The direct reason is that the performance of the subsequent processing steps would be affected by the segmentation result. A major difficulty for medical image segmentation is the intensity inhomogeneity, which is caused by many factors but I think the following reasons are more important, such as imperfections of imaging devices.^{[1],[2],[3],[4],[5]} The primary objective of image segmentation is to detect the edge of the objects in the given image. Image segmentation is often used to locate objects and boundaries (lines, curves, etc.) in an image. More precisely, image segmentation is a process of tagging each pixel in an image. At present, the method of image segmentation can be divided into two classifications: based on the curve methods and the method based on region. The methods based on curve include the Snake model and geometric active contour model proposed by Kass et al.^{[6]} and Caselles et al.,^{[7]} respectively. In the past decades, a lot of methods have been proposed for image segmentation. The wellknown models are as follows:^{[8],[9],[10],[11],[12],[13]} each of them has its own advantages and disadvantages. In practice, these kinds of models only can segment the images, which clearly defined objects. These models have another disadvantage that is very sensitive to initialization. The MumfordShan ^{[14]} and ChanVese (CV)^{[15]} model represent the method based on region, which has been successfully used in the binary phase segmentation with the assumption that each image region is statistically homogeneous. However, the PC model does not work well for medical images with intensity inhomogeneity.
Recently, a kind of models has been proposed to improve the performance in segmenting images with intensity inhomogeneity. Wang et al.^{[16]} also incorporated the multiscale segmentation idea into the level set formulation to segment inhomogeneous images.
Materials and Methods   
Methods
Let C be a simple closed curve in the image domain Ω, the curve can be continuously contracted and expanded at a certain speed along the normal direction, and the time variable t is introduced, a curve family with timevarying t is formed as C(t):
Curve family can be regarded as zero level set of highdimensional function ø. At this time, ø (x, t) is called a timevarying level set function, and a zero level set is a set whose function value is zero on the surface.
Given a level set function ø (C [t]), t), the ø (C [t]), t) = 0 is the zero level set of the evolution curve C(t) at the time t, according to the derivation rule of function, the partial derivative of time t, and the evolution equation of level set is obtained as follow:
Because the curve evolution is in accordance with the normal direction, then:
Where ∇ is the gradient operator, so we can get the basic equation of level set function such as follows:
In fact, the evolution of the level set function is a partial differential equation with the time change, describing the initial evolution curve of the level set function, application (4) to update the zero level set function, and realize the evolution curve of movement.
According to the theory of partial differential equations, solving (4), we need to give the definite conditions, generally speaking, the boundary conditions for Neumann boundary conditions, i.e., in the boundary on ∂Ω given boundary conditions:
The initial level set function is defined as the signed distance function as follows:
where d (C (t), x) refers to x to C (t) the Euclidean distance.
Finally, we can get the level set evolution formula like this:
Statistical analysis used
Let be an image, C (t) be a simple closed curve in the image domain for a tworegion segmentation, the curve C (t) defines a partition of the image domain Ω into two regions. R_{1}= R_{c}, which corresponding to the inside of the curve. And refers to the outside of the curve.
Generally speaking, various image segmentation models can be defined as the following general form:
Where D is a data fidelity term, which measures the statistical information of an image, and λ is a positive weight parameter, which balancing the contribution of the curve length term against the data term.
Most of the variational image segmentation algorithms optimize this functional. In practice, the region models usually are considered to a parametric distribution. In the proposed model, we assume each region obey normal distribution as follows:
where μ_{k} and σ_{k} denote, respectively, the mean and standard deviation of the image data information within R_{k} where k = 1, 2, 3, .... refers to each region.
In this section, we will give a brief review of some related models. Next, those methods can be recommended primitively as follows. The Snake or active contour model is put forward by Kass et al. The curve evolution is received by energy functional minimization. Let Ω be a boundary and open subset of R^{2}, with ∂Ω as its boundary. Let I(x) be the given image, as a boundary function defined on with real values. Usually, is a rectangular domain. Let be a parameterized curve. Then, the Snake method is minimizing the following energy functional:
The first two terms refer to the membrane energy and the elasticity energy, controlling the smoothness of the curve. They are called the internal energy. The third term is the external term and depends on the image data. It is easy to say that the external energy term is small, when the gradient of I has a large magnitude, thus pushing the curve toward edges. Such functions are usually called edge detectors. The active contour model was further developed by Caselles et al., Caselles and Coll, Kichenassamy et al., Sapiro, and Sandberg and Vese ^{[17],[18],[19],[20],[21]} using different edge detectors. For example, Caselles et al. introduced a geodesic active contour model by minimizing the functional energy:
where is a decreasing function such that:.
In the geodesic active contour model, is replaced by . In addition, considering that the Snake model does not allow topology change in the curve evolution, consequently can only detect the object in the image. Caselles et al.^{[7]} employed a level set representation, based on the pioneering work of Osher and Sethian.^{[22]} Let ϕ (t,.) be a level set function such that C(q) is the zero level set of ϕ, i.e., then the level set function would be evolved instead of the curve. Caselles proposed the following equation for ϕ:
where K is the curvature and c is a constant.
Mumford and Shah addressed an active contour model by minimizing the following energy functional:
where the boundary curve C is exactly the discontinuity set of u. CV formulated a piecewise constant variant of this model, and the boundary curve C was represented by a level set function, ϕ satisfying ϕ > 0 inside C and ϕ < 0 outside C. By defining the Heaviside function and onedimensional Dirac measure is , and the energy functional becomes:
The gradient decent equation for CV active contour model is:
with boundary condition .
The CV model uses global region information to detect objects. Its process is dependent of the initial position of the active contour. Basically, it is resistant against noise and can segment objects without apparent contour or texture. However, it often misclassifies objects into background when a lowcontrast or nonuniform illumination image is occurred. Assume that the parameters in (16) are λ_{1} = λ_{2} = 1.0 and μ = 0.005×255×255. The CV model has good performance in image segmentation due to its ability of obtaining a larger convergence, range and being less sensitive to the initialization. However, the CV model is only adapted for 2phase image. If the intensities with inside C are outside C are not homogeneous. The constant C_{1}and C_{2} will not be accurate. As a result, the CV model generally fails to segment images with intensity inhomogeneity. Similarly, more general piecewise constant models in a multiphase level set framework ^{[23]} are not good at such image either. To overcome the difficulty caused by intensity inhomogeneities, Li et al. proposed the local binary fitting (LBF) model,^{[11],[24],[25]} which can utilize the local intensity information. In the LBF model, two spatially varying fitting function F_{1}(x) and F_{2}(x) are introduced to approximate a given point x ∈ Ω, the local intensity fitting energy is defined by:
where g (z) is a Gaussian kernel function: .
F_{1}(x) and F_{2}(x) are two values that approximate image intensities inside and outside contour E_{x}(C,F_{1}(x),F_{2}(x)), which are defined for a center point x. For all the center point x in the image domain Ω, the energy functional can be defined by:
Inserting (17) into (18) and using the level set function ϕ to represent the contour C, then have:
To ensure stable evolution of level set function ϕ, Li et al. added the distance regularizing term in variational level set function, to penalize the distance function.^{[24],[25],[26],[27]} The deviation of the level set function ϕ from a signed distance function is characterized by the following integral:
To regularize the zero level contour of ϕ, they also need the length of the zero level set curve of ϕ, which is given by:
hence the new energy functional is designed to be as
The proposed model
We aim to solve the erroneous detection problem of LBF model and CV model when dealing with lowcontrast and nonuniform illumination images. Therefore, another regularization term, called global information function,^{[28],[29],[30],[31]} is proposed to detect the gradient amplitude at the outer region of the contour. Hence, we introduce the regularization functional as follows:^{[32]}
Where is an even and strictly convex function, nondecreasing on R^{+} and let φ^{∞} be the recession function φ of defined by, , then we set:^{[32]}
for details of the definition, referring to Aubert G and Kornprobst P's book.^{[32]} Moreover, there exists constants c > 0 and b > 0 such that:
where η ≥ 0 is a weight parameter. The aim of the term of the energy functional (25) is to keep the level set function smooth. Hence, now we define the entire energy functional as follows:
For the evolving curve C, we use an implicit representation given by the level set method C ∈ Ω, is represented implicitly by Lipschitz function , such that:
Let H is the Heaviside function and use unknown variable ϕ to replace the unknown variable C, then the energy functional's level set formulation of (18) can be expressed as:
In practice, the Heaviside function H in (28) is approximated by a smooth function H_{ε}defined:
The derivative of H_{ε}is the following smooth function:
By replacing H and δ in (28) with H_{ε} and δ_{ε}, thus, the energy functional E^{total} can be approximated by:
This is the energy functional, which we will minimize to find the object boundary.
By calculating variations,^{[32]} it can be shown that the function F_{1}(x) and F_{2}(x) minimize the energy functional in (31) satisfy the following EulerLagrange equations:
From (32) and (33), F_{1}(x) and F_{2}(x) are as follows:
To minimize the energy functional in (31) for a fixed ϕ, keeping F_{1}(x) and F_{2}(x) fixed, and minimizing the energy functional E^{total} (ϕ,F_{1}(x),F_{2}(x))with respect to ϕ, we solve the gradient flow equation as , where is the Gâteaux derivative ^{[32]} of the energy functional E (ϕ). We need to derive the Gâteaux derivative of the regularization term in (31) by calculating variations; it is straightforward to derive the Grateaux derivative of the energy E ^{REG} as:
Hence, we can get the gradient flow equation such as:
where e_{1} and e_{2} are the functions:
Numerical method
In this subsection, we briefly present the numerical scheme and procedure to solve evolution PDE (38) numerically:
To discrete the equation in ϕ, we use finite difference explicit schemes. We denote the space variable by (x,y) ∈ Ω, and the usual notations for finite differences. Let △x, △y be the steps of space, and fixed space steps are , and the time step △t = 0.1, and (x_{i},x_{j})=(ih,jh) be the grid points, for . Let with , the finite difference schemes are:
The curvature is discretized by computing the gradient of ϕ and its partial derivatives at the current iteration using central differences:
where operators △^{xy}, △^{xx},△^{yy} are applied to at i, j and iteration n.
The algorithm is as follows: first, we compute and using (34) and (35), respectively, then we compute by the following, discretization and linearization of the equation in ϕ:
Finally, the principal steps of the algorithm are as follows:
Step 1: Initialize the level set function
Step 2: Compute and using (34), (35), respectively.
Step 3: Compute and using (38) and (39), respectively.
Step 4: Solve discrete equation (43) to obtaining .
Step 5: Check whether the zero level set of ϕ is stationary, i.e., set a maximum iteration number N, a tolerance tol., and n = n + 1, If, , then stop the calculation. If not, set n = n + 1, and go to Step 2.
Convergence and stability analysis
Definition 1 BV (Ω)^{[32]} we define BV (Ω) the space of functions of bounded variation, as
Its norm defined like this: , so BV (Ω) is the Banach space.
Lemma 1 for arbitrary a given image , such that E^{New}(ϕ).
Proof: whitch satisfy the following growth conditions (25), so is bounded because the Heaviside is also bounded, we can get is bounded in L^{2} (Ω), namely:
represents the length of the level set function, so it also bounded, namely,
So we can easily get:
where, e_{1}, e_{2} are like (38), (39), and C is a constant. So,
Lemma 2,^{[32]} the regularized energy function:
So ∃ unique u ∈ BV (Ω) ∃ (46) hold.
Proof: “Existence” Let ϕ_{n}: is Minimization sequence of (46),
Meanwhile, .(46) uniformly boundness.
We can also proof that ϕ_{n} is also bounded in the L^{1}(Ω)space.
Let, , and then we can get , so Dv_{n}< M,
We can easily get ∃C > 0, such that the following inequalities hold true, namely:
Because , we can get , let , available through variable substitution, thereby,
So . Because of the norm of BV space is like this: So we can get:
Due to is ϕ_{n} both bounded in L^{1}(Ω), L^{2}(Ω). Furthermore, is bounded and ϕ_{n} is bouned in BV space. Therefore, to ϕ_{n}, it must exist u ∈BV (Ω), such that:
So we can get:
and,
So ϕ is the minimal point of the energy functional .
“Uniqueness”: If both of ϕ_{1}, ϕ_{2}are the minimal points of the energy functional , because is nondecreasing, strict convex, and even function therefore, Dϕ_{1}= Dϕ_{2}; furthermore, we can get , because is strict convex
Function, we can get ϕ_{1}= ϕ_{2}, so the solution is unique.
Theorem 1 Let so the energy functional (31) exists the minimum in BV space, and the minimum is unique.
Proof: Because the LemmaLemma 2, we can easily proof the Theorem 1 is true.
Results   
In this section, we show some numerical results for our image segmentation with the (44). We also test our proposed model on medical images. The first row is a vessel image with intensity inhomogeneity. The second row is an image of leaf with corners. The third row is different medical MR images. The forth row is the texture images. Like other papers, the stopping criterion is defined as follows:
where ϕ^{k} denotes the iteration of the scheme, and ε is a very small number. For most methods, ε = 10^{6} for the twophase segmentation and ε = 10^{4} for the multiphase segmentation.
To compare segmentation results quantitatively, we consider the segmentation accuracy (SA) defined the Dice coefficient,^{[33]}
Where S_{1} and S_{2} are twosegmentation accuracy. The results are presented in [Table 1].  Table 1: Segmentation accuracy (Dice coefficient) for the synthetic image
Click here to view 
Image of blood vessel
[Figure 1] shows the segmentation result for a real image of blood vessel with inhomogeneous intensity through using of the proposed numerical method. The computational domain set is Ω = (0, 1) × (0, 1) with a 64 × 64 mesh. The image is successfully segmented after 450 iterations.  Figure 1: Segmentation on vessel, which shows the segmentation result for a real image of blood vessel with inhomogeneous intensity through using of the proposed numerical method
Click here to view 
Image of with corners
[Figure 2] shows that our proposed model can segment object with corners successfully.  Figure 2: Segmentation on object with corners, which shows that our proposed model can segment object with corners successfully
Click here to view 
Images of brain magnetic resonance
[Figure 3] shows that our proposed model can be used to analyze medical images to provide necessary information for medical. The segmentation of brain MR image is shown on the computational domain Ω = (0, 1) × (0, 1) with a 256 × 256 mesh. As can be observed from [Figure 3], the agreement between the area of the brain, solid tumor, and the segmentation of image is good.  Figure 3: Segmentation on different magnetic resonance images, which shows that our proposed model can be used to analyze medical images to provide necessary information for medical
Click here to view 
The texture images
Texture images, based on local spatial variations of intensity or color to identify these types of homogeneous image regions, are an important attribute in analysis and pattern recognition. [Figure 4] shows that the proposed model can be very useful in detecting texture image segmentation. The computational domain is set to Ω = (0, 1) × (0, 1).  Figure 4: Segmentation on different texture images, texture images, based on local spatial variations of intensity or color to identify these types of homogeneous image regions, is an important attribute in analysis and pattern recognition
Click here to view 
The first and the third columns are the CV model's experimental effect images, the second and the forth columns are the LBF model's experimental effect images.
Discussion   
Settings and design
We test our method on some synthetic images, real MR medical images and some images with corners. The proposed model was implemented by MATLAB2016 on a computer with Intel Core2 Duo2.2GHz CPU, 2G RAM.
Material
We quantitatively compare our method with other two stateoftheart algorithms, namely, CV model and LBF model in segmenting synthetic MR images with the ground truth from BrainWeb.
The intensity nonuniformity is a common phenomenon which often occurs in image processing, especially medical images, and is caused by imperfect image acquisition,^{[34],[35],[36]} for example, radiofrequency field inhomogeneities appear due to either a nonuniform field itself or a nonuniform sensitivity of the receiver and transmitter coils. In medical images, artifacts lead to undesired intensity variations in the same tissue types across the image. When the segmentation algorithms rely only on image intensities, the artifacts affect successfully automated image segmentation.^{[37]}
Due to the problems above, so an efficient and reliable method deal with for intensity nonuniformity phenomenon is required. We propose to formulate the clustering energy functional with explicit level set function ϕ regularization because the property of function so the regularization term is strict convex, which results in the energy functional E^{total} exists and unique solution. In fact, we have presented the existence and uniqueness solution in the part of convergence and stability analysis. Where is the convex total variation ^{[38],[39]} of ϕ, and η is a constant, which refers to the weight between the data terms and the regularization term, so the convex total variation term serves for explicit smoothing of the level set function ϕ in the process of evolution.
Sensitivity to initialization: Due to our proposed model is a local image segmentation model, so the initial contour must be placed inside the object region but must not be closed to the true boundaries.
We also conduct comparisons of proposed method against CV model and LBF model, using the synthetic images shown in [Figure 5]. Comparisons with CV Model and LBF Model, the parameters for CV model are: λ_{1} = 1.0, λ_{2} = 1.0, μ = 0.005 × 255 × 255, ν = 1 and for LBF model are: λ_{1} = 1.0, λ_{2} = 1.0, μ = 0.001 × 255 × 255, σ = 3. The parameters for our proposed model are: λ_{1} = 1.0, λ_{2} = 1.0, μ =0.001 × 255 × 255, ν = 1, η = 2, σ = 3. From the discussion above, we can get the conclusion as follows: by introduce a new regularization term (this functional is to keep the let set function smooth), the proposed model can segment intensity inhomogene brain images, and other natural images; however, we can see from the numerical results, the proposed model is not suitable for segmenting the text image.  Figure 5: Segmentation on different images with ChanVese and local binary fitting model, the first and the third columns are the ChanVese model's experimental effect images, the second and the fourth columns are the local binary fitting model's experimental effect images. We also conduct comparisons of proposed method against ChanVese model and local binary fitting model, using the synthetic images shown in Figure 5. Comparisons with ChanVese Model and local binary fitting Model
Click here to view 
The main advantages of our model can be presented as follows. First, the proposed model is built the local intensity information in the formulation of the level set method. Second, a novel regularized energy objective functional is proposed by minimizing the images, which makes it more efficient in segmenting medical images with intensively inhomogeneity than the state of the art methods. Third, we demonstrate the performance of the proposed method as well as compare it to other stateoftheart about the SA. The method was tested on synthetic and real images and evaluated qualitatively and quantitatively. However, our method still exists some limitation, for example, it can be seen from the experimental results; our method is failed for the texture image. Moreover, the segmentation results are sensitive to the initialization.
Conclusions   
In this article, we utilize the Gaussian kernel function to detect the intensity difference on the term. We also propose a new segmentation local information of an image and introduce a new regularization functional is to keep the level set function smooth. Finally, various experimental results on real and lowcontrast image, showing which is a powerful type of images, including some that would be difficult to segment with gradientbased methods. In addition, the advantages of the proposed model are better than CV model and the LBF model. Our new model can effectively segment a real image.
Acknowledgment
Funding was provided by the Natural Science Foundation of China under Grant NO.61361126011, No. 90912006; the Special Project of Informatization of Chinese Academy of Sciences in “the Twelfth FiveYear Plan” under Grant No. XXH12504106. The authors wish to gratefully thank all anonymous reviewers who provided insightful and helpful comments.
Financial support and sponsorship
Funding was provided by the Natural Science Foundation of China under Grant NO.61361126011, No. 90912006;the Special Project of Informatization of Chinese Academy of Sciences in “the Twelfth FiveYear Plan” under Grant No. XXH12504106.
Conflicts of interest
There are no conflicts of interest.
References   
1.  Wong WC, Chung AC. Bayesian image segmentation using local isointensity structural orientation. IEEE Trans Image Process 2005;14:151223. 
2.  Chen Y, Zhang J, Yang J. An anisotropic images segmentation and bias correction method. Magn Reson Imaging 2012;30:8595. 
3.  Awate SP, Tasdizen T, Foster N, Whitaker RT. Adaptive markov modeling for mutualinformationbased, unsupervised MRI braintissue classification. Med Image Anal 2006;10:72639. 
4.  Zhou S, Wang J, Zhang S, Liang Y, Gong Y. Active contour model based on local and global intensity information for medical image segmentation. Neurocomputing 2016;186:10718. 
5.  Zhou S, Wang J, Zhang M, Cai Q, Gong Y. Correntropybased level set method for medical image segmentation and bias correction. Neurocomputing 2017;234:21629. 
6.  Kass M, Witkin A, Terzopoulos D. Snakes: Active contour models. Int J Comput Vision 1988;1:32131. 
7.  Caselles V, Kimmel R, Sapiro G. Geodesic active contours. Int J Comput Vision 1997;22:6179. 
8.  Chen Y, Zhang J, Mishra A, Yang J. Image segmentation and bias correction via an improved level set method. Neurocomputing 2011;74:352030. 
9.  Chen Y, Zhang J, Macione J. An improved level set method for brain MR images segmentation and bias correction. Comput Med Imaging Graph 2009;33:5109. 
10.  Zhan T, Zhang J, Xiao L, Chen Y, Wei Z. An improved variational level set method for MR image segmentation and bias field correction. Magn Reson Imaging 2013;31:43947. 
11.  Li C, Huang R, Ding Z, Gatenby JC, Metaxas DN, Gore JC, et al. Alevel set method for image segmentation in the presence of intensity inhomogeneities with application to MRI. IEEE Trans Image Process 2011;20:200716. 
12.  Li C, Huang R, Ding Z, Gatenby C, Metaxas D, Gore J, et al. Avariational level set approach to segmentation and bias correction of images with intensity inhomogeneity. Med Image Comput Comput Assist Interv 2008;11:108391. 
13.  Lewis EB, Fox NC. Correction of differential intensity inhomogeneity in longitudinal MR images. Neuroimage 2004;23:7583. 
14.  Mumford D, Shah J. Optimal approximation by piecewise smooth function and associated variational problems. Commun Pure Appl Math 1989;425:77685. 
15.  Chan T, Vese L. An Active Contour Model without Edges, ScaleSpace Theories in Computer Vision. Second International Conference; 1999. 
16.  Wang XF, Min H, Zou L, Zhang YG, Tang YY, Chen CL. An efficient level set method based on multiscale image segmentation and hermite differential operator. Neurocomputing 2016;188:90101. 
17.  Caselles V, Catte F, Coll T, Dibos F. A geometric model for active contours in image processing. Numer Math 1993;66:131. 
18.  Caselles V, Coll B. Snakes in movement. SIAM J Numer Analy 1996;33:244556. 
19.  Kichenassamy S, Kumar A, Olver P, Tannenbaum A, Yezzi A. Gradient flows and geometric active contour models. Proceedings of the Fifth International Conference on Computer Vision, ICCV; 1995. p. 8105. 
20.  Sapiro. Vector snakes: A geometric framework for color, texture and multiscale image segmentation. Proceedings IEEE International Conference on Image Processing; 1996. 
21.  Sandberg B, Vese L. Active contours without edges for vector valued images. J Vis Commun Image Represent 2000;11:13041. 
22.  Osher S, Sethian JA. Fronts propagating with curvaturedependent speed: Algorithms based on HamiltonJacobi formulations. J Comput Phys 1988;79:1249. 
23.  Vese LA, Chan TF. A multiphase levelset framework for image segmentation using the munford and Shahn model. Int J Comput Vis 2002;50:27193. 
24.  Li CM, Kao CY, Gore JC, Ding ZH. Implicit active contours driven by local binary fitting energy. IEEE Conference Computer Vission and Pattern Recognition; 2007. p. 17. 
25.  Li C, Kao CY, Gore JC, Ding Z. Minimization of regionscalable fitting energy for image segmentation. IEEE Trans Image Process 2008;17:19409. 
26.  Li CM, Xu CY, Gui CF, Fox MD. Levelset evolution without reinitialization: A new variational formulation. IEEE Computer Society Conference Computer Vision and Pattern Recognition; 2005. p. 4306. 
27.  Zhang K, Song H, Zhang L. Active contours driven by local image fitting energy. Pattern Recogn 2010;43:1199206. 
28.  Chan TF, Esedoglu S, Nikolova M. Algorithms for finding global minimizing of image segmentation and denoising models. SIAM J Appl Math 2006;66:163248. 
29.  Chambolle A. An algorithm for total variation minimization and applications. J Math Imaging Vision 2004;19:8997. 
30.  Hua G, Liu Z, Zhang Z, Wu Y. Iterative localglobal energy minimization for automatic extraction of objects of interest. IEEE Trans Pattern Anal Mach Intell 2006;28:17016. 
31.  Lankton S, Tannenbaum A. Localizing regionbased active contours. IEEE Trans Image Process 2008;17:202939. 
32.  Aubert G, Kornprobst P. Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations. Applied Mathematical Sciences; 2 ^{nd} ed., Vol. 147. New York: Springer; 2006. 
33.  Dice LR. Measures of the amount of ecologic association between species. Ecology 1945;26:297302. 
34.  Ivanovska T, Laqua R, Wang L, Schenk A, Yoon JH, Hegenscheid K, et al. An efficient level set method for simultaneous intensity inhomogeneity correction and segmentation of MR images. Comput Med Imaging Graph 2016;48:920. 
35.  Wang L, Shi F, Yap PT, Lin W, Gilmore JH, Shen D, et al. Longitudinally guided level sets for consistent tissue segmentation of neonates. Hum Brain Mapp 2013;34:95672. 
36.  Zhou Y, Shi WR, Chen W, Chen YL, Li Y, Tan LW, et al. Active contours driven by localizing region and edgebased intensity fitting energy with application to segmentation of the left ventricle in cardiac CT images. Neurocomputing 2015;156:199210. 
37.  Gonzalez R, Woods RE. Digital Image Processing. Englewood Cliffs: Prentice Hall International; 2002. 
38.  Jung M, Kang M. Efficient nonsmooth nonconvex optimization for image restoration and segmentation. J Sci Comput 2015;62:33670. 
39.  Kang M, Kang M, JungM. Total generalized variation based denoising models for ultrasound images. J Sci Comput 2017;72:17297. 
[Figure 1], [Figure 2], [Figure 3], [Figure 4], [Figure 5]
[Table 1]
