{'fraction_non_alphanumeric': 0.0411759073063296, 'fraction_numerical': 0.013980656899358422, 'mean_word_length': 4.51487414187643, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Local Active Contour Model for Image Segmentation with Intensity Inhomogeneity

A novel locally statistical active contour model (ACM) for image segmentation in the presence of intensity inhomogeneity is presented in this paper. The inhomogeneous objects are modeled as Gaussian distributions of different means and variances, and a moving window is used to map the original image into another domain, where the intensity distributions of inhomogeneous objects are still Gaussian but are better separated. The means of the Gaussian distributions in the transformed domain can be adaptively estimated by multiplying a bias field with the original signal within the window. A statistical energy functional is then defined for each local region, which combines the bias field, the level set function, and the constant approximating the true signal of the corresponding object. Experiments on both synthetic and real images demonstrate the superiority of our proposed algorithm to state-of-the-art and representative methods.


I. INTRODUCTION
Image segmentation is an important procedure in many computer vision and pattern recognition applications.
The edge-based ACM methods are applicable to images with intensity inhomogeneity but, in general, they are sensitive to the initialization of the level set function. Moreover, they easily suffer from serious boundary leakage for images with weak boundaries [19]. Many region-based ACMs are based on the global region information; they assume that the image intensity is homogeneous [13] [22], and thus are not suitable for segmenting images with intensity inhomogeneity.
Intensity inhomogeneity caused by the imperfection of imaging devices or by illumination variations often occurs in real-world images, and it can lead to serious misclassifications by intensity-based segmentation algorithms that assume a uniform intensity [13] [22]. Statistically, misclassifications are caused by the prolonged tails of the intensity distribution of each object so that it is difficult to extract the desired objects accurately based on their respective intensity distributions (refer to Fig. 1 please). The well-known Mumford-Shah (MS) model [12], which assumes that an image is piecewise smooth, is suitable for modeling images with intensity inhomogeneity. The MS model uses a set of contours C to separate different regions.
However, it is difficult to minimize the energy functional of the MS model because the set C of low dimension is unknown and the problem is non-convex [14]. Some simplified versions of the MS model have been proposed, such as the seminal Chan-Vese (CV) model [13] and the piecewise-smoothing (PS) model [14][15]. All these methods represent the contour C using the zero level of a function called the level set function, and then segmentation proceeds by evolving a level set evolution equation, obtained by minimizing some energy functional. However, the CV model is not applicable to images with intensity inhomogeneity because it models images by means of the piecewise constant. The PS model can yield a desirable segmentation for images with intensity inhomogeneity. However, the iteration of two partial differential equations (PDE) is needed in the PS model to approximate the original image, and this limits the practical application of PS model because it is very time-consuming.
Recently, some local region-based ACMs have been proposed for images with intensity inhomogeneity, such as the local region descriptors (LRD) method [20], the local region based (LRB) method [21], the local binary fitting (LBF) model [17][18], the local intensity clustering (LIC) method [19], and the local region model (LRM) [41], etc. However, there exist some drawbacks with these local region based models. The LRD model needs to tune a balloon force parameter for images with strongly overlapping intensity distribution, but how to define the degree of the overlap is not mentioned in [20]. Therefore, it is difficult to select a proper balloon force. The LRB method has two drawbacks. First, the Dirac functional used there is restricted to a neighborhood around the zero level set, which makes the level set evolution act locally. Therefore, the evolution can be easily trapped in local minima [13]. Second, the region descriptor in LRB is only based on the region mean information without considering the region variance, and this may lead to inaccurate segmentation. The second drawback of LRB also holds for the LBF model, because they use a similar energy functional. The LIC method can be considered as a locally weighted K-means clustering method [19]. It does not consider the clustering variance, which may cause inaccurate segmentation; similar drawback exists for the K-means clustering based method in [26]. The LRM method relates the local region statistics, i.e., local region means and variances, in interpreting the MS model. However, the local region means and variances are only defined empirically, but not derived from minimizing the MS energy.
Intensity inhomogeneity is usually ascribed to a smooth and spatially varying field multiplying the constant true signal of the same object in the measured image [11] [19]. This spatially varying and smooth field is called a bias field. This paper presents a novel statistical ACM for simultaneous segmentation and bias correction. By exploiting the image's local redundant information, we define a mapping from the original image domain to another domain such that the intensity probability model is more robust to noise and the overlapping of intensity is suppressed to some extent. We then devise a statistical energy functional for the distribution of each local region in the transformed domain, which combines the bias field, the level set function, and the constant approximating the signal of the corresponding object. Analysis of the proposed approach shows that it is a soft classification model, which means that each pixel can be assigned to more than one class. In contrast, the hard classification used in previous methods [13] [45] assigns each pixel to only one class. Therefore, the proposed approach can achieve a better segmentation result. In addition, the proposed method can be applied to simultaneous tissue segmentation and bias correction for magnetic resonance (MR) images. Our preliminary work on this has been briefly presented in [23].
The rest of the paper is organized as follows. Section II introduces the background of our research, including the seminal MS model, the CV model, and the PS model. Section III describes the statistical model for intensity inhomogeneity and our proposed model in detail. Section IV presents the proposed algorithm.
Section V shows extensive experimental results. In Section VI, we discuss the relationships between our method and the CV model, the LBF model and the LIC model. Section VII concludes the paper.

II. BACKGROUND
Let Ω ⊂ R N , where N = 2 or 3, be the image domain and I(x): Ω→R be an input image. Mumford and Shah [12] approximated an image with a piecewise-smooth function u(x): Ω→R, such that u varies smoothly within each of the sub-regions, and abruptly across the boundaries of the sub-regions. Let C(p): R→Ω approximate the edges of the sub-regions. The energy functional is given as follows [12]: where μ, ν > 0 are two fixed parameters, and |C| represents the length of the contour. The image segmentation can be performed by minimizing Eq.(1) with respect to u and C. However, it is difficult to minimize the above functional in practice, due to the unknown set C of lower dimension and the non-convexity of the functional [14]. Many methods have been proposed to simplify or modify the functional [13][14][15], which will be reviewed in the following. Here we mainly discuss two-phase cases, where the region Ω is separated by a contour C and The results can be readily extended to multi-phase cases by utilizing two different level set functions to represent different regions (refer to [14] for details).
Chan and Vese [13] proposed an active contour model which can be seen as a particular case of Eq.(1), and they represented u(x): Ω→R in Eq.(1) as a piecewise-constant function. Chan and Vese proposed to minimize the following energy functional (i.e., the CV model) where c 1 and c 2 are two constant functions which approximate the average intensities inside and outside the contour C, respectively. This energy functional can be represented by a level set formulation, and by evolving the level set, we can obtain its minimum.
The piecewise-smooth (PS) models [14][15] where H(•) is a Heaviside functional defined as , u + and u − must be obtained by solving the two damped Poisson equations before each iteration of the level set function, and the computational cost is very high. Moreover, u + and u − must be extended to the whole image domain, which is very difficult to implement in practice. More details about the implementation of the PS model can be found in [14].
The CV and PS energy functionals are non-convex with respect to the level set function φ, and therefore they may fall into the local minima [29]. The global CV and PS models have been proposed in [28] by adding a total variation term into the energy functional and by utilizing Chambolle's method [30] to yield the global minimizer. However, the global CV and PS models, like the original CV and PS models, still have some drawbacks, such as failures to segment an image with intensity inhomogeneity for the CV model and large computational cost for the PS model, because the two smooth functions u + (x) and u − (x) still need to be iteratively solved.

A. Statistical Model of Intensity Inhomogeneity
Let Ω be the image domain, b(x): Ω→R be an unknown bias field, I(x): Ω→R be the given image, J(x): Ω→R be the true signal to be restored, and n(x): Ω→R be noise. We consider the following model of intensity inhomogeneity [11] [19] where μ i (x) is the spatially varying mean, and σ i is the standard deviation. Since b(x) varies slowly, it can be assumed to be a constant in a small window [19]. Thus we can assume that , N} in our following discussions. As shown by the blue solid curves in Fig. 1, when I(x) falls into the tail of the distribution, misclassification will occur. One approach to reducing the overlapping tail is to compress the profile of the distribution while keeping its mean fixed [22] (see the red dashed curves in Fig. 1).

