Active Contour Driven by Local Region Statistics and Maximum A Posteriori Probability for Medical Image Segmentation

This paper presents a novel active contour model in a variational level set formulation for simultaneous segmentation and bias field estimation of medical images. An energy function is formulated based on improved Kullback-Leibler distance (KLD) with likelihood ratio. According to the additive model of images with intensity inhomogeneity, we characterize the statistics of image intensities belonging to each different object in local regions as Gaussian distributions with different means and variances. Then, we use the Gaussian distribution with bias field as a local region descriptor in level set formulation for segmentation and bias field correction of the images with inhomogeneous intensities. Therefore, image segmentation and bias field estimation are simultaneously achieved by minimizing the level set formulation. Experimental results demonstrate desirable performance of the proposed method for different medical images with weak boundaries and noise.


Introduction
Image segmentation is an important procedure in various image processing and computer vision applications.Many algorithms such as histogram thresholding, region growing, fuzzy c-means (FCM), and active contour models (ACM) have been proposed for image segmentation.However, due to the presence of noise, low contrast, and intensity inhomogeneity in medical images, the segmentation is still a challenging problem in the majority of applications.
Since being introduced by Kass et al. [1], active contour ideas have been widely used in the computer vision field.There are mainly two kinds of active contour models: edgebased methods [2][3][4] and region-based methods [5][6][7][8].Edgebased methods utilize image gradients to drive the curve evolution process, which is usually sensitive to image noise and weak edges.On the other hand, instead of employing gradient information, region-based methods utilize image region information to capture targets of interest.One of the early efforts towards region-based active contours was made by the M-S model [9], which approximates a given image by a piecewise smooth image.Compared to the edge-based model, the region-based models are less sensitive to initial contours and have better performance for images with weak object boundaries.The region-based model is insensitive to initialization of the level set function and can recognize boundaries of the object efficiently.Therefore, region-based models, especially the C-V model [5], have been widely applied in image segmentation.
Recently, Li et al. [10] proposed a variational level set approach to segment images and estimate the bias field, which is motivated by the weighted K-means clustering.The main advantages of this model are robustness to initialization, accurate segmentation of images with intensity inhomogeneity, and estimate of the bias field.However, only using intensity means as a criterion of classification, this method may not be able to precisely extract the object boundaries, especially when the object texture distributions have identical means but different variances.Wang et al. [11] proposed an active contour model driven by local Gaussian distribution fitting (LGDF) energy with a level set function and local means and variances as variables for more accurate segmentation.The LGDF model can distinguish regions with similar intensity means but different variances by using the local Gaussian distribution fitting energy.However, the drawbacks of LGDF model are that it may introduce many local minima and cannot estimate the bias field so as to correct the intensity inhomogeneity in the original image.The rest of the paper is structured as follows.In Section 2, a review on some well-known region-based models is given.Section 3 shows the proposed model in detail.In Section 4, experimental results are given.Finally, conclusions are drawn in Section 5.

Background
2.1.The Description of C-V Model.The Chan and Vese model, integrating the Mumford-Shah model and level set, utilizes the region information for segmentation [5].The level set formulation of the C-V mode, regarding the time evolution of the level set function , can be described as where () is the original image, , ],  1 ,  2 are the positive constants of each item and the   () is the Dirac function.
The first term keeps the level set function smooth, while the third and fourth terms are the internal and external forces, respectively, that drive the contour towards object boundaries.In C-V model,  1 and  2 are defined as follows: Obviously,  1 and  2 are related to the global properties of the image intensities inside and outside the contour.The C-V model has good performance in image segmentation due to its ability of obtaining a larger convergence range and being less sensitive to the initialization.But because this model does not include local information, contour curve will deviate from the target boundary and it may cause error segmentation results when the energy function reaches to the minima.

Method of Li et al. Li et al. proposed a variational level
set approach to estimate the bias field and to segment the images with inhomogeneous intensity, which is motivated by a weighted K-means clustering.It is based on the following model to describe images in [10]: where  is the measured image,  is the true image,  is the bias field,  is the noise.In [10], the assumptions on the true image  and the bias field  can be stated as follows.
(I) The bias field  is slowly varying in the entire image domain.(II) The true image intensities  are approximately a constant within each class of tissue; that is, () ≈   for  ∈ Ω  , with {Ω  }  =1 being a partition of Ω.Let   = { : ‖ − ‖ ≤ } be the circular neighborhood with a small radius  centered at each point  in the image domain Ω.Then, in each small region, an approximation is given by Considering the task of classifying the intensity () in the neighborhood   into  classes, a local intensity energy function is defined as where ()  are the cluster centers to be optimized and (− ) is a nonnegative weighting function, which is defined as follows: where  is a constant.Thus, ( 5) can be rewritten as In the method of Li et al., the case of  = 2 is mainly considered.And the image domain is divided into two regions by the zero level contour of a function ; that is,  > 0 and  < 0. Using the Heaviside function   , we can get the energy function as follows: where  1 (()) =   (()) and  2 (()) = 1 −   (()).
Since local image intensity information is embedded into the energy function, the method can deal with some kinds of images with intensities inhomogeneity and estimate of the bias field.However, only using intensity means as a criterion of classification, this method is hard to correctly segment the image when the mean values of intensities of background and object are rather close to each other.

