A Variational Level Set Approach Based on Local Entropy for Image Segmentation and Bias Field Correction

Image segmentation has always been a considerable challenge in image analysis and understanding due to the intensity inhomogeneity, which is also commonly known as bias field. In this paper, we present a novel region-based approach based on local entropy for segmenting images and estimating the bias field simultaneously. Firstly, a local Gaussian distribution fitting (LGDF) energy function is defined as a weighted energy integral, where the weight is local entropy derived from a grey level distribution of local image. The means of this objective function have a multiplicative factor that estimates the bias field in the transformed domain. Then, the bias field prior is fully used. Therefore, our model can estimate the bias field more accurately. Finally, minimization of this energy function with a level set regularization term, image segmentation, and bias field estimation can be achieved. Experiments on images of various modalities demonstrated the superior performance of the proposed method when compared with other state-of-the-art approaches.


Introduction
Image segmentation has always been a crucial step in image understanding and computer vision. However, due to the limitations in imaging instrument and other external effects, bias field is often seen in many medical images which can be ascribed to a spatially varying field. Therefore, it is often a fundamental step to correct the bias field before performing quantitative analysis of the image data. Methods of bias field correction can be divided into two groups: prospective methods [1,2] and retrospective methods [3,4]. Prospective methods aim to avoid intensity inhomogeneity by using the shim techniques, special imaging sequences, or special hardware. However, the retrospective methods rely exclusively on the information of the acquired images and thus can extract information on intensity inhomogeneity.
Among all kinds of bias field correction methods, those based on segmentation are the most popular. In these methods, the tasks of segmentation and bias field correction are interleaved in an iterative process, thereby allowing their optimal results to be simultaneously achieved. In [5], Wells III et al. proposed an expectation-maximization (EM) algorithm for image segmentation and bias field estimation. However, such algorithm is sensitive to the choices of initial conditions, which limits its applications in automatic segmentation. Based on the EM algorithm, Meena and Shantha [6] presented a fully automated method for MR brain image segmentation by introducing the Fuzzy C-means (FCM) with spatial information. Also, a method of pixel relabeling is included to improve the segmentation accuracy. So, this method can very well be extended to segmentation of clinical MR brain images and the identification of pathologies. In [7], Xie et al. propose an approach joining a modified MRF classification and bias field estimation in an energy minimization framework, whose initial estimation is based on -means algorithm in view of prior information on MRI. Thereby, this algorithm is also more reliable and effective.
Recently, the level set method has been applied to simultaneously segment images while estimating the bias field [8][9][10][11][12][13][14][15]. For example, Li et al. [8] introduced a local weighted -means clustering-based variational level set approach to estimate the bias field and segment the images with intensity nonuniformity. A unique feature of this method is that the calculated bias field is essentially guaranteed by the data term in the variational formulation, without any extra effort to maintain the smoothness of the bias field. In [9], Zhang et al. presented a statistical and variational multiphase level set (SVMLS). This approach used the Gaussian distribution with spatially varying mean and variance to describe the image model. Thereby, it can distinguish regions with similar intensity means but different variances. In [10], Wang and Pan proposed a new image-guided regularization to restrict the level set function. In this method, tissue segmentation and bias field estimation are unified into a single Bayesian inference framework and are simultaneously achieved by minimizing the objective energy functional. So the method can be used for accurate segmentation and bias correction of medical images in the presence of severe intensity inhomogeneity. However, all of them are sensitive to initialization of the contour to some extent.
In this study, we propose a novel LGDF model in terms of noise and robustness for simultaneous image segmentation and bias field estimation. Firstly, according to the observed signal model of the image with intensity inhomogeneity, the LGDF energy function that includes the local entropy is defined for driving the evolution contour of the level set toward the desired boundary. The means of this objective function have a multiplicative factor that estimates the bias field in the transformed domain. Furthermore, by incorporating the bias field prior into a variational level set formulation with a regularization term, our method can be used for segmentation and bias field correction. Finally, we have also made some comparisons with several state-of-the-art models to show the superiority of our method over the traditional local region-based methods.
The contributions of this paper may be summarized as follows: (1) Our model is built by simultaneous segmentation and bias field correction within a single framework. Hence, it can make full use of a priori knowledge on the bias field and level set function. (2) By incorporating the local entropy information, our model can estimate the bias field more accurately.

Li's Method.
In order to overcome the problem of intensity inhomogeneity, Li et al. [8] introduced a variational level set approach to estimate the bias field and segment the images with intensity nonuniformity. The method is based on the following model to describe an observed image: where is the measure image, is the bias field, is the true image, and is the additive noise. The assumption on the true image and the bias field has the following properties: (i) the bias field is slowly varying in the image domain and (ii) the true image intensities are approximately a constant within each class of tissue; that is, ( ) ≈ for ∈ Ω , with {Ω } =1 being a partition of Ω. According to the above assumption, this method applies a circular neighborhood with a small radius centered at each point in the image domain Ω, defined by = { : | − | ≤ }. Then, the value ( ) for all in the circular neighborhood can be well approximated by ( ). Therefore, the intensities ( ) ( ) in each subregion ∩ Ω are approximately the constant ( ) . Considering that the intensities ( ) in the neighborhood can be classified into classes, the local clustering criterion can be described as where {Ω } =1 denotes a partition of the image domain Ω, ( ) is the cluster center to be optimized, and ( − ) is a nonnegative weighting function.
where |Ω| is the number of pixels in ∩ Ω . By using maximum a posteriori probability and Bayes' rule [9], the local Gaussian distribution fitting energy can be described as follows: (4) where are the positive constants and , ( ( )) is the probability density in region ∩ Ω , which is defined as