B. Principle of the Proposed Method
For each position x in the image domain, we denote by Ο x a neighboring region center, i.e., Ο x = {y | ||y − x|| ≤ ρ}, whereρ is the radius of the region Ο x . Since there are N non-overlapping objects in the image with Ω i being the domain of the i th object, the whole image domain Ω can be represented as to another domain R(T) as follows: where m i (x) = ||Ω i ∩O x ||. The intensity of pixel x is assumed to be independently distributed [22]. Thus, whereby the corresponding probability density function (PDF) is still a Gaussian [27], i.e., Referring to the red dashed curves in Fig. 1, the overlapping tails of the distributions are suppressed to some extent.
Since the intensity inhomogeneity manifests itself as a smooth intensity variation across an image [19], Therefore, Let D = {I (x|θ i ), x∈Ω, i = 1,…,N}, we have the following likelihood function for the i th object [27]: We construct the following joint likelihood function Put Eq.(7) into Eq.(11), then Eq. (11) can be re-written as a univariate Gaussian distribution as follows: Obviously, the joint likelihood function in Eq.(10) of each pixel (voxel) is composed of multiple classes of intensities; thus, by using Eq.(13), our model can yield a soft classification, satisfying the condition of the partial volume effect [33] (i.e., the intensity of each volume voxels is mixed from multiple classes [33]).
Moreover, as can be seen from Eq. (6), the intensity in the transformed domain exploits the information about neighboring pixels belonging to the same class, so its classification result is less sensitive to noise and can result in a smoother object border.
Using Eqs. (5) and (15), and eliminating the trivial constant term, l(θ) can be re-written as