The
LGDF Model.Wang et al. [11] proposed an active contour model driven by local Gaussian distribution fitting (LGDF) energy with a level set function and local means and variances as variables for more accurate segmentation.The energy functional for the LGDF model for each pixel using the level set function  is defined as follows: where ] ≥ 0 and  ≥ 0 are fixed parameters, ( − ) is a nonnegative weighing function, and  , (()) is where   () and   () are local intensity means and standard deviation, respectively.The LGDF model can distinguish regions with similar intensity means but different variances by using the local Gaussian distribution fitting energy.Thus the LGDF model is able to deal with both noise and intensity inhomogeneity.However, it still easily gets stuck in local minima and cannot estimate the bias field so as to correct the intensity inhomogeneity in the original image.

Proposed Method
3.1.Energy Functional Formulation.Recently, a kind of segmentation ordinations based on Kullback-Leibler distance, also named entropy, such as the minimum entropy method [12] and the maximum posterior entropy method [13], has been proposed.When comparing with other segmentation methods, these methods can achieve better segmentation accuracy, especially for image with the speckles and noises.KLD provides a mathematical formulation for measuring the difference between two probability density functions (PDFs).And the KLD-based image segmentation method significantly depends on the PDF.For discrimination between two probability density functions  1, (()) and  2, (()), the KLD can be defined as follows: We take the logarithm for  1, (()) in (11).So, the improved KLD formulations can be rewritten as In this paper, we give the definition of the KLD in a spatially continuous domain.For a fixed regular region Ω  ⊂ Ω, it is a neighborhood centered at a point  ∈ Ω.The KLD of the point  is defined by In ( 13), a kernel function is introduced.We define the following two objective functions and a local energy function: According to the Bayesian formulation and maximum a posteriori probability method [14,15], we define the following objective function for each point  in the image domain Ω: where  , (()) is given by From (5) in the above section, we can prove that   () can be approximated to ()  : where |Ω| is the number of pixels in   ∩   .
Then (16) can be rewritten as ) ,  = 1, 2. (18) Remark.The proposed model is different from the method in [11] for the following reason.We introduced the Kullback-Leibler distance with likelihood ratios as measure of classification into this model.It can enhance the accuracy of segmentation.Moreover, our model can estimate the bias field so as to correct the intensity inhomogeneity in the original image.From (15), we can see that the energy function has the following properties: (I)  = 0 denotes that two region distributions are equal; (II) if  is minimization, two region distributions are not equal completely; that is, the reparability is better.

Level Set Formulation.
Using the Heavside function   , then the energy functional (15) can written as follows: In order to obtain accurate evolution of the level set function, we need a signed distance function (SDF) to regularize the level set function.The SDF can be characterized as follows: To regularize the level set contour, we need its length to derive a smooth contour during evolution: The regularized Heaviside function   is defined by Using the following regularized Dirac function, Thus, the entire energy functional is defined as where ,  1 ,  2 , , ] are positive constants.

Gradient Descent Flow.
Minimization of the energy function () in ( 24) with respect to  can be achieved by solving the gradient descent flow equation: where By fixing , , , we find an optimal bias field that minimizes (): Using the similar optimal process, the global means   and local variance  2  can be written as Equation ( 25) is the level set evolution equation to be solved in the proposed model.The term −  ()( 1 −  2 −  3 ), which is derived from the first term of the energy (), plays a key role in the proposed model since it is responsible for driving the active contour towards object boundaries.

Experimental Results
In this section, the proposed method has been tested on both synthetic and real medical images.The results are all carried out by Matlab 2011a in the PC with Pentium CPU 2.50 and 4 GB of RAM.In our experiments, the parameters are set as follows: iteration time step Δ = 0.1, standard deviation of the Gaussian kernel  = 3, and weighting coefficients  1 = 1.05,  2 = 1.0,  = 1.0.

Application on Synthetic
Images.We first apply our method to segment two synthetic images, which are displayed in the first column of Figure 1.These two images are corrupted by strong noise and intensity inhomogeneity.And the initial contours and the final contours are plotted as red contours and blue contours, respectively.The final contours obtained by running our method are shown in the second column.The computed bias corrected images and bias fields are simultaneously obtained, as shown in the third and fourth columns.It can be seen that our method yields satisfactory segmentation results owing to considering and exploiting the image local region information, which can successfully extracts the object boundaries for these images.

