Abstract
This paper presents a variational level set method for simultaneous segmentation and bias field estimation of medical images with intensity inhomogeneity. In our model, the statistics of image intensities belonging to each different tissue in local regions are characterized by Gaussian distributions with different means and variances. According to maximum a posteriori probability (MAP) and Bayes’ rule, we first derive a local objective function for image intensities in a neighborhood around each pixel. Then this local objective function is integrated with respect to the neighborhood center over the entire image domain to give a global criterion. In level set framework, this global criterion defines an energy in terms of the level set functions that represent a partition of the image domain and a bias field that accounts for the intensity inhomogeneity of the image. Therefore, image segmentation and bias field estimation are simultaneously achieved via a level set evolution process. Experimental results for synthetic and real images show desirable performances of our method.
1. Introduction
Image segmentation is an important and necessary step in various image processing and computer vision applications. However, due to the imperfection of the image acquisition process, intensity inhomogeneity (or bias field) is often seen in many real-world images, especially in medical images [1]. For example, the intensity inhomogeneity in magnetic resonance (MR) images usually manifests itself as a smooth intensity variation across the image [2]. Thus the resultant intensities of the same tissue vary with the locations of the tissue within the image. This can cause serious misclassifications when intensity-based segmentation algorithms are used. Therefore, intensity inhomogeneity has been challenging difficulty in image segmentation.
The level set method, originally used as numerical technique for tracking interfaces and shapes [3], has been increasingly applied to image segmentation in the past decades [4, 5]. Compared with the classical image segmentation methods such as edge detection, thresholding, and region growing, level set methods have three desirable advantages. First, they can achieve subpixel accuracy of object boundaries [6]. Second, they allow incorporation of various prior knowledge, for example, shape and intensity distribution, so as to get more robust segmentation [7, 8]. Third, they can provide smooth and closed contours as segmentation results, which are necessary and can be readily used for further applications such as shape analysis and recognition [9]. In general, the existing level set methods can be categorized into two classes: edge-based models [6, 10–12] and region-based models [9, 13–17].
Edge-based models typically use image gradient as an image-based force to attract the contour toward object boundaries. These models have been successfully used for general images with strong object boundaries, but they are generally sensitive to the initial conditions and may suffer from boundary leakage problem for medical images which typically contain weak boundaries. These drawbacks greatly limit their utilities for medical images. Region-based models use a certain region descriptor to guide the motion of the active contour. Therefore they are less sensitive to initial contours and have better performance for images with weak object boundaries. A typical example is the piecewise constant (PC) model proposed in [13]. This model assumes that image intensities are statistically homogeneous in each region and thus always fails to segment images with intensity inhomogeneity. Intensity inhomogeneity can be dealt with by more complicated models than PC model. Tsai et al. [16] and Vese and Chan [17] independently proposed two similar region-based models, widely known as piecewise smooth (PS) models, for segmentation of more general images. The PS models cast image segmentation as a problem of finding an optimal approximation of the original image by a piecewise smooth function. Although the PS models do not assume homogeneity of image intensities and therefore are able to segment images with intensity inhomogeneity to some extent, they are computationally too expensive.
Recently, local intensity information has been incorporated into level set methods to effectively handle intensity inhomogeneity [1, 9, 18–23]. For example, Li et al. [9] defined a region-scalable fitting (RSF) energy in terms of a contour and two fitting functions that locally approximate the image intensities on the two sides of the contours. Based on the multiplicative model of images with intensity inhomogeneity and a derived local intensity clustering property, Li et al. [21] presented a variational level set framework for simultaneous segmentation and bias correction. Similarly, Chen et al. [19] adopted localized -means-type clustering to define an energy which contains the bias field as a variable, and thus the energy minimization can also implement image segmentation and bias field estimation simultaneously. These models essentially draw upon local intensity means, which enable them to cope with intensity inhomogeneity. However, the local intensity means do not provide enough information for accurate segmentation, especially in the presence of strong noise and intensity inhomogeneity [1]. Therefore, more complete statistical characteristics of local intensities recently have to be taken into account. For instance, Rosenhahn et al. [18] used both local intensity means and variances to characterize the local intensity distribution in their proposed level set method. However, the local intensity means and variances are defined empirically in their models. Wang et al. [1, 20] proposed a local Gaussian distribution fitting (LGDF) energy with a level set function and local means and variances as variables, in which the local intensity means and variances are strictly derived from a variational principle, instead of being defined empirically. However, LGDF model cannot estimate the bias field so as to correct the intensity inhomogeneity in the original image. An improved LGDF model has been proposed by Chen et al. [22] to implement image segmentation and bias field correction simultaneously. Although these two LGDF models were derived based on maximum a posteriori probability (MAP) rule, they assumed the a priori probability of each partition among all partitions is equiprobable; that in other words, the MAP rule actually reduces to maximum likelihood (ML) rule. In [23], based on MAP rule, Shahvaran et al. proposed a variational level set combined with Markov random field model for simultaneous segmentation and bias filed correction of MR images, in which the a priori probability was approximated by the normalized pseudolikelihood [24].
In this paper, we propose a new level set method based on local region statistics and MAP rule for simultaneous image segmentation and bias field estimation. Based on the additive model of images with intensity inhomogeneity, we describe image intensities belonging to each different tissue in local regions by Gaussian distributions with different means and variances. By using maximum a posteriori probability (MAP) and Bayes’ rule, we first define a local energy for image intensities in a neighborhood around each pixel. The local energy is then integrated over the entire image domain to form a global energy, which is converted into a level set formulation. Minimization of this global energy is achieved by an interleaved process of level set evolution and estimation of the bias field. As an important application, our method can be used for accurate segmentation and bias correction of medical images in the presence of severe intensity inhomogeneity and noise.
The remainder of this paper is organized as follows. Section 2 reviews three different models of image with intensity inhomogeneity and some relevant level set methods. Section 3 gives the energy formulation of our method and its minimization in level set framework is presented in Section 4. In Section 5, experimental results are presented and analyzed using a set of synthetic and real images, followed by some discussions in Section 6. This paper is summarized in Section 7.
2. Background
2.1. Models of Image with Intensity Inhomogeneity
According to [2, 25], three commonly used models of image with intensity inhomogeneity have been proposed in the literature, depending on how the true image , bias field , and noise interact. Let be the observed image; the first model assumes that the noise only arises from the scanner and is therefore independent of the bias field, which is defined as
In the second model, only biological noise is considered, which is scaled by the bias field , so that the signal to noise ratio (SNR) is preserved
The third model is based on log-transformed intensities, by which the multiplicative bias field becomes additive. However, if two main sources of noise, namely, the biological noise and the scanner noise, should be considered in log-space, the noise formation will be nontrivial and unknown. Therefore, for methodological convenience, the noise is still assumed to be additive Gaussian [2]. The image model is expressed as
2.2. Li’s Method
A generally accepted assumption on the bias field is that it is smooth or slowly varying. Ideally, the intensity belonging to the th tissue should take a specific value , which represents the measured physical property [21]. Based on these two properties and the image model (1), Li et al. [21] firstly defined a local objective function to classify the data in the circular neighborhood centered on into clusters using a -means clustering method, and then this local objective function is integrated with respect to the neighborhood center over the entire image domain O to formulate the proposed energy function as follows. where denote a partition of the image domain O, is the th cluster center within the neighborhood , is a truncated Gaussian kernel which controls the size of the neighborhood . Note that the term inside the outer brackets is the local objective function in the neighborhood .
2.3. Chen’s Method
Similar to Li’s method, Chen et al. [19] applied the -means clustering to classify the log-transformed intensities and thus proposed the following energy function: where , represent and , respectively. is the log-transformed value of the intensity belonging to the th tissue. Correspondingly, is the th cluster center within the neighborhood .
In essence, both Li’s and Chen’s methods utilize the local intensity means, and hence they can handle intensity inhomogeneity. However, the local intensity means do not provide enough information for accurate segmentation, especially in the presence of strong noise and intensity inhomogeneity. Besides, these two methods ignore the noise, which is also a key factor to influence the segmentation accuracy.
2.4. Chen’s Improved Method
For the image model (2), taking the logarithmic transformation of both sides, we have Let , , and denote , , and , respectively. Equation (6) can be written as Assuming that the intensities of different tissue regions within the neighborhood have different Gaussian distributions, that is, means and variances, Chen et al. [22] utilized MAP rule (strictly speaking, ML rule) to derive the local Gaussian distribution fitting energy as follows: where is the probability density of the th region within the neighborhood , which is defined as where , are intensity mean and standard deviation of the th region within the neighborhood , respectively.
3. Our Proposed Energy Formulation
In this paper, we adopt the image model (3), which was used in [26–28]. Let , , and denote , , and , respectively, (3) can be rewritten as where noise is assumed to be zero-mean Gaussian noise with variance within the th tissue domain.
To effectively exploit information on local intensities, we need to characterize the distribution of local intensities via partition of neighborhood as in [1]. For each point in the image domain O, we consider a circular neighborhood with a small radius , which is defined as . The partition ( is the total number of segmented regions) of the entire domain O induces a partition of the neighborhood , that is, forms a partition of [1, 19]. For a slowly varying bias field , the values for all in the circular neighborhood can be well approximated by its value at the center of , that is
If we assume that the intensity belonging to the th tissue should take a specific value as in [19], the intensities within subregion can be described as a Gaussian distribution with mean () and variance , according to (10). In other words, the probability density in subregion is
Now we consider the segmentation of the circular neighborhood based on maximum a posteriori probability (MAP) rule. Let be the a posteriori probability of the subregion given the neighborhood gray values . According to Bayes’ rule: where is the probability density in region , which is denoted by in (12). is the a priori probability of the partition among all partitions of , which can be denoted by . The a priori probability is independent of the choice of the region. Therefore, (13) can be simplified as
Assuming that the pixels within each region are independent, the MAP can be achieved by finding the maximum of . Taking a logarithm, the maximization can be converted into the minimization of the following energy:
Generally, pixels far away from the neighborhood center are expected to have less influence than pixels close to . Therefore, we incorporate a nonnegative and monotonically decreasing weighting function into the energy (15) to constrain the influence of different pixels. As in [1, 9, 19–23], we choose a truncated Gaussian kernel as the weighting function ; that is, where is the standard deviation of the Gaussian kernel and is a constant to normalize the Gaussian kernel. With this weighting function, the above objective function can be rewritten as as for .
The ultimate goal is to minimize for all the center points in the image domain , which directs us to define the following double integral energy:
Substituting (12) into (18), we can obtain the proposed energy formulation as follows:
Note that our method is different from the method proposed in [23]. First, the image model of intensity inhomogeneity we use is based on log-transformed intensities. Second, we incorporate the a priori probability as a variable into the energy function; thus the optimal value will be found by the variational principle.
4. Energy Minimization in Level Set Framework
The proposed energy in (19) is expressed in terms of the regions . It is difficult to derive a solution to the energy minimization problem from this expression of . Alternatively, we can use one or multiple level set functions to represent the disjoint regions . Thus this energy can be converted into an equivalent level set formulation, which can be solved by using well-established variational methods [29]. Besides, the energy is subject to the constraint ; therefore, its minimization can be derived using the Lagrange multiplier method.
4.1. Two-Phase Level Set Formulation
We assume that the image domain can be partitioned into two regions corresponding to the object and background; that is, . These two regions can be represented as the regions outside and inside the zero level set of function ; that is, and . Using the Heaviside function , the energy in (19) can be expressed as where and . Heaviside function is usually approximated by a smoothing function defined by with as in [13]. The derivative of is the following smoothing function:
As in typical level set methods, we need to regularize the zero level set by penalizing its length to obtain a smooth contour during evolution. This can be realized by the length term
To avoid the expensive reinitialization procedures, we regularize the level set function by penalizing its deviation from a signed distance function as in [30], characterized by the following energy:
Therefore, the entire energy functional is defined as where , are positive constants.
The minimization of energy can be achieved by an iterative process: in each iteration, we minimize the energy with respect to each of its variables , , , , and , given the other four updated in previous iteration. We give the solution to the energy minimization with respect to each variable as follows.
For fixed , , , and , the minimization of in (25) with respect to can be achieved by using standard gradient descent method, namely, solving the gradient flow equation where is the Gâteaux derivative of the energy .
By calculus of variations, we can compute the Gâteaux derivative and express the corresponding gradient flow equation as where is the gradient operator, is the divergence operator, is the Laplace operator, and
For fixed , , , and , the optimal that minimizes the energy is given by
Similarly, we can obtain
4.2. Multiphase Level Set Formulation
For the case of , we can use multiple level set functions to represent multiple regions with as in [17]. For notational simplicity, we denote , , , and . The energy for general multiphase formulation of our method can be defined as where are membership functions of the regions , , such that
For example, in the case of , we use two level set functions and to define , , and to give a three-phase level set formulation of our method. For , can be defined as , , , and .
The minimization of the energy in (33) with respect to the variable can be performed by solving the following gradient flow equations
The energy minimization with respect to , , , and can be achieved by the same procedure as in the two-phase case. And it is easy to show that optimal , , , and are given by (29), (30), (31), and (32), respectively, only by replacing with .
The implementation procedure of our proposed method can be summarized as follows:
Step 1. Initialize the level set function (or ), the bias field , and the a priori probability .
Step 2. Update and using (29) and (30), respectively.
Step 3. Update the level set functions (or ) according to (27) (or (35)).
Step 4. Update the bias field using (31).
Step 5. Update the a priori probability using (32).
Step 6. Repeat Steps 2–5 until the convergence criteria are met.
5. Experimental Results
In this section, the proposed method has been tested on both synthetic and real images from different modalities. The level set function can be simply initialized as a binary step function as in [9], which takes a negative constant value (e.g., ) inside a region and a positive constant value outside it. The bias field is initially equal to zero for all . The a priori probability is initially equiprobable. For the experiments in this paper, we set the iteration time step , standard deviation of the Gaussian kernel , neighborhood radius of the Gaussian kernel , weighting coefficients and for all images on a trial basis.
5.1. ?Application to 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. The initial contours superposed on the original images are chosen manually. The intermediate contours obtained by running our method are shown in the second and third columns, and the final contours obtained after the convergence of our algorithm are shown in the fourth column. The two images in the fifth column are the estimated bias fields. It is revealed from Figure 1 that the new contours gradually emerge during the evolution process. In the final segmentation results, the complete boundaries of multiple objects can be successfully extracted despite of the impact of noise and intensity inhomogeneity.
(a)
(b)
(c)
(d)
(e)
5.2. ?Application to Real Medical Images
We also evaluate the performance of the proposed method on a set of real medical images. The first two rows of Figure 2 show two 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. A coronal slice of brain MR image is shown in the third row of Figure 2, in which intensity inhomogeneity is quite obvious. The original images and initial contours, the intermediate results, the final results, and the estimated bias fields are shown from the left column to the right column of Figure 2. These satisfactory segmentation results demonstrate desirable performance of our method for these challenging medical images.
(a)
(b)
(c)
(d)
(e)
5.3. ?MR Image Segmentation and Bias Correction
We then focus on the application of the proposed method, which is implemented using three-phase model, to segmentation and bias correction of brain MR images. The first column of Figure 3 shows two 3T MR images with obvious intensity inhomogeneity. Some parts of white matter (WM) and gray matter (GM) boundaries are quite fuzzy due to low contrast. In addition, the noise also causes certain difficulties. However, from the final contours and the segmentation results, respectively, shown in the second and third columns of Figure 3, it is obvious that our method successfully extracts the WM and GM in spite of the problems. The estimated bias fields and the bias corrected images are shown in the fourth and fifth columns. It can be seen that the intensities within each tissue become quite homogeneous in the bias corrected images. The improvement of the image quality in term of intensity homogeneity can be also demonstrated by comparing the histograms of the original images and the bias corrected images in Figure 4. Obviously, there are three well-defined and well-separated peaks in the histograms of the bias corrected images, each corresponding to a tissue or the background in the images. In contrast, the histograms of the original images do not have such well-separated peaks due to the mixture of the intensity distribution caused by the bias field.
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
5.4. ?Comparison with Other Methods
In this subsection, we first compare the proposed method with Chen’s method, Li’s method, and Chen’s improved method.
Figure 5 shows the results of four algorithms for two synthetic images with intensity inhomogeneity. The first image contains a square object and background of the same intensity means but different variances. Since only the local intensity means are taken into account, Chen’s and Li’s methods fail to locate the object boundaries. By contrast, Chen’s improved method and our method, which utilize local intensity means and variances efficiently, can differentiate the object and background using local intensity variance information. The segmentation task for the second image is to extract three objects with different shape. Due to the complicated intensity inhomogeneity and strong noise, Chen’s and Li’s models fall into a local minimum and thus obtain unsatisfactory segmentation results. Nevertheless, both Chen’s improved method and our method successfully implement the segmentation task.
(a)
(b)
(c)
(d)
(e)
(f)
Figure 6 presents the results for two ultrasound images. The original images have severe intensity inhomogeneity and strong noise. Chen’s and Li’s methods are unable to overcome these handicaps, and hence the segmentation results have many false contours. Because both Chen’s improved method and our method fully exploit intensity distribution information in local regions, their segmentation results successfully separate the lesion from normal tissues.
(a)
(b)
(c)
(d)
(e)
(f)
?Next, in order to quantitatively compare the proposed method with the above-mentioned other three methods, we use 20 T1-weighted simulated brain MR images with ground truth from Brain-Web in the link http://www.bic.mni.mcgill.ca/brainweb/. All the MR images are selected from the same volume data with 40% image intensity inhomogeneity and 7% noise. The task of segmentation is to partition the brain MR images into four regions, that is, white matter (WM), gray matter (GM), cerebral spinal fluid (CSF), and background. The segmentation results obtained by applying these methods to two sample images are shown in Figure 7. It can be observed that in contrast with Chen’s improved method and our method, both Chen’s and Li’s methods misclassify more WM into GM, especially in the lower part of the image in the second row of Figure 7. We adopt the Jaccard similarity () [2] as a measurement of the segmentation accuracy. The between two regions and is defined as . We compute the between the segmented region obtained by the algorithm and the corresponding region given by the ground truth. The closer the value to 1, the better the segmentation result. The resulting average values over 20 images for these four methods are listed in Table 1. It demonstrates that both Chen’s improved method and our method, which use local means and variances, achieve higher accuracy than the other two methods which only utilize local means. However, Chen’s improved method virtually derived from ML rule ignores the a priori information with respect to each partition and therefore has lower accuracy than our method which is really based on MAP rule.
(a)
(b)
(c)
(d)
(e)
(f)
Finally, we also quantitatively evaluate the performance of our method for 20 real MR images selected from the Internet Brain Segmentation Repository (IBSR). The data sets are available at http://www.nitrc.org/projects/ibsr. The test images are selected from sequence 100_23. The expert manual segmentation results provide the basis for a “ground truth” set to be used for comparative study. Figure 8 shows the segmentation results of our method in two sample images and corresponding ground truth. We compare the accuracy of our method with that of six segmentation methods provided by the IBSR, including the MAP, adaptive MAP, bias MAP, FCM, maximum likelihood, and tree-structure -means algorithm. Furthermore, two additional methods, FSL-FAST [31] and SPM [25], which report average results on IBSR data set, are also included. The average values of nine methods are shown in Table 2. It can be seen that our method outperforms other eight methods in terms of segmentation accuracy.
(a)
(b)
(c)
6. Discussion
Our proposed method is based on log-transformed data, and hence low intensity pixels have to be excluded to avoid numerical problems [25]. As in other level set methods, there are parameters which need to be adjusted to obtain appropriate segmentation results. The parameter in (25) or (33) influences the regularization term of the level set, and its effect has been detailed in [30]. As shown in [22], the parameter can be adjusted to smooth the curves, in a way that the smoothness of the curve increases when increases. The truncated Gaussian kernel is usually constructed as a mask (or for 3D data), where is a small odd number typically within . In our experiments, we find that can appropriately complete segmentation task for all the images.
In addition, it can be observed from (28)–(32) that for every iteration, a few convolution operations are performed. The convolution operations require a high computational complexity , where and are the size of the image and the Gaussian kernel, respectively. For the images with the size , the average CPU time using Matlab codes on a PC with Intel Dual-Core 2.93?GHz Processor and 2G RAM is about 20?s. For the 3D data with the size from the IBSR, the average CPU time is about 20 minutes. Obviously, our method must be computationally expensive to segment the whole brain. However, with the initial curve near the edge of the objects, the proposed method may need less iterations and less CPU time.
7. Conclusion
In this paper, we have proposed a new level set method for segmentation and bias field estimation of images with intensity inhomogeneity. Based on local region statistics and maximum a posteriori probability, we define an energy in terms of the level set functions that represent a partition of the image domain and a bias field that accounts for the intensity inhomogeneity of the image. Our method can estimate intensity inhomogeneity, handle noise efficiently, and also allow flexible initialization. In addition, the regularity of the level set function is intrinsically preserved by the level set regularization term to ensure an accurate computation avoiding expensive reinitialization procedures. Comparisons with other approaches demonstrate the superior performance of the proposed method. However, our proposed method does have some limitations such as high computation cost. In our future work, we plan to reduce the computational scheme by using the Split Bregman method [32], instead of gradient descent method used in this paper.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to thank the anonymous reviewers for their helpful comments. This work was supported by the National Natural Science Foundation of China (60903127, 61202314), the NPU Foundation for Fundamental Research (JCY20130130), the “Top Paper” Project (13GH013308) and the Ao-Xiang Star Project (11GH0315) at the Northwestern Polytechnical University, Xi’an, Shaanxi, China.