C. Energy Functional Formulation using the Level Set Method
One level set function φ can only represent two regions, inside and outside the contour C, as Ω + = inside(C) = {φ > 0} and Ω − = outside(C) = {φ < 0}, respectively. This is called the Two-Phase model. If there exist more than two different regions, two level set functions, φ 1 and φ 2 , can be used to represent different regions based on the Four-Color Theorem [14] such that any two adjacent regions can be represented by different colors. This is called the Four-Phase model. We define the phase indicators as follows: where H(•) is the Heaviside functional, Φ represents the set of the level set functions such that Φ={φ} for the Two-Phase model, and Φ={φ 1 ,φ 2 } for the Four-Phase model. The energy functional l(θ) can then be re-written
We only consider the Two-Phase level set method in the following discussions. Similar discussions can be readily extended to the Four-Phase case. Let's revisit the data-fitting term of our proposed energy functional in Eq. (18). When the variance This data-fitting term is similar to that of the well-known CV model [13] except for the trivial constant πρ 2 /2.
Thus, our energy functional is a generalized version of the CV model, which uses two constants c 1 and c 2 which are the average intensities in different regions to approximate the image intensity in the regions Ω + ={φ(x)>0} and Ω − ={φ(x)<0}, respectively.
However, the average intensity in a large region cannot fit the intensity well when the intensity of the image is inhomogeneous. The boundary information can handle inhomogeneous intensity well. The GAR model [48] [49] combines the region and boundary information into an energy functional. Different from the CV model only using the mean of region in Eq. (19), the region part in GAR model also considers the variance of region.
LRM [41] can handle the inhomogeneous intensity by using spatially varying means and variances to replace the constant mean and variance. However, LRM directly introduces a Gaussian kernel to compute the varying means and variances which is inconsistent with theory as pointed by [46].
The data-fitting term of the energy functional in the LBF model is as follows: ( ) where K σ (•) is a truncated Gaussian kernel with standard deviation σ, which satisfies ∫ K σ (x)dx=1, and f 1 and LBF data we can obtain that f 1 and f 2 are the weighted average image intensities in a Gaussian window inside and outside the contour, respectively. This is why the LBF model can handle images with intensity inhomogeneity well.
It is worth noting that our proposed model can also be seen as a generalization of the LBF model [17] [18] and LIC model [19] [47]. However, their intrinsic principles are different. If we set 2 Eq. (18), and use a truncated Gaussian kernel K ρ with standard deviation ρ satisfying ∫ K ρ (x)dx=1 to replace the original constant kernel, then Eq. (18) can also be written as This is the same as 1 2 ( , , ) LBF data E f f φ , except for the trivial constant 1/2. If we set f i = bc i , i = 1, 2 in Eq. (21), it is the energy functional of the LIC model [19] [47]. It should be noted that in the theoretical analysis of our method in Section III.B, the constant kernel is reasonably used as a local region indicator. However, the LBF model and LIC model use Gaussian kernel as the locally spatially weighted function to relate the pixel x and its neighboring pixel y. The closer y is to x, the larger the weight is assigned, representing the higher similarity between the intensities of pixels y and x. Therefore, the principles among our method, the LBF model [17] and the LIC model [19] are different.
The energy functional of LGD model [46] is slightly different from the LBF model which replaces the L2-norm terms in Eq. However, the Gaussian kernel is used in LGD model whose physical meaning is similar to the LBF model as we discussed above. Third, our model can be used for simultaneous segmentation and bias correction while the LGD model can be only used for segmentation.

