Abstract
Image segmentation is still an open problem especially when intensities of the objects of interest are overlapped due to the presence of intensity inhomogeneities. A bias correction embedded level set model is proposed in this paper where inhomogeneities are estimated by orthogonal primary functions. First, an inhomogeneous intensity clustering energy is defined based on global distribution characteristics of the image intensities, and membership functions of the clusters described by the level set function are then introduced to define the data term energy of the proposed model. Second, a regularization term and an arc length term are also included to regularize the level set function and smooth its zero-level set contour, respectively. Third, the proposed model is extended to multichannel and multiphase patterns to segment colorful images and images with multiple objects, respectively. Experimental results and comparison with relevant models demonstrate the advantages of the proposed model in terms of bias correction and segmentation accuracy on widely used synthetic and real images and the BrainWeb and the IBSR image repositories.
1. Introduction
Image segmentation is a fundamental but one of the most important problems in pattern recognition and computer vision [1, 2]. In general, it is aimed at separating an image into several parts corresponding to the objects of interest. Inner elements (i.e., pixels for 2D images or voxels for 3D images) of each part, recognized as components of a desired object, are considered as having an identical characteristic in terms of shape, structure, or texture [3, 4]. As well known, due to the decades of efforts of numerous researchers, a variety of impressive image segmentation methods have been proposed in the literature [5–11]. However, it is still a challenging task to extract objects of interest accurately from a complex image [12, 13]. In particular, if the image is corrupted by noises and inhomogeneous intensities (generally viewed as intensity biases), intensity homogeneity of the image will be destroyed due to intensity overlaps between different objects caused by the noises and biases, which certainly brings challenges to classical segmentation methods that are based upon edge detection or thresholding [14–16]. Unfortunately, intensity inhomogeneities exist in most of real-world images inevitably. Figure 1 gives an example to demonstrate negative effects of inhomogeneities on intensity distribution of a camera captured image and a medical brain image.
As mentioned earlier, a variety of segmentation methods have been proposed in the literature. As one class of the most popular segmentation methods, active contour models (ACMs) have been extensively studied and have proven to be specially effective for image segmentation due to their ability to elastically deform and delineate object boundaries with smooth and closed contours in subpixel accuracy [17, 18]. The fundamental idea of ACMs is to introduce a contour to represent boundaries of interested objects and then drive the contour moving toward its interior normal direction under some constraints [19]. The constraints are generally contained in a predefined energy function which will finally get its minimal value when the contour stops on true boundaries of the desired objects [20]. However, there are inherent drawbacks of traditional ACMs, i.e., initialization sensitivity and difficulties associated with topological changes in merging and splitting of the evolving contour. Therefore, since the active contour model was proposed by Kass et al. [21], many efforts have been devoted to developing improved methods to overcome the inherent drawbacks [22]. As one of the most important improvements of ACMs, level set methods regard the active contour as the zero-level set contour of a predefined one-dimension higher function named as level set function in the literature. Motion of the contour is implied in evolution of the entire level set function under a principled energy minimization framework instead of directly driving the contour itself [23]. Therefore, interesting elastic behaviours of the active contour are preserved with topological changes of the contour efficiently being handled implicitly during the evolution of the level set function. In addition, level set methods are easily extended to a higher dimension and prior knowledge of objects of interest can be incorporated into their energy framework to guide the zero-level set contour moving close to the desired boundaries [24].
Existing level set-based image segmentation methods are usually divided into two categories, which are edge-based models and region-based models, according to whether an edge indicator or a region descriptor is used to guide the motion of the zero-level set contour [25, 26]. Edge-based level set methods are good at identifying boundaries from images with strong intensity gradients. Therefore, they inevitably suffer from a weak boundary problem caused by the presence of intensity inhomogeneities and noises [27]. This drawback restricts their applications while in turn promotes the passion of researchers in this field to develop region-based level set models which take statistical information of the image intensities in general as guide descriptors to identify each region of interest [28]. In this paper, a region-based level set model is proposed where the image intensities are clustered globally and the bias correction is embedded as well. Specifically, the intensity bias, which causes inhomogeneous intensities, is estimated in the model by orthogonal primary functions based on the principle that any of a smooth function can be estimated by a linear combination of primary functions in enough order. To illustrate the problem visually, a demonstration of orthogonal Legendre functions in fitting smooth two-dimensional functions is given in Figure 2. Another contribution of this work is that the proposed model is further extended to segment multichannel images and images with multiple objects.
The paper is organized as follows. Related work and some typical and highly relevant ACM models are reviewed briefly in Section 2. Details of the proposed model are presented in Section 3. Experimental results of the proposed model on synthetic and natural images that are widely used in the literature and comparison with greatly relevant models on the BrainWeb and IBSR image repositories are given in Section 4. We analyse and discuss relationship and improvement of the proposed model with representative models, its robustness to initialization, and coefficient impact in Section 5, which is followed by conclusion and future work in Section 6.
2. Related Work
Let be a 2-dimensional continuous domain and be an image defined on the domain. Thus, the observed image can be viewed as a mapping from to . In general, the problem of image segmentation using ACMs is in fact to find an optimal contour to divide the image domain into disjoint subregions, which are generally denoted by and for binary segmentation problems in the literature.
2.1. Classical Mumford-Shah Functional Model
To find the optimal contour , Mumford and Shah propose an energy model-based image segmentation method in [29]. The basic idea of this classical model is to find a pair of the two tuples for a given image , where is a piecewise smooth approximation of . The general form of this model can be written into the following energy function, defined by where is the modulus operation and and are positive weighting coefficients. Note that unless otherwise specified, integrations are all performed on the entire image domain in this paper.
When the contour is exactly located on the desired boundaries and is piecewise smooth enough to approximate , this functional takes its minimal value and vice versa. However, it is not easy to find the optimal solution of the above-defined energy functional due to different natures of the unknown and and the nonconvexity of the functional as well.
2.2. Chan-Vese’s Piecewise Constant Model
To overcome the difficulties in solving Equation (1), Chan and Vese propose a piecewise constant case of the Mumford-Shah model in [30], which has proven to be particularly influential in binary segmentation problems. In the well-known CV model, the contour is considered as the zero-level set contour of a one-dimension higher level set function defined on the image domain , i.e., . The level set function takes positive and negative values, respectively, in the subregions and which are separated by the contour , i.e., and . Thus, membership functions and can be, respectively, used to represent these two regions by making for , for , and otherwise, both of them are . Note that is the Heaviside function. Then, the energy functional of the CV model is defined by where is the area of the subregion which is enclosed by the 0-level set contour , is the length of the 0-level set contour , and are positive weighting coefficients, and and are two constants that are used to approximate average intensities of the given image on either side of the 0-level set contour . It is obvious that and are related to the global properties of the image intensities in and , respectively. This model has also been further extended to segment images into multiple parts using multiphase level set functions [31]. But the CV model and its multiple phase extension are both on account of the assumption that intensities of the image are statistically homogeneous in each part and use different constants to estimate intensities of these parts. They are therefore well known as piecewise constant (PC) models, which will fail to segment images with intensity inhomogeneity when disordered intensity distribution introduces overlaps between objects of interest. That is to say that even though the CV model is robust to some extent with respect to noise and is also less sensitive to the initialization, it generally fails to segment images with intensity inhomogeneity [30].
2.3. The Piecewise Smooth Model
In addition to introducing a local energy term as proposed in [32] or improving original global energy by means of image local characteristics in [33], two similar ACMs are proposed by Vese and Chan [31] and Tsai et al. [34] instead under the framework of minimization of the Mumford-Shah functional to overcome the difficulty of the CV model in segmentation of images with intensity inhomogeneity. These models are widely known as piecewise smooth (PS) models where the image intensities are considered as two piecewise smooth functions instead of constants to represent intensities on either side of the contour [31] by minimizing where and are positive weighting coefficients.
Although intensity inhomogeneities are handled to some extent by the piecewise smooth models, it is obvious that the involved update of and at each iteration will certainly increase the computational burden due to solving of two partial differential equations on the entire image domain [33]. In addition, the level set function of the above model has to be periodically reinitialized using signed distance functions, which not only introduces problems like when and how it should be performed but also affects numerical accuracy in an undesirable way [35].
2.4. Region-Scalable Fitting Model
To resolve undesirable effects caused by reinitialization, a distance regularization term is proposed to regularize the level set function preserving its sign distance property in [27], defined by
In [36], besides the above-defined regularization term, the local region information is incorporated into the region-based level set method as a data term relying on the assumption that intensities are locally homogeneous. Specifically, for a given point , two fitting functions and are used to approximate image intensities in the local subneighborhoods of which come from separating its local neighborhood by the 0-level set contour of and are obviously subsets of and , respectively. Let where is a normalized even function that satisfies , if , and with the local size being implied and scalable. The level set model proposed in [36] is named as the region scalable fitting model, defined by where and are positive constants and and are positive weighting coefficients.
Wang et al. further extended the RSF model to distinguish regions with similar intensity means but different variances by introducing Gaussian distributions to describe the local image intensities [35]. This improvement is in fact based on the assumption that intensities of the image obey normal distribution. Nevertheless, the image intensities are not necessarily described by a specific distribution, i.e., the intensities vary in any positions and directions and so do the intensity inhomogeneities. Therefore, histogram of the intensities and local statistics regarding the intensity and the magnitude of gradient is used to drive the evolution of the zero-level set contour in [37, 38].
Although the above-mentioned RSF model and its improved models perform well in segmenting images with intensity inhomogeneity, the size of local scalable regions determined by the kernel function and the initial placement of the contour have a great influence on segmentation performance [39]. In fact, the final zero-level set contour of this kind of model will not be able to stop at desired locations, in particular, when all of the internal pixels of local regions only belong to either or due to an inappropriate scalable size are set. Moreover, it is obvious that the models discussed above are unable to estimate the bias field and further remove it from the inhomogeneous image to be segmented.
2.5. Local Intensity Clustering Model
In contrast, the local intensity clustering (LIC) model is proposed in the purpose of performing a satisfactory segmentation result and estimating the bias field existing in the image at the same time [40], which is based on the following two assumptions. That is, the bias field and the true image are viewed as two multiplicative components of the image and the bias field is smooth and varies slowly. That is to say, in a small enough circular neighborhood of a given point , the bias field can be seen as a constant and the standard -means clustering can be used to classify intensities in the neighborhood. Taking all the center points of the entire image into account, the energy functional of the LIC model is defined by where , and are positive weighting coefficients, and are the membership functions of and , and is a normalized function with properties described in Section 2.4.
The LIC model shows that the image is well segmented and the bias field existing in the image can be estimated at the same time [40]. However, this model is still sensitive to the size of the local scalable region and the initial location of the zero-level set contour, as the RSF model [39]. In addition, there is no explicit constraint on the bias field to guarantee its slowly varying property.
2.6. Local Inhomogeneous Intensity Clustering Model
In [41], a constraint on the bias field is therefore established to ensure its smoothly and slowly varying property by linearly combining a given set of orthogonal primary functions, denoted by , , …, and . The model proposed in [41] is named as the local inhomogeneous iNtensity clustering (LINC) model, defined by where , and are column vectors defined by and , respectively, and , ,…, are weighting coefficients of the primary functions. Note that is the transpose operator and all the other symbols represent the same with those in Section 2.5.
As demonstrated in [41], the LINC model has the capability in extracting desired objects accurately from noisy images and correcting the intensity biases simultaneously, and it is robust to initialization. Furthermore, the LINC model converges in less iterations than RSF and LIC [41]. However, convolution operation in the evolution results in a heavy computational burden.
3. The Proposed Model
In fact, the size scalable kernel function acts as a role of localization in the RSF, LIC, and LINC models. Due to the approximatively homogeneous characteristics of local image intensities, it also plays an important part to ensure smoothness of the bias field in the LIC and LINC models. That is to say that its function overlaps with the linear combination of primary functions to some extent for the LINC model except for localization. To eliminate this overlap and promote the computational performance, we propose a global inhomogeneous intensity clustering-based active model in this paper and give details of the proposed model in this section.
3.1. The Image Model and Problem Representation
Given an intensity inhomogeneous image defined on , its intensities are generally in the literature viewed as where and are, respectively, the observed and true intensities at location of the image, is the bias field accounting for the intensity inhomogeneities in the observed image, and is additive zero-mean noise. It has been pointed out that the image model defined in Equation (8) is famous in the literature. In this paper, we propose a new model to estimate the bias field rather than define a new image model. In the literature, the true image is generally regarded as constant piecewise which takes a specific intensity value for pixels belonging to the same object. We denote the number of all objects of interest by and the intensity value of the th type of objects by in this paper. Thus, the problem of image segmentation and bias correction is naturally converted into finding the optimal estimation of and the true intensity values for each type of the objects.
3.2. Representation of the Bias Field
As well known, smooth functions are usually used to approximate the bias field in the literature due to the slowly varying property of . And a function can be theoretically approximated by a linear combination of a given number of primary functions up to arbitrary accuracy, only if the number of the primary functions is sufficiently large [42]. Therefore, given a set of primary functions with a specific cardinality , denoted by , , …, and , we represent the bias field as a linear combination of these functions with , , …, and as weightings, i.e., where is the transpose operator and and are column vectors defined by and , respectively. Note that the primary functions used in this paper are orthogonal and estimation of the bias field is performed by finding the optimal coefficients .
3.3. Global Inhomogeneous Intensity Clustering (GINC)
As mentioned earlier, the true image consists of nonoverlapping regions, each corresponding to one type object of interest with a constant intensity, denoted by . Therefore, the product of image components and can be rewritten as where . Taking Equation (9) into account, the above equation is rewritten as
In consideration of the image model given in Equation (8), we have
As mentioned earlier, is additive zero-mean noise. Therefore, we define the following inhomogeneous intensity clustering energy where are positive constants to indicate the preference of the proposed model to the corresponding class. Note that when boundaries of the regions for exactly agree with the image structure, i.e., the territory distribution of all objects of interest, the above-defined energy takes its minimal value.
3.4. Two-Phase Level Set Formulation GINC2
It is obvious that the proposed energy in Equation (13) is expressed in terms of the regions , , …, and , which makes it difficult to derive a solution to minimize the energy from this expression. In the case that the image domain is separated into two disjoint regions and , i.e., , the energy defined in Equation (13) can be converted to a level set formulation by representing the two disjoint regions with a given level set function defined on . Then, the energy minimization problem can be solved by using well-established variational methods [40]. Let the level set function take negative and positive signs on either side of the 0-level set contour denoted by , which can be used to represent a partition of the domain with two disjoint regions. The disjoint regions separated by the contour can be represented by and . In consideration of properties of the Heaviside function , the regions are further represented by the following member functions and , respectively. Note that the above definition is in accordance with the one given in Section 2.2. Thus, for the case , we rewrite the energy described in Equation (13) into the following level set formulation: which is obviously a function of variables , intensity vector , and . This energy function is considered as the data term of the proposed model GINC, defined by where the representations of and are the same ones as given in Equations (2) and (4), respectively. They are used here to smooth the zero-level set contour and regularize the entire level set function .
The proposed model is essentially different from both LIC and LINC. First, there is no normalized even convolution kernel function in the proposed model and the integral is therefore one layer which is less than either LIC or LINC. Second, compared with the LIC model where there is no any specific constraint on the bias field, the linear combination of a given set of primary functions is introduced as an additional constraint to ensure its smoothness. Third, the proposed method is global based whereas the LIC and LINC models are both local based.
3.5. Extension to Multichannel Case
It is obvious that the above model defined in Equation (15) is applicable in extracting objects of interest from gray images. But multichannel images of the same scene that come from different imaging modalities or colorful images are becoming more and more common in our life. To extend the proposed model to be able to extract objects of interest from multichannel images, we first denote a given multichannel image by where is the channel number of . Let where are positive weighting coefficients that are used to control influence of the th channel and are all set to be 1 in this paper unless otherwise specified. We then rewrite Equation (14) as follows: where is an matrix defined by and is a matrix with elements defined by . We finally rewrite Equation (15) as follows:
3.6. Further Extension to Multiphase Case
Since one level set function can only be used to represent 2 subregions of image domain denoted by membership functions and , which are in fact inside and outside of the zero-level contour of , level set functions are required to represent subregions where . Thus, the subregion can be represented by the member function , i.e., for and otherwise where and . To extend the proposed model to segment multiple objects from images with intensity inhomogeneity, we first further rewrite Equation (16) as follows:
We then define and where and , respectively. Finally, we rewrite Equation (17) as follows:
3.7. Energy Minimization
In the proposed model, image segmentation and bias correction are determined by the final level set function and the optimal weighting coefficients that are obtained by minimizing the energy functional defined in Equation (19). The energy minimization is achieved by an iterative process. That is to say, the energy functional is minimized with respect to each of its variables , , and using an iteratively alternating update strategy, i.e., trying to update each of the variables in each iteration by considering the other two with known values obtained from the last iteration.
For fixed and , we minimize the energy functional with respect to using the standard gradient descent method and obtain where .
For fixed and , we minimize the energy functional with respect to by solving the equation where is a matrix with constant value and obtain where and .
For fixed and , we minimize the energy functional with respect to by solving the equation where is a matrix with constant value and obtain where and is a matrix with elements and is an -dimensional column vector, given by
3.8. Implementation
In our numerical implementation, the Heaviside function is approximated by a smooth version with , defined by The derivative of which can be deduced and written as
In this paper, the bias field is theoretically approximated by a linear combination of 15 orthogonal Legendre polynomial functions up to four order precision in our implementation, i.e., . The implementation of the proposed model is straightforwardly expressed as follows in Algorithm 1.
| ||||||||||||||
Note that the convergence criterion used in this paper is , where is the cluster center updated at the th iteration and is the Euclidean distance of ⋆.
The main additional computational cost in the proposed model is for computing in Equation (22) compared with representative models reviewed in Section 2. However, we notice that and are independent of the level set functions and clustering centers which indicate that we can compute for Equation (23) and for Equation (24) in advance and keep the results fixed during the iteration to accelerate the proposed model.
4. Results
In this section, we first evaluate effectiveness of the proposed model GINC on synthetic images that are widely used to verify ACMs and selected natural images from public datasets. We then evaluate the proposed model on two public MR brain image repositories qualitatively and quantitatively. Unless otherwise specified, we set , , , , and in this paper.
4.1. Effectiveness of GINC
In this subsection, we qualitatively evaluate the effectiveness of the proposed model on synthetic images and selected natural images from public datasets and give the validation results in the following paragraphs. Note that the synthetic and natural images are either widely used in the literature to verify active contour models or appropriate for application of the proposed model GINC.
We first apply the proposed model to three synthetic gray images (widely used to evaluate active contour models in the literature), one cardiac X-ray image, and one brain MR image in this subsection. The original images with initial 0-level set contour in green are given in the first column of Figure 3. From the red curves used to indicate final segmentation results of the proposed model on the corresponding image, also given in the first column of Figure 3, it is obvious that segmentation results of the proposed model on the images are agreed with contents contained in the images even though intensities of the images are not homogeneous due to existing of severe intensity biases as shown in Figure 3. That is to say that it is difficult to extract objects of interest from the images because intensity ranges of objects (including the background) in the images are overlapped due to severe intensity inhomogeneities existed in the images which manifests as there are no well-separated peaks in intensity histograms of the images as shown in Figure 3. However, there are well-defined and separated peaks in histograms of the bias corrected images, each corresponding to one object or the background. This demonstrates the capability of the proposed model in correcting bias fields from images with intensity inhomogeneity. Meanwhile, the biases estimated by the proposed model with orthogonal primary functions are all slowly (not sharply) varying as shown in Figure 3 which meets properties of the bias field as described in Section 1.
We then apply the proposed two-phase level set model to segment four selected natural images with three color channels from BSD database [7], namely, 135069, 42049, 3096, and 86016, respectively. The reason we selected these images is that each of the images contains only one object besides the background which can therefore be distinguished with one level set function, i.e., and . Results of the proposed model on segmentation of the images with two-phase level sets are given in Figure 4. It is obvious that the estimated biases are smoothly varying and the corrected images are more homogeneous than the originals. Furthermore, the energy functional of the proposed model defined in Equation (17) is converged (generally in less than 50 iterations) as shown in the right column of Figure 4.
We thirdly apply the proposed three-phase level set model to segment two MR brain images which are corrupted by severe intensity inhomogeneities and two selected natural images from MSRCORID database [43], namely, 164_6484 and 112_1204. The first two images are widely used to evaluate multiple-phase active contour models in the literature, and the last two images are selected because three kinds of objects are contained which are suitable for three-phase segmentation. From the results given in Figure 5, we can see that the estimated biases are smooth and the corrected images are much more homogeneous. Moreover, the extracted objects are coincided with the images.
We fourthly evaluate energy convergence of the proposed model on all above-mentioned images and show the iteration process of the proposed model on four of them in Figure 6. The images are appropriate to evaluate the proposed model in the sense of one-channel-two-phase, one-channel-multiple-phase, multiple-channel-two-phase, and multiple-channel-multiple-phase, respectively. It can be seen that the proposed model is convergent and satisfactory results can be generally obtained in less than 20 iterations. Note that three kinds of color are used to show the results clearly.
We finally compare results of the proposed model with greatly relevant models on one synthetic image and one natural image from BSD qualitatively and show the result in Figure 7. Note that to be fair, initializations on either image are all the same for each of the comparable models. And we set the parameters , , and . As the CV, RSF, LIC, and LINC models are short of the capability to extract interested objects from color images directly, we first convert the color image to a gray image using the rgb2gray function of matlab and then input the image to the models. However, the proposed model can be directly used to deal with color images (three channels). Therefore, segmentation contour of the proposed model on the natural image given is marked on the original colorful image whereas results of the others are marked on the gray images. As shown in Figure 7, due to the absence of dealing with intensity inhomogeneity, segmentation results of the CV model include other regions besides geometrical shapes really exist in the synthetic image and eagles in the natural image. Segmentation results of the RSF model are a little better than those of the CV model because it can handle intensity inhomogeneity to some extent. But the RSF model lacks the capability of bias estimation and correction. As shown in Figure 7, the bias fields estimated by the LIC model are obviously not smooth enough and segmentation results are certainly wrong. Although segmentation results and bias estimations of the LINC model are desirable, color images cannot be directly input into the model before being converted to gray ones. In addition, as mentioned in Section 1, convolution operation in the evolution results in a heavy computational burden for LINC which we will further discuss in Section 5.2. It is obvious that the proposed model achieves the best segmentations, bias estimations, and corrections.
4.2. Evaluation on Public Image Repositories
In this subsection, we evaluate the effectiveness of the proposed model quantitatively on one simulated MR dataset and one real MR image dataset. The first one consists of 9 cases of MR images with three different levels of noise and intensity inhomogeneity, respectively. Resolutions of the images are with 1 mm in-plane pixel size and 1 mm slice thickness. For more information about the dataset, interested readers are referred to the website http://brainweb.bic.mni.mcgill.ca/brainweb/ and the reference [44]. To construct a much more challenging dataset for segmentation methods, three more levels of nonlinear intensity inhomogeneities are added to the original image with noises. Therefore, there are totally 18 image cases for the first image dataset. The second image set is known worldwide as the Internet Brain Segmentation Repository (IBSR) which contains 18 cases of T1-weighted brain MR image cases with skull-removed masks and manually guided expert segmentation results. Resolutions of the images are all . Interested readers are referred to https://www.nitrc.org/projects/ibsr for detail. Note that for each image case, the segmentation task is to extract white matter (WM), gray matter (GM), and cerebrospinal fluid (CSF) from the background. As intensities of the background are all zero for the images, two level set functions are used to partition the images into three regions that is and . To compare the performance of the proposed model with representative models like CV, LIC, and LINC on these image datasets, we first extend the comparable models to three phase (matlab codes will be released on our personal homepage if this paper got published). We then define membership functions , , and to represent WM, GM, and CSF, respectively. For a fair comparison, we first extend comparison models to three phases and then use the same parameter set and the same strategy to initialize the level set functions for all the comparison models. The initialization strategy is that areas separated by a predefined threshold are adopted to initialize by considering the areas as inside and outside of the zero-level contour. Areas separated by another predefined threshold are adopted to initialize the level set function . The thresholds are adaptively defined as 0.8 and 0.3 times of maximal intensity of pending to be segmented images. We have to point out that the proposed model is robust to initialization which will be discussed in Section 5.3. Note that we applied the proposed model and comparable representative models only on image slices that really contain WM, GM, and CSF.
4.2.1. Qualitative Comparison
Segmentation results of the proposed model and three representative methods on the 90th slice of selected BrainWeb cases and the 128th slice of selected IBSR image cases are given in Figures 8 and 9. The corresponding bias estimation and correction results are given in Figures 10 and 11, respectively. The reason we select these images is that they are the most noisy and biased and they are therefore challengeable. It can be seen that (1) the proposed model is much more robust to noises and bias fields and (2) segmentation results of the proposed model are much closer to corresponding ground truth. Quantitative evaluation will be given in Section 4.2.2.
4.2.2. Quantitative Evaluation
To quantitatively evaluate segmentation results of the proposed framework with the representative method, false positive ratio (FPR), false negative ratio (FNR), and dice similarity coefficient (DSC) are metrics used in this paper. Let NFP and NFN be the number of FP (false positive) and FN (false negative) and be the ground truth, FPR and FNR can then be defined by and respectively. Pairwise vertical mouldings denote the size of the contained region. As well known, values of FPR and FNR are both in with a smaller value indicating a better match between the segmentation and the ground truth. On the other side, the definition DSC can be written as where is the intersection operator. Values of DSC are in the interval of with a higher value indicating a better match between the segmentation result and the ground truth .
Quantitative comparison of segmentation results of the proposed model with representative models on the BrainWeb and IBSR images in terms of FPR, FNR, and DSC is given in Figure 11 on the left column and right column, respectively. For the BrainWeb dataset, it can be seen that boxes of WM, GM, and CSF of the proposed model in terms of FPR and FNR are much more compacted and the mediums are lower than CV, LIC, and LINC which indicates segmentation results of the proposed model match better with corresponding ground truth than representative models. On the contrary, boxes of WM, GM, and CSF of the proposed model in terms of DSC are also compacted but the mediums are higher than CV, LIC, and LINC which indicates a better match of the segmentation results with corresponding ground truth. On the other side, for the IBSR dataset, FPR boxes of WM and FNR boxes of GM and CSF are more compacted and lower than representative models. DSC boxes of WM, GM, and CSF are more compacted than other models with medium values similar to CV but higher than LIC and LINC. As shown in Figure 12, biases of IBSR images are weak than BrainWeb and ground truths in IBSR images consider more nonzero area as gray matter and therefore decrease areas of WM and CSF. This is the main reason that performance of the proposed model on IBSR is worse than that on BrainWeb images. It has to be pointed out that we set to suppress the areas considered as WM by the proposed model and impact of weighting coefficients will be discussed in Section 5.4.
5. Discussion
5.1. Relationship with CV and PS
We can figure out that when , the proposed energy in Equation (14) can be simplified as the first term of the CV model as shown in Equation (2). To ensure the bias field to be 1 on the entire image domain , the necessary and sufficient condition is that (1) and (2) and for . This indicates that the proposed model is a generalization of the Chan-Vese model. In a similar way, if we define , the energy defined in Equation (14) will reduce to the first term of the PS model as shown in Equation (3) and the smoothness of are ensured by the orthogonal primary functions , , …, and implied in . Therefore, no further regularization term like the second term in Equation (3) are needed to smooth .
5.2. Improvement to LINC
As described in [41], in the case of two-phase implementation of the LINC model, there are 7 convolutions in the size of normalized kernel for each iteration of the level set function, which are the main factor causing computational burden of LINC. As smoothness of bias fields existing in images with inhomogeneous intensities can be guaranteed by orthogonal primary functions, the proposed model GINC removes the convolution kernel from the LINC model and therefore there is no convolution in iterations of the level set function anymore.
5.3. Robustness of GINC to Initialization
As mentioned above, the proposed model is a generalization of CV and a simplification of LINC. It is well known that the intensity constants in CV can be seen as the global average of inside and outside regions separated by the 0-level set contour. Therefore, CV is greatly nonsensitive to local intensities and robust to initialization [30]. On the other side, as pointed out in [41], LINC is also robust to initialization. Thus, as a generalization of CV and a simplification of LINC, the proposed model is robust to initialization. We give a demonstration of the proposed model on one vessel image in four initialization strategies in Figure 13 to verify robustness of the proposed model to initialization. It is obvious that there are no obvious differences between any two strategies in terms of bias estimation and final segmentations, which proves that the proposed model is robust to initialization.
5.4. Impact of Weighting Coefficients
For three-phase segmentation of the proposed model on BrainWeb and IBSR datasets where , , and , the formulation in Equation (20) can be rewritten into
It is obvious that in Equation (20) and where . Therefore, the first term on the right hand of Equation (30) is monotone increasing for and and decreasing for , respectively, only if they take positive values. Thus, given a positive increment on and , the level set function will be increased much harder in each iteration. On the contrary, given a positive increment on , will be decreased much harder. As described in Algorithm 1, we let the level set functions take negative and positive values inside and outside the 0-level set contours, respectively. Hence, for all the others fixed, the greater the coefficient and are, the smaller the region inside the 0-level set contour is, and vice versa. Similarly, the greater the coefficient is, the smaller the region outside the 0-level set contour is, and vice versa. The same analysis can be applied to Equation (31) to conclude that the greater the coefficient and are, the smaller the regions inside and outside the 0-level set contour are, and vice versa. As mentioned earlier, the regularization term and arc length term are used to maintain regularity of the level set function and smooth 0-level set contour. Thus, the greater the parameters and are, the level set function is closer to sign distance function and the smoother the 0-level set contour is.
6. Conclusion and Future Work
The proposed model is effective in segmenting images with inhomogeneous intensities and provide a smooth bias estimation of the inhomogeneity. We will further improve the proposed model to extract brain tissues in 3D on public image repositories in our future work.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (grant No. 61602101) and the Fundamental Research Funds for the Central Universities of China (grant No. N2016003).