Application on Real Medical
Images.We also evaluate the performance of the proposed method on a set of real medical images in Figure 2. The first row shows an MR image of a human bladder.The second row shows the results of our method on infrared images, which include some rather weak boundaries.The third row shows an X-ray vessel images with intensity inhomogeneity.In addition, parts of the vessel boundaries are quite weak.Such properties render it a nontrivial task to segment the vessels in the images.The original images and initial contours, the final results, the bias corrected images, and the bias fields are shown from the left column to the right column in Figure 2. Our method can get satisfactory segmentation results because it considers the statistical information in a transformed domain, where the intensity of the object and background are less overlapping than that in the original domain, making our method have a very strong discriminative capability for the object and background.

Application on MR Image Segmentation.
In order to further demonstrate the capability of our method in dealing with intensity inhomogeneity, we also use our method to segment images on two MR brain images with intensity inhomogeneity due to the bias field effect, which is usually seen in medical images.The segmentation results, the bias corrected images, and the bias fields are shown from the left column to the right column in Figure 3.It can be clearly seen that the image quality is significantly improved by our method.Some regions whose intensity contrast is too low to be identified are able to be distinguished clearly after the bias correction.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.The LGDF model fully exploits intensity distribution information in local regions.However, it still easily gets stuck in local minima and cannot estimate the bias field.Our method can produce much better segmentation results because it can not only distinguish regions with similar intensity means but different variances, but also estimate the bias field so as to correct the intensity inhomogeneity in Mathematical Problems in Engineering the original image, making our method have a very strong discriminative capability for the object and background.

Results for Synthetic Noisy Images.
To further demonstrate the ability of the proposed model when dealing with noisy images, we create images by adding Gaussian noise to a synthetic image in Figure 5.The two images are contaminated with multiplicative Gaussian noise with standard deviations of (0.1, 0.15, 0.2, and 0.25).It can be seen from the bottom row that the results appear similar visual quality for the original and noisy images.By quantitative comparison, we can see that our model do produce perfect results for these original and noisy images.

Conclusion
In this paper, we present a novel region-based active contour model in a variational level set formulation for simultaneous segmentation and bias field estimation of medical images.
Based on objects of images with intensity inhomogeneity, we characterize the statistics of image intensities belonging to each different object in local regions as Gaussian distributions with different means and variances.Benefiting from the local variance, our method can obtain more accurate segmentation results.Experiments on synthetic and real medical images show that our model is desirable.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 with intensity inhomogeneity.

Figure 1 :
Figure 1: Application of our method on synthetic images.Column 1: original images and initial contours.Columns 2: final contours.Column 3: bias corrected images.Column 4: estimated bias fields.

Figure 2 :Figure 3 :
Figure 2: Application of our method on real medical images.Column 1: original images and initial contours.Columns 2: final contours.Column 3: bias corrected images.Column 4: estimated bias fields.

Figure 5 :
Figure 5: Results for the images corrupted by additive Gaussian noise.

4. 4 .
Comparison with Other Methods.Figure4shows the comparison of the proposed method with C-V model[5], Li's model[10], and LGDF[11] model in real medical images with intensity inhomogeneity.The first row shows the results for ultrasound images of left ventricle with severe intensity inhomogeneity and noise problems.The second row shows the left ventricle in a tagged MR image.It is evident that the image is corrupted by noise, severe intensity inhomogeneity, and weak boundaries.The third row shows the CT images of heart.It can be observed that the C-V model assumes that the image intensity is piecewise constant and uses the global intensity means to fit the image intensity.Therefore, they do not perform well in images with intensity inhomogeneity.Li's model uses the local intensity mean to fit the measured image, and thus it can yield better segmentation results than C-V model on images with intensity inhomogeneity.However, if

2
Mathematical Problems in EngineeringIn this paper, a new region-based active contour model based on improved Kullback-Leibler distance and local neighborhood information is presented for medical image segmentation.In our model, the statistics of image intensities belonging to each different object in local regions are characterized by Gaussian distributions with different means and variances.According to maximum a posteriori probability and Bayes' rule, we first derive a local objective function for image intensities in a neighborhood around each pixel.Then this local energy function is integrated with respect to the neighborhood center over the entire image domain to give a global criterion.In a level set formulation, this global criterion defines an energy in terms of level set functions that represent a partition of the image domain and a bias field that accounts for the intensity inhomogeneity of the image.Experimental results demonstrate desirable performance of the proposed method for different medical images with weak boundaries and noise.