V. ENERGY MINIMIZATION AND LEVEL SET EVOLUTION FORMULATIONS
The minimization of l(θ,Φ) with respect to (w.r.t.) each variable in θ = {c i , b,σ i , i = 1,…,N, N = 2, or 4} can be obtained by fixing other variables, yielding the closed forms of solutions, as described below.

A. Closed Form Solutions for Different Variables
Minimization w.r.t. c i . By fixing the other variables in Eq. (18), we obtain the minimizer of c i , denoted by i c , as follows: where * denotes the convolution operator.
Note that b is actually the normalized convolution [34], which naturally leads to a smooth approximation of the bias field b.
For an explanation of how to derive the above solutions, please refer to Appendices A and B.

B. Two-Phase Level Set Evolution Formulation
Minimizing the energy functional ( , ) l Φ θ w.r.t. φ, we have the corresponding gradient descent formulation as follows: where δ(φ) is the Dirac functional.
In order to keep the numerical implementation stable, the level set function should be regularized during the iteration of Eq.(25). Li et al. [35] proposed a signed distance-regularization formulation to regularize the level set function. However, as indicated by Xie [37], Li et al.'s method can produce some unnecessary valleys and peaks, which makes the level evolution easy to fall into some local minima. Li et al. [35] also proposed another improved signed distance-regularization formulation, but we found that the level set evolution is unstable if we use it in our formulation. The reason may be that we use different approximations for the Heaviside functional and the Dirac functional. In this paper we propose a new and simple method to regularize the level set function during iteration, which can make the evolution stable. After each iteration of level set evolution, we diffuse the level set function using the following formulation: 1 2 2 n n n t φ φ φ + = +Δ ⋅∇ (26) where φ n represents the level set function yielded by Eq. (25) during the n th iteration, ∇ 2 represents the Laplacian operator, and Δt 2 represents the diffusion strength. The φ n+1 in Eq. (26) can also be approximated by φ n+1 =K * φ n , where K is either a Gaussian kernel [24] or a constant kernel [23].

C. Four-Phase Level Set Evolution Formulation
Minimizing the energy functional ( , ) l Φ θ with respect to φ 1 and φ 2 , respectively, we can have the corresponding gradient-descent formulations as follows: Similar to the Two-Phase case, the level set functions are regularized by the following formula after each iteration:

D. Numerical Implementation
In the level set evolution of Eqs. (25) and (27), we only need to approximate the temporal derivative as a forward difference because there are no partial derivatives. The Laplacian operator ∇ 2 φ is approximated by where * is a convolution operator, and K is a kernel defined as K= 0 1 0 [37]. Therefore, the solution of the diffusion equation Eq.(26) can be discretized as follows: The standard Von Neumann analysis [43] can be used to analyze the stability for the time step Δt 2 . Putting Therefore, we have 1-8Δt 2 ≤ r ≤1. By solving the inequality |1-8Δt 2 |≤1, we have: The Heaviside functional H(z) is approximated by a smooth function H ε (z) as where ε > 0 is a constant. The Dirac functional δ(z) is approximated by δ ε (z) as follows The profile of δ ε (z) has a large support, which acts on all level curves of the level set function. This can allow the level set evolution to reach a global minimum easily [13].
Finally, based on the above description of our algorithm, the procedures for our proposed algorithm are summarized in Table I.

VI. EXPERIMENTAL RESULTS
In this section, we compare our method with the CV model [13], the global CV (GCV) model [28] (the code was downloaded from [39]), the LRB model [21] (the code was downloaded from [40]), the LIC model [19], and the LBF model [17] [18]  , for which we set Δt 2 = 0.01. Our method is stable for a wide range of ρ, e.g., 5< ρ < 25. In most cases, we set ρ = 6. A small ρ will make the computation in each iteration more efficient, but the convergence is slower. On the other hand, a large ρ will increase the computational burden in each iteration, while the convergence rate can be increased because information from larger regions is used. Therefore, the total computation burden is comparable for different ρ.
In Section V-A, we will compare our method with the CV model and the GCV model for segmenting two synthetic images and two real vessel images with intensity inhomogeneity. In Section V-B, we compare our method with the LBF model, the LIC model and the LRB model for segmenting a synthetic image and two real MR brain images with severe intensity inhomogeneity. Sections V-A and V-B focus on the Two-Phase level set model, and in Section V-C we test our method for simultaneous segmentation and bias correction.
We first use a synthetic image with the ground-truth bias field to qualitatively evaluate our method in comparison with the LIC model. Then, we compare our method with the LIC model on a real 3T MR brain image with different initial contours, and apply our method to two 7T MR brain images. At last, in Section V-D, we quantitatively compare our method with the CV model, the GCV model, the LRB model, the LBF model, and the LIC model by segmenting five synthetic images with gradually increasing strength of intensity inhomogeneity.

A. Comparisons with the CV Model [13] and the GCV Model [28]
Both the CV model [13] and the GCV model [28] assume that the image intensity is piecewise constant and use the global intensity means to fit the image intensity. Therefore, they do not perform well in images with intensity inhomogeneity. In this section, we compare our method with the CV model and the GCV model by applying them to some real and synthetic images with intensity inhomogeneity in order to demonstrate the superior performance of our method to them for images with intensity inhomogeneity. Fig. 2 demonstrates the segmentation results on two synthetic images and two real vessel images with intensity inhomogeneity by using the CV model [13], the GCV model [28], and our method, respectively. It can be seen that our method yields satisfactory segmentation results because we consider and exploit the image local region information, which can better separate the object from background. For the CV model and the GCV model, they use the global region information for segmentation and thus result in severe misclassification on these images because there exist severe overlaps of intensity between the objects and background. Using the global mean information cannot discriminate the overlapping intensity. . From top row to bottom row: segmentation results by the CV model [13], the GCV model [28] and our method. The blue circle represents the initial contour and the red lines represent the final segmentation contour.
We set ρ = 12 for the left image and ρ = 6 for the other images. [19] and the LRB Model [21] In this section, we use one synthetic image and two real MR images with severe intensity inhomogeneity to demonstrate the superior performance of our method to the LBF model, the LIC model and the LRB model, which all use local means to fit the image intensity and perform well in images with light intensity inhomogeneity. The LBF model [17][18] uses the local intensity mean to fit the measured image, and thus it can yield better segmentation results than CV and GCV models on images with intensity inhomogeneity.

B. Comparisons with the LBF Model [17][18], LIC Model
However, if the intensity inhomogeneity is severe, using only the local mean information may fail to discriminate the intensity between an object and its background, leading to inaccurate segmentation. Some experimental results are shown in the second column of Fig. 3. Similar drawbacks exist for the LIC model [19] (refer to the third column of Fig. 3) and LRB model [21]. Moreover, the localized version of the Heaviside functional and the Dirac functional used by the LRB model makes the level set evolution easily fall into local minima [13], as can be seen in the fourth column of Fig. 3. Our method can produce much better segmentation results because it considers the statistical information in a transformed domain, where the intensity of the object and background is less overlapping than in the original domain, making our method have a very strong discriminative capability for the object and background (refer to the left column of Fig. 3).

Figure 3:
Segmentation on a synthetic image (1 st row) and two real MR images (2 nd and 3 rd rows) with severe intensity inhomogeneity. From left to right: segmentation results using our method, the LBF model [17][18] , the LIC model [19] and the LRB model [21], and the histograms of image intensity. The blue circles represent the initial contours, and the right lines represent the final segmentation. We set ρ = 6 for the first two images and ρ = 20 for the last image.

C. Application to Simultaneous Segmentation and Bias Correction
In this section, we apply our ACM for simultaneous segmentation and bias correction, especially for MR images. We compare our method with the state-of-the-art LIC model [19] because only the LIC model is applicable for simultaneous segmentation and bias correction while the other models (i.e., the CV model, the GCV model, the LBF model and the LRB model) can only be used for segmentation.  [19], the ground-truth bias field, the estimated bias field using our method, and the estimated bias field using the LIC model [19]. We set ρ = 10 for the experiments.  Because one level set function can only represent two classes of tissues, we need to evolve two level set functions φ 1 and φ 2 according to Eq. (24) for four classes of tissues [14]. It can be seen that the tissue segmentation results by our method (the 1 st and 3 rd rows in Fig. 5) are much more accurate than those by the LIC method (the 2 nd and 4 th rows in Fig. 5). It should be noted that it is very easy to initialize the level set functions in our method. The initial contours can be set inside, outside or across the object boundary. The two initializations in the 1 st and 3 rd rows are very different, but the final segmentation results by our method are very similar. This demonstrates the robustness of our method to initialization.
Finally, we show in Fig. 6 the bias-correction results using our method on two 7T MRI images. This experiment aims to obtain the bias corrected MRI images. The original images, estimated bias fields, and the bias-corrected images are shown in the left, middle, and right columns, respectively. It can be clearly seen that the image quality is significantly improved by our method. Some regions (inside the red circles) whose intensity contrast is too low to be identified are able to be distinguished clearly after the bias correction.  GCV [21] LRB [14] LBF [11] LIC [12] Figure 7: Quantitative comparisons among our method, the CV model [13], the GCV model [28], the LRB model [21], the LBF model [17] [18], and the LIC model [19]. Top row: tested images, where the strength of intensity inhomegeneity is gradually increased from left to right. Bottom row: the corresponding JS values yielded by the competing methods on the five images. We set ρ =10.5 for all the experiments. Figure 8. From left to right: segmentation results on image with strong intensity inhomegeneity by our method, the CV model [13], the GCV model [28], the LRB model [21], the LBF model [17] [18], and the LIC model [19]. We set ρ = 10.5 for all the experiments.

D. Quantitative Evaluation
In this section, we use the Jaccard similarity (JS) [44] as an index to measure the segmentation accuracy and to quantitatively evaluate the segmentation performance of competing methods on images whose intensity inhomogeneity has different strength. The JS index between two regions S 1 and S 2 is calculated as JS(S 1 ,S 2 ) = |S 1 ∩S 2 |/|S 1 ∪S 2 |, which is the ratio between the intersectional area of S 1 and S 2 and their union. Obviously, the closer the JS value is to 1, the more similar S 1 is to S 2 . In our experiments, S 1 is the segmented object region produced by the six competing methods (the CV model [13], the GCV model [28], the LRB model [21], the LBF model [17][18], the LIC model [19], and our proposed method), and S 2 is the ground-truth.
We test the six competing methods on five synthetic images with different intensity inhomogeneity. The JS values of the six methods are shown in the bottom image of Fig. 7. Obviously, the JS values obtained by our method have little difference for intensity inhomogeneity with different strength, which demonstrates that it is very robust to image intensity inhomogeneity. For the global region-based CV model and GCV model, when the strength of intensity inhomogeneity is not strong (refer to images 1 and 2 in Fig. 7), both of them can yield a high JS value. However, when the strength of inhomogeneity is not low (refer to images 3, 4, and 5), the performance of the CV model and the GCV model degenerates rapidly. For the local region-based models (i.e., the LRB model, the LBF model and the LIC model), when the strength of intensity inhomegeneity is strong (refer to images 4 and 5), the performance of these methods also degenerates severely.
This is because using only the local region means cannot discriminate the object and the background satisfactorily when the intensity between the object and background overlaps severely. Our method gives consistently the highest JS values because it pursues the segmentation in a transformed domain where the intensity overlapping is much suppressed. The final segmentation results for the top right-most image in Fig.   7 by the six competing methods are shown in Fig. 8.

E. Robustness to Initializations and Region Scale Parameters
We use the tested image in Fig. 8 to demonstrate the robustness of our method to different level set initializations and region scale parameter ρ. We again use the JS index to measure segmentation accuracy.
The top row in Fig. 9 demonstrates the segmentation results using different initial contours (the blue lines).
We can see that there are no obvious visual differences for these segmentation results.

VII. CONCLUSION
In this paper, we presented a locally statistical active contour model for segmenting images with intensity inhomogeneity. Our method combines the information about neighboring pixels belonging to the same class, which makes it strong to separate the desired objects from background. Moreover, the proposed method yields a soft segmentation, which can, to some extent, satisfy the condition of the partial volume effect. In addition, the segmentation results are very robust to the initialization of the level set function, making it useful for automatic applications. Comparisons with five representative and state-of-the-art methods on synthetic and real images have demonstrated the effectiveness and the advantages of the proposed algorithm.