Local Entropy.
The definition of entropy has first been introduced by Shannon [16] and has been further developed by the information theory community. As far as segmentation is concerned, a region may be characterized using the average amount of information, namely, the entropy, carried by the intensity of the region, or using joint entropy for features combination. The entropy of an image is expressed as follows [17]: where is the probability of the given images . In this study, we give the definition of the local entropy in a spatially continuous domain Ω ⊂ Ω; then the local entropy of the point can be written as where ( , Ω ) is the grey level distribution. It is given by We first apply (7) and (8) to compute the local entropy on images, which are displayed in Figure 1. As discussed in our previous work [18], the local entropy has good robustness for noise and intensity inhomogeneity.

New
LGDF Energy Based on Local Entropy. As mentioned above, the energy in (4) is sensitive to initialization and noise. In order to handle these problems, we use the local entropy ( , Ω ) defined in (7) to describe the intensity variation in a neighbourhood Ω of a point . Therefore, we redefined the new LGDF energy as follows: where ( ) = ( , ( , )) is the local entropy of , The ultimate goal is to minimize NLGDF for all the center points in the image domain Ω, which directs us to define the following energy: Substituting (5) into (10), we can obtain the proposed energy formulation as follows:

Level Set Formulation.
In this section, level set formulation is used to solve our energy NLGDF in (11). Let : Ω → be a level set function; then its signs partition the domain into two disjoint regions Ω 1 = { : ( ) > 0} and Ω 2 = { : ( ) < 0}. We assume that the image domain Ω can be separated into two regions Ω 1 and Ω 2 . These two regions can be represented by their membership functions defined as 1 ( ( )) = ( ( )) and 2 ( ( )) = 1− ( ( )), respectively, where (•) is the Heaviside function. Thus, the energy in (11) can be written as the following level set formulation: In our implementation, Heaviside function is usually approximated by a smoothing function defined by The derivative of is the following smoothing function: To derive an accurate and smooth contour, we need to add a length term ( ) and a regularization term ( ). Therefore, the entire energy functional is where ( ) = ∫ ∇ ( ( )) , and ] and are positive constants. By minimization of the energy ( , , , ) in (15), image segmentation and bias field estimation can be simultaneously achieved. The procedure is as follows: in each iteration, we minimize the energy ( , , , ) with respect to each of its variables , , , and , given the other three updated in previous iteration. The energy minimization with respect to each variable can be obtained as follows.
For fixed , , and , minimization of the energy in (15) with respect to can be obtained by solving the following gradient flow equation: where / is the Gâteaux derivative of the functional . The Gâteaux derivative can be obtained by using calculus of variation [19]; hence the gradient flow equation is expressed by where 1 and 2 are the functions as follows: For fixed , , and , the optimal that minimizes the energy is given by By fixing the other variables in (15), we obtain the minimizer of as follows: By fixing the other variables in (15), we find an optimal as follows: Remark. The minimization problem in (18) is nonconvex, so we need to make the convergence criteria. The terminal condition is similar to | ( ) − ( +1) | < 0.001 or the number of iterations that is set in advance. If the convergence criteria are reached, stop the iteration.

Description of Algorithm
Steps. For a deep understanding of our model, the iterative procedure is summarized below in this section.
Step 6. Update the level set function according to (18).
Step 7. Check the convergence criteria and iteration number. If the iteration number reached a predetermined maximum number or | ( ) − ( +1) | < 0.001 ( ( ) is the th iteration result of ), stop the iteration; otherwise, return to Step 5.

Experimental Results
In this subsection, Li's model [8], Zhang's model [9], Wang's model [10], RSF model [20], LGDF model [21], and the method of this paper are applied on a variety of synthetic images and medical images. The algorithm is implemented

Segmentation of Synthetic and Real Images.
We firstly apply our method to segment five synthetic images, which are displayed in Figure 2(a). These images are corrupted by strong noise and intensity inhomogeneity. The final segmentation results obtained after the convergence of our algorithm are displayed in Figure 2(b). The computed bias field and the bias-corrected images are shown in Figures 2(c) and 2(d).
As the local regional difference is considered, incorrect estimations of the true image in local region can be corrected in every iteration step. It can be seen from Figure 2 that the new contours gradually emerge during the evolution process. In the final segmentation results, the complete object boundaries can be effectively extracted despite the impact of intensity inhomogeneity or heavy noise. We also test our method on real images as shown in Figure 3. It reveals that the proposed model can segment multiple objectives successfully in both real images via driving the contours to desirable boundaries.

Segmentation of Medical Images.
We also evaluate the performance of our model on five medical images with obvious intensity inhomogeneity and high noise. correctly. Taking more statistical characteristics into account, LGDF method and Zhang's method yield similar visual quality to our model in some images. In order to test the performance of our model, the computational time and iterations for segmentation are presented in Table 1. Compared to LGDF method and Zhang's method, the proposed model is much easier to converge. The reason is that, with the local entropy, only a simple alternating optimization is needed in every iteration step. Experiments have proven that our method has higher computing efficiency besides the accurate segmentation.

Segmentation of Noise Images.
In order to evaluate the sensitivity to noise, we apply our method to the images corrupted by various levels of Gaussian noise, as shown in Figure 5. Figure 5(a) shows the original images with the ground truth. Figure 5(b) shows the images with Gaussian noise levels {0.05, 0.1, 0.15, 0.2, 0.25, 0.3}, respectively. We observe from Figure 5 that RSF, Li's, LGDF, and Wang's models successfully extract the object when images are corrupted by noise of lower strength, while the other two models can segment both objects from all noisy images successfully. To further illustrate the effectiveness of our method, we utilize the dice similarity coefficient (DSC) [22][23][24] to evaluate the performances of both models quantitatively. If 1 and 2 stand for the areas enclosed by contours obtained by the model and the manual method, respectively, then the DSC metric is defined as follows: where (•) indicates the number of pixels in the enclosed region. The DSC indices by applying the three compared methods are reported in Figure 6, from where it is seen that the DSC value of our model is larger than those of other methods. This demonstrates the superior performance of our method in segmentation of the images with high noise. By quantitative comparison, we can see the proposed model has good robustness to Gaussian noise. Figure 7 shows five stars images with different degree of intensity inhomogeneity. Figure 7(a) represents original image with initial contours, whereas the segmentation results obtained by RSF method, Li's method,

Quantitative Evaluation.
LGDF method, Zhang's method, Wang's method, and our method are shown from Figure 7(b) to Figure 7(g), respectively. Visual inspection clearly shows that all algorithms can segment the images precisely when intensity inhomogeneity is not strong, as shown in the first and second columns. With the increasing of intensity inhomogeneity, segmentation results of Li's method, LGDF method, and Wang's method show that they are not able to strictly find the object boundary. To further measure the quality of the extracted objects, Jaccard similarity (JS) [25,26] is used as a quantitative measure to evaluate the segmentation results of six methods. The JS index is the ratio between two regions 1 and 2 , which can be defined as JS = | 1 ∩ 2 |/| 1 ∪ 2 |. The corresponding JS values for Figure 7 are shown in Figure 8. It can be seen that our method is superior in terms of accuracy than the other models even if strong intensity inhomogeneity exists. This means that our method is very robust to image intensity inhomogeneity.
As shown in Figure 9, we carry out the experiments with Gaussian noise levels ranging from 0.1% to 0.6%. The JS values of these algorithms are listed in Table 2 and the best results of these algorithms are in bold. From this table, we can see that the proposed method obtained the highest JS values in spite of high level noise. This analysis indicates that our method also has good robustness toward the levels of noise.
To further measure the quality of the extracted objects, we test these methods on a strong intensity inhomogeneity image with different initial contours. The results of the compared   methods are reported in Figure 10. The corresponding JS values for Figure 10 are shown in Figure 11. It can be seen that our method is superior in terms of accuracy than the other models even if strong intensity inhomogeneity exists. Moreover, the JS values have few differences for different initial contours. This means that our method is very robust to initial contours.

Analysis for the Parameters.
In this section, we simply discuss the parameters that need to be manually given to obtain appropriate segmentation results. Generally, time step Δ , penalty term coefficient , binary value 0 , weighting coefficients 1 and 2 , and the positive constant of Heaviside function are relatively stable for all the experiment images. However, the parameters and ] seem to be sensitive. Therefore, it is necessary to discuss the relationship between the segmentation results and these parameters, fixing the other parameters and only changing one parameter each time, using the images in Figure 5. From the segmentation results shown in Figure 12, it is illustrated that the scale parameter is the standard deviation of Gaussian kernel. Increasing the value of will introduce more local image information.  Hence, higher may lead to oversmoothed segmentation of the images with abundant details and textures. The regularity term coefficient ] can be adjusted to smooth the curves in a way that the smoothness of the curve increases when ] increases. On the contrary, when ] is too small, the results may be smooth enough and the obtained contour is sensitive to noise.

Conclusion
In this paper, we propose a novel LGDF model based on local entropy to simultaneously correct the bias field and segment the images. The local Gaussian distribution fitting term is responsible for attracting the contour toward object boundaries. By including the local entropy, our method can handle noise and intensity inhomogeneity efficiently. The experimental results on synthetic and medical images show the superiority of our method over several state-of-the-art active contour models. However, our model cannot segment images with different tissue types, such as brain MRI or tumor PET images. In the future, our model will be extended from two-phase to multiphase level set formulation, which would further enhance its capability in processing more complex medical images.