Hybrid two-stage active contour method with region and edge information for intensity inhomogeneous image segmentation

This paper presents a novel two-stage image segmentation method using an edge scaled energy functional based on local and global information for intensity inhomogeneous image segmentation. In the first stage, we integrate global intensity term with a geodesic edge term, which produces a preliminary rough segmentation result. Thereafter, by taking final contour of the first stage as initial contour, we begin second stage segmentation process by integrating local intensity term with geodesic edge term to get final segmentation result. Due to the suitable initialization from the first stage, the second stage precisely achieves desirable segmentation result for inhomogeneous image segmentation. Two stage segmentation technique not only increases the accuracy but also eliminates the problem of initial contour existed in traditional local segmentation methods. The energy function of the proposed method uses both global and local terms incorporated with compacted geodesic edge term in an additive fashion which uses image gradient information to delineate obscured boundaries of objects inside an image. A Gaussian kernel is adapted for the regularization of the level set function and to avoid an expensive re-initialization. The experiments were carried out on synthetic and real images. Quantitative validations were performed on Multimodal Brain Tumor Image Segmentation Benchmark (BRATS) 2015 and PH2 skin lesion database. The visual and quantitative comparisons will demonstrate the efficiency of the proposed method.


Introduction
Image segmentation is the well-defined scheme of partitioning an image into several regions and it has a great significance in image processing and computer vision [1]. Intensity variations and inhomogeneities often occur in medical images mostly due to the imperfect image acquisition process and influence of the external surroundings. The existence of inhomogeneities and complicated boundaries create difficulties in image segmentation and eventually, it leads to uncertainty for doctors and radiologists during the decision-making process. Active contour models have been popular techniques for those type of images [2]. These methods provide smooth contours for image segmentation by using some energy minimization principle [3]. Existing active contour models are mainly categorized as edge-based [3][4][5][6][7][8] and regionbased [9,10,[12][13][14][15] methods. Edge-based methods exploit image gradient information to deploy a balloon force used to capture object boundaries in curve evolution [5]. Nevertheless, these methods are very sensitive to noise and unable to segment objects in the presence of fuzzy or blurred boundaries. Alternatively, region-based methods essentially intend to distinguish each region by using the statistical and curvature information to control the movement of contour, thus they have better execution over noisy and blurred images. These methods are further categorized into a global [9][10][11][12] and a local [13][14][15][16] region based active contour methods.
One of the most popular global region method is Chan-Vese method [12], which is a special case of a piecewise constant Mumford and Shah method [10]. This method relies on constant intensity information across the regions therefore, it fails to segment images with intensity inhomogeneity. In order to overcome the limitation of the piecewise constant methods, Li et al. [14] proposed a local region-based method. This method has introduced a kernel function that considers the local image information known as LBF (local binary fitting) energy to segment images with intensity inhomogeneity. Segmentation of inhomogeneous image is shown in Fig 1, where Fig 1(a) shows segmentation of traditional Chan-Vese (global region method) and Fig 1(b) shows segmentation of LBF (local region method). From Fig 1, it can be deduced that LBF method is capable of handling intensity variations across the regions and give better result in terms of intensity inhomogeneity.
Recently, hybrid models [2,[17][18][19][20][21][22] have gained popularity among region based methods. These methods either combine region (local and/or global) and edge information or both local and global intensity information in their energy formulations. In [17], Zhang et al. proposed a method which combines the advantages of an edge-based and region-based active contours. It defines a SPF (signed pressure force) function using global means from Chan-Vese [12] method. A Gaussian kernel is also adapted to regularize level set and to avoid re-initialization. Oh et al. [20] has also formulated a hybrid method using scaled edge and region information. It combines global region information from Chan-Vese [12] method and coupled geodesic edge term from GAC [5] method to prevent boundary leakage problem.
Alternatively, A region-based active contour methods are proposed [22,23] in the context of local and global intensity information. These methods integrate local and global intensity information in their energy formulations. Global force drives the movement of contour while local force captures the inhomogeneities across the regions.
A variational two-stage image segmentation approach for intensity inhomogeneity is proposed by Wang et al. [21]. Two stages of this method were based on global and local region active contours described by Gaussian distributions. The final segmentation of the first (global) stage is served as the initial contour of the second (local) stage. In this way, this method is able to segment images corrupted by intensity inhomogeneity.
This paper presents, a novel two-stage hybrid edge and region(local and global) active contour method for image segmentation with a level set formulation. Firstly, we obtain the coarse segmentation result of the first stage and the final segmentation result in the second stage. In the first stage, we incorporate global region force term and geodesic edge term for statistical and edge information. Global intensity used to segment homogeneous regions while edge term used to prevent boundary leakage problem. In the second stage, we develop a modified energy function which integrates local binary fitting term and geodesic edge term. Minimization of an energy functional in the second stage mutually segments the objects with intensity inhomogeneity and restrict the contour at object boundary. Since the local fitting methods [14,15] are very sensitive to the initial position of contour. Proposed method takes final contour of the first stage as initialization of the second stage thus, it provides a better choice for initial contour position and give satisfactory performance for intensity inhomogeneous image segmentation.
In level set methods, it is often required to initialize the level set function to SDF (signed distance function) and avoid it to be too steep or too flat near the interface. In this regard, reinitialization is required to maintain the stability of level set function. There are numerous operational methods to implement re-initialization schemes. Traditional variants of the level set piecewise constant methods [11,24,25] eliminate the re-initialization by solving the constrained optimization problems. Liu et al. [26] proposed Lagrangian constrained optimization method to eliminate the re-initialization in level set function.
Proposed method adapts a Gaussian kernel from Zhang et al. [17] which not only regularizes the level set but also avoid an expensive re-initialization. Proposed method has been executed over several real and synthetic images including brain and skin lesion datasets and compared with previous methods. Experimental results will determine that proposed method offers an improved scheme and achieve better segmentation results than previous methods.
Rest of this paper is organized as follows. The background study is discussed in section 2. The proposed method is explained in section 3. Final results and comparisons are illustrated in section 4 using both real and synthetic images. The quantitative results are shown in section 5 using the public dataset from Multimodal brain tumor image segmentation benchmark (BRATS) 2015 [27] and PH 2 skin lesion [28] database. In the end, conclusion and future work are described in section 6.

Background and theoretical justifications
Geodesic active contour model GAC method [5] is regarded as a standard active contour method which uses image gradient information from the boundary of the object. Let I: O & R 2 is an image domain, I: O ! R is an input image and C(q) is a closed curve. They proposed the following energy functional: where g is the edge stopping function defined as: where rG σ Ã I is the convolution of an image I with a Gaussian kernel whose standard deviation is σ. By minimizing the above energy functional we get following Euler-Lagrange equation: where k denotes the curvature of the contour andÑ is inward normal of the curve. The final level set equation is defined as follows: This method relies on edge based contour evolution which can only capture objects with edges defined by its gradient. Therefore, this method doesn't support regional information and fall into local minimum when the initial contour is not placed near object boundaries.

Mumford and Shah model
Mumford and Shah [10] proposed a region-based image segmentation method. This method finds an optimal piecewise approximation function μ of the image I, which keep changing smoothly within a sub-region of image domain I: O & R 2 . The proposed energy formulation of this method is written as: where L(C) is the length of the curve, μ and v are positive parameters. Practically, the energy functional of this model owns non-convex behavior and the non-regularity of the edge term causes difficulty during energy minimization.

Chan-Vese model
In [12], Chan-Vese proposed a simplified energy functional based on Mumford and Shah model [10] by approximating the image intensities inside and outside of contour known as c 1 and c 2 respectively. Let I: O & R 2 be an input image, ϕ: O & R 2 a level set function and C a closed curve corresponding to the zero level set: C = {x 2 O|ϕ(x) = 0}. The energy functional of Chan-Vese method is defined as follows: where μ ! 0, v ! 0 and λ 1 , λ 2 ! 0 are scaling constants, constant μ ! 0 scales the Euclidian length of the curve C and constant v scales the area term inside the contour C. H ε (ϕ) is the regularized Heaviside function defined as: where handles the smoothness of Heaviside function. In Eq (6), two constants c 1 and c 2 , represent global region intensities inside and outside of contour C, respectively. Minimizing Eq (6), with respect to ϕ using gradient descent method [29], the corresponding level set formulation is obtained as follows.
δ (ϕ) is a smooth version of a Dirac delta function, which is defined as: beside handling the smoothness of a Heaviside function in Eq (7), also controls the width of the a Dirac delta function in Eq (9). Keeping ϕ fixed and minimizing energy function in Eq (6), with respect to c 1 and c 2 , we have: The curve evolution process in Chan-Vese method is related to global characteristic of an image region inside and outside curve C, respectively. Therefore, this method produces unacceptable result if the image has local or inhomogeneous regions inside and outside of the curve C.

Local binary fitted model
In, [14] Li et al. proposed local binary fitted (LBF) method by encapsulating local image information into their energy functional to deal with intensity inhomogeneity problem. Let, I: O & R 2 be an input image, ϕ: O & R 2 is a level set function, and C is a closed curve. They propose following energy formulation: where λ 1 , λ 2 ! 0 are positive parameters, H ε (ϕ) is the regularized Heaviside function defined in Eq (7). f 1 (x) and f 2 (x) are locally approximated intensities inside and outside of contour C defined as.
In order to ensure stable solution, the distance regularization term from [4] is fused in to Eq (12). Moreover, the Euclidean length term is also added to regularize the zero curve of ϕ. The total variational formulation of this method is defined as: where μ is a scaling parameter for distance regularized penalty term and v is scaling parameter for length term which drives the movement of contour towards object boundaries. K σ is a Gaussian kernel with a standard deviation which is defined as: where σ is a standard deviation of the Gaussian kernel function which controls the degree of localization from small neighborhood to whole image domain. The introduction of Gaussian kernel function considers an image local information inside and outside of contour C and enable this method to segment images with intensity inhomogeneity. however, local image information is not always sufficient to carry out an accurate segmentation. Moreover, this method is extremely sensitive to the position of initial contour and stuck into local minima if we place initial contour far from boundaries.

ESRAC model
Oh et al. [20] proposed a hybrid region and edge method known as edge scaled region-based active contour method. This method redefines Chan-Vese [12] region information weighted by edge stopping function and integrate geodesic edge term [5] into their formulation. They propose following energy formulation: where g(x) is an edge stopping function scaled with region force term. parameter v used to scale the area term inside and outside the contour C and α is scaling parameter whose values ranges between 0 α 1, it controls the balance between region and geodesic edge term. Minimizing Eq (17), using gradient descent method [29], the corresponding level set formulation is obtained as follows.
In above equation g is monotonically decreasing edge function defined as: The regional term in above equation is weighted with edge indicator function therefore, this method can segment object boundaries by using gradient information. In the case of insufficient gradient information this method can also segment object by using statistical region information. However, this method is not available for inhomogeneous intensity regions.

Proposed model
Inspired by the work of Oh et al. [20] and wang et al. [21] we propose a novel two-stage active contour segmentation method by integrating the region-based global and local methods with geodesic active contour method. In the first stage, we acquire segmentation result from the global region-edge active contour method (GREAC). This method rapidly captures all homogeneous regions and spurious boundaries in an image and achieves its coarse segmentation result. Later in the second stage, we use final contour of the last stage as initialization curve and continue segmentation process by minimizing the energy functional of the local region-edge active contour method (LREAC). The initial stage gives a partial segmentation result and provide a suitable initialization for second stage, in this way this method eliminates the common problem of initial contour existed in local region based methods [13,14]. Finally, the second stage considers an image local information and gets the final segmentation result.

First stage
Global region-based segmentation has been regarded as a fast and stable segmentation for homogeneous regions and it takes global image intensities inside and outside of contour. Furthermore, it is not sensitive to initialization of contour. In the first stage, we propose an energy functional of the global region-edge active contour method (GREAC) defined as: In Eq (20), c 1 and c 2 are edge scaled intensity means defined as:: where g(|rI|) is a monotonically decreasing edge stopping function defined in Eq (2). By using steepest gradient descent method [29], minimization of Eq (20), with respect to ϕ leads to the following gradient decent flow: The energy functional in above equation is scaled with a parameter w whose value ranges between 0 w 1. This parameter manages the weight of the global region and geodesic edge term depending on the type of image we are dealing with. When w is large the global region term act as the main force with small influence of edge term to detect global regions and boundaries inside an object. Similarly, when w is small geodesic force term act as the main force with small influence of global region information to attract the curve towards the desired boundary of the object. This process is shown in Fig 2(a) which describes the role of global and geodesic energy terms during the first stage.

Second stage
The first stage of the proposed method uses global image information across the regions. Therefore, it is unable to capture objects with intensity inhomogeneity. In this stage we extend global region-edge active contour method (GREAC) to local region-edge active contour method (LREAC) by using local binary fitting energy from LBF [14] method. The energy functional of the local region-edge active contour method (LREAC) is defined as: þ ð1 À wÞð f 1 (x) and f 2 (x) are locally approximated intensities inside and outside of contour C defined as.
where g(|rI|) is a monotonically decreasing edge stopping function defined in Eq (2). By using steepest gradient descent method [29], minimization of Eq (24), with respect to ϕ leads to the following gradient decent flow: The energy functional in above equation is also scaled with a parameter w whose value ranges between 0 w 1. This parameter manages the weight of the local region and geodesic edge term depending on the type of image we are dealing with. When the w is large the localregion term act as the main force with small influence of edge term to detect inhomogeneous regions and boundaries inside an object. Similarly, when w is small geodesic force term act as the main force with small influence of local-region information to attract the curve towards the desired boundary of the object. This process is shown in Fig 2(b) which describes the role of local and geodesic energy terms during the second stage.
Traditional let set methods need to initialize their level set function (ϕ) periodically to SDF (signed distance function) during the contour evolution. Li et al. [4] proposed a penalization term which is also computationally expensive. The proposed method uses a Gaussian filter which regularizes a level set function and also avoids an expensive re-initialization. The initial level set function ϕ 0 for the proposed method is defined as: where ρ ! 0 is a constant, O 0 is region inside an initial contour, O is the image domain and @O is an initial contour. Finally, the iterative steps of the proposed method are summarized in following algorithm: Algorithm

Results and quantitative comparisons
All the experiments were implemented in MATLAB on a personal computer with Intel core i5, 3.2 GHz, and 8 GB RAM. We have used parameters for all experiments which are listed in Table 1.  [14] method and LCV [16] method since these methods are very sensitive to initial contour therefore, these methods have not properly delineated the required object boundary. Column (e) shows the result of Wang et al. [21] method. This method is able to properly segment the second image, however, some unwanted contour has been produced in a fourth and fifth image. Moreover, this method is unable to detect obscured edge in the first image. Column (f) shows the result of the proposed method which is able to capture and segment all the objects properly.   achieved desired segmentation results. On the other proposed method has achieved an accurate segmentation results with less number of iterations and CPU time compared to LBF [14], LCV [16] and Wang et al. [21] methods respectively. The Local region based active contour methods [14,16] are able to distinguish small changes between the background and the foreground. Therefore, they are suitable for intensity inhomogeneous images. However, there exists an intrinsic drawback of initial contour in local region based active contour methods. Therefore such methods are not sufficient for segmenting inhomogeneous information alone. Proposed method eliminates the initial contour problem by providing an ideal contour position in the first stage from global region based segmentation. In second stage local method gets the advantage of an ideal contour position and capture all inhomogeneous regions efficiently.
In Fig 5, we have shown the influence of initial contour position on the segmentation result. Column (a) shows five different level set initializations on the same image, Chan-Vese [12] result is shown in column (a), LBF method [14] results are shown in column (c), LCV method [16] results are shown in column(d), Wang et al. [21] method results are shown in column (e)  Two-stage active contours and proposed method results are shown in column (f) respectively. Results demonstrate that proposed method has not intruded by the position of initial contour, in contrast previous methods have shown their sensitivity to the initial contour and produced undesired results.
In Table 3, we have recorded the number of iterations and CPU time of each method based on Fig 5. The proposed method have conceded the second least number of iterations and CPU time after Chan-Vese [12] method and achieved desired segmentation results, whereas the Chan-Vese method was unable to get required results. Though, Wang et al. [21] method has also produced good results in the first and second row however it consumed a large number of iterations and CPU time for obtaining segmentation result, while proposed method has obtained the same results in all four initialization examples within less number of iterations and CPU time.
The experiment has also been performed on a real set of images with inhomogeneous and complex region properties. Fig 6 shows segmentation of three different real images, where column (a) shows the initial contour with original images. Column (b), (c), (d), (e) and (f) shows the result of Chan-Vese [12], LBF [14], LCV [16], Wang et al. [21] and proposed method respectively. All previous methods either capture uninterested regions or stuck in the background and thus have failed to get the desired segmentation results. On the contrary, proposed method gets the required results smoothly. Table 4     proposed method in second, third, fourth, fifth and sixth column respectively. The CPU times of these images are listed in Table 5. By comparing the computational time of the proposed method and previous methods, it is clear that proposed method is also modest and proficient when it comes to same accuracy level with previous state-of-the-art methods.

Quantitative analysis
The frequency of skin cancer has been consistently expanding over the past decade. Melanoma has been considered as one of the dangerous forms of skin cancer, reports show that survival rate of this cancer is very less during its last stage [30]. Therefore, detection of melanoma is very important and it helps dermatologists for initial treatment and diagnosis. Active contours have become a hot topic for melanoma image segmentation therefore, proposed method has been also validated over skin lesion dataset and compared with the previous state of art methods. In this regard, we use publically available skin lesion dataset PH 2 [28]    applied the proposed method and previous state of art methods on 200 skin images and calculated the average accuracy of each method. The formula for accuracy analysis is defined as: where A is the segmented region and B is the ground truth. Fig 9 shows that proposed method gets the high accuracy compared to the previous state of the art methods. MRI (Magnetic Resonance Imaging) is regarded as most commonly used noninvasive modality for brain tumor detection [31]. Manual segmentation of brain tumor by radiologists and experts suffer from errors mostly due to tiredness and fatigue. Segmentation of abnormal regions in MRI is essential because it assists doctors and specialists to study the growth of a tumor and quantitatively analyze it over a period of time. Many segmentation approaches have been proposed for brain tumor segmentation. Recently active contours have gained much popularity among brain tumor segmentation [32]. Therefore, we have also validated our method on the BRATS 2015 challenge dataset [27]. Four MRI sequences are available for each patient known as: T1-weighted(T1), T1 with gadolinium-enhancing contrast (T1c), T2-weighted (T2) and FLAIR. We have tested our method on BRATS 2015 [27] training dataset which comprises of 200 HGG(high-grade glioma) and 44 LGG(low-grade) patient volumes. We use three validation metrics for quantitative evaluations which are Dice coefficient (DSC), sensitivity(Sen) and specificity (Spe). The obtained result will be considered as good when the measured values of these metrics are close to 1. Sensitivity metric value defines that all detected regions(tumors) are correct and belongs to the ground truth. similarly, specificity specifies that none of the healthy tissue is considered as tumor tissue. Moreover, Dice coefficient measure how much detected tumor region overlaps the ground truth. These metrics are defined as: Spe  [27]. The first column show the original images with initial contour and the second column shows the ground truth followed by the result of (from left to right) Chan-Vese [12], LBF [14], LCV [16], Wang et al. [21] and proposed method. It has be observed that final segmentation of the proposed method is robust and accurate against the complex background, where almost all previous methods fail to delineate the right boundaries. Moreover, Fig 10 show the segmentation results for three planes: axial, coronal and sagittal as shown in ground truth(second column). Results show that proposed method has accurately segmented brain tumor, Necrosis, and surrounding edema tissues successfully. Fig 11 summarize the comparative results of each method on BRATS 2015 challenge dataset, we compute the average sensitivity, specificity, and Dice index of each method against the manual ground truth. It shows that proposed method has demonstrated the capability and best rate of segmentation. The Chan-Vese method could not manage to get desired results, LBF and LCV methods have good Dice index values but their sensitivity and specificity values are not satisfactory, Wang et al method perform better with values close to the ground truth in some cases however it is unable to delineate proper boundaries in majority cases finally, proposed method on the other hand consistently outperformed traditional methods and achieved good results.

Conclusion
In this paper, we have presented a hybrid novel two-stage segmentation method with application to synthetic and real (including medical) images combining global, local and edge intensity information. The first stage gets the intensity homogeneous segmentation result by employing global and edge information, after reaching to maximum number of iterations in the first stage we get the contour near the object boundaries. The second stage takes the final contour from the first stage and employs the local intensity and edge information to get the final accurate result for intensity inhomogeneous image segmentation. By using the concept of the two stages, proposed method provides an appropriate solution to the contour initialization which is a very common problem in traditional local active contour methods. Proposed method scales global and local region terms with edge function which prevents the boundary leakage and over segmentation problem. Finally, a Gaussian kernel is used to regularize the level set function and to avoid an expensive reinitialization.
The main contribution of our method is the formulation of the local and geodesic energy functional based on ESRAC [20] method for capturing intensity inhomogeneous regions. Initially, experimental results have been taken from a set of synthetic and real images in order to show the robustness of the proposed method. Finally, for qualitative and quantitative analysis we have evaluated the proposed method on publically available databse PH 2 [28] for skin lesion and BRATS 2015 [27] for brain tumor segmentation. The comparative analysis shows that proposed method has obtained better results than previous state-of-art methods.