Segmentation of Intensity Inhomogeneous Brain MR Images Using Active Contours

Segmentation of intensity inhomogeneous regions is a well-known problem in image analysis applications. This paper presents a region-based active contour method for image segmentation, which properly works in the context of intensity inhomogeneity problem. The proposed region-based active contour method embeds both region and gradient information unlike traditional methods. It contains mainly two terms, area and length, in which the area term practices a new region-based signed pressure force (SPF) function, which utilizes mean values from a certain neighborhood using the local binary fitted (LBF) energy model. In turn, the length term uses gradient information. The novelty of our method is to locally compute new SPF function, which uses local mean values and is able to detect boundaries of the homogenous regions. Finally, a truncated Gaussian kernel is used to regularize the level set function, which not only regularizes it but also removes the need of computationally expensive reinitialization. The proposed method targets the segmentation problem of intensity inhomogeneous images and reduces the time complexity among locally computed active contour methods. The experimental results show that the proposed method yields better segmentation result as well as less time complexity compared with the state-of-the-art active contour methods.


Introduction
Image segmentation is a fundamental problem in the areas of computer vision and image processing. It is used to partition an image into two or more than two nonoverlapping regions based on textual, intensity, or gradient information. Image segmentation is a particularly difficult task for numerous reasons. Firstly, partitioning the image into nonoverlapping regions and extracting regions of interest requires a tradeoff between the simplicity of algorithm, selection of parameters, computational efficiency of algorithm, and accuracy of the segmentation result. Secondly, image artifacts, such as noise, intensity inhomogeneity, artifacts involved with the image acquisition, and poor contrast of image, are very difficult to account for in segmentation algorithms without high level of interactivity from the user. Different methods are devised in a context of the segmentation problem and each of them has their own advantages and disadvantages. Some of the common techniques used for image segmentation are thresholding based segmentation, segmentation based on image classification, and edge based and region based (region growing) image segmentation.
Active contour is one of the devised techniques for image segmentation problem, which segments an image by evolving a level set curve. In late 1980s, Kass et al. introduced one of the image segmentation techniques based on active contour [1]. In this method, a curve evolves toward the object boundary under a force, until it stops at the boundary. To be more specific, the curve moves toward the object boundaries by minimizing the energy. The energy functional is based on different image characteristics, for example, image gradient, curvature, and image statistical properties.
The existing active contour models are categorized into two groups: edge-based models [1][2][3][4] and region-based models [5][6][7][8][9][10][11][12][13][14]. Both of these types have their own paybacks and drawbacks, and the choice between them to use in applications depends on the different characteristics of the images. The edge-based model builds an edge indicator function using image edge information, which can drive the contour towards the object boundaries [2]. The edge indicator function based on the image gradient can hardly stop at the right boundaries, for the images with intense noise or a weak edge.
On the other hand, a region-based model uses statistical information to construct a region stopping function that can stop the contour evolution between different regions. One of the early efforts towards a region-based active contours was made by the Mumford and Shah segmentation model [5], which approximates a given image by a piecewise smooth image. Compared to the edge-based model, the region-based model can perform better for images with blurred edges. The region-based model is not sensitive to initialization of the level set function and can recognize the object's boundaries efficiently. Therefore, region-based models, especially the Chan and Vese (CV) model [6], have been widely applied for image segmentation.
Although the region-based model is better than edgebased model in some aspects, it still has limitations. The traditional region-based models [5,6], which were proposed in the context of binary images with the assumption that each image region is statistically homogeneous, do not work perfectly for images with intensity inhomogeneity. Figure 1(a) shows an image with white background which contains intensity inhomogeneous region in it and Figure 1(b) shows the ineffectiveness of traditional region-based active contour method in case of image with intensity inhomogeneity.
The traditional region-based models [5,6], which were proposed in the context of binary images, do not work well if the target image contains intensity inhomogeneous regions in it. Li et al. [11,12] proposed the LBF model by embedding the local image information. LBF model is able to segment images with intensity inhomogeneity and is much more accurate than the previously formulated methods. The basic idea of LBF was to introduce a Gaussian kernel function in the energy functional formulation. Although it segments well the images with intensity inhomogeneity, it has quite high computational time complexity. Therefore, segmentation process takes quite a time as compared to old segmentation methods. Zhang et al. [13] proposed an active contour method driven by local image fitting (LIF) energy, which provides almost same segmentation results and has less time complexity as compared to LBF model.
In this paper, we proposed a region based active contour method which works well under the intensity inhomogeneity problem. The proposed region based active contour model utilizes both edge and region information to segment an image into nonoverlapping region. It is implemented by replacing the edge indicator function in the area term of the edge-based level set method [3] with a new regionbased signed pressure force (SPF) function that utilizes the image local information obtained using the local binary fitted (LBF) energy model. By introducing the SPF function based on local fitted image (LFI), the formulated method is able to segment images with intensity inhomogeneity. In the proposed model, region information used in the area term helps to stop the contour at weak or blur edges while edge information in the length term accelerates the detection of those weak or blur edges in corporation of that area term. As the introduced model contains both edge and region information it works better than the traditional edge-based and region-based methods.
Reinitialization, a technique for occasionally reinitializing the level set function to a signed distance function (SDF) during the evolution, has been extensively used as a numerical remedy for maintaining stable curve evolution and ensuring desirable results. From a practical viewpoint, the reinitialization process can be quite convoluted and expensive and has subtle side effects [15]. Zhang et al. [7] proposed the active contour with selective local or global (ACSLG) segmentation method which uses a Gaussian kernel to regularize the level function after each iteration. It not only regularizes the level set but also removes the need of reinitialization. In Computational and Mathematical Methods in Medicine 3 the proposed algorithm we also use the Gaussian kernel to eliminate the need of reinitialization. Regularization using Gaussian kernel has better smoothing results and no energy leakage as compared to the area smoothing and penalization terms used by Li et al. [3].
The proposed segmentation algorithm is applied to synthetic and brain MR images in order to demonstrate the accuracy, effectiveness, and robustness of the algorithm. A comparison is shown with previous related methods to show advantages of the proposed method.

Active Contour Method Driven by Locally Computed Signed Pressure Force (LCSPF) Function
Li et al. [11,12] proposed the LBF energy model by employing the local image information. LBF is able to segment image with intensity inhomogeneity and provides much more accurate results than the traditional region-based methods. The basic idea is to introduce a kernel function to define the LBF energy functional. Let : Ω → be an input image and let be a closed curve, for that the LBF energy functional, LBF ( , 1 , 2 ), is defined as follows: Minimizing the above energy functional with respect to 1 and 2 using steepest gradient descent method [16] we get the following formulations: In image segmentation, active contours are dynamic curves that move toward the object boundaries to partition an image into nonoverlapping regions. To achieve this goal, we explicitly define an energy functional that can move the zero level curve toward the object boundaries. Figure 2 illustrates the above assumptions and notations on the level set function defining the evolving curve , where at boundary of curve value of = 0 and our level set function moves inwards or outwards, based on the signs of the SPF for the further evolution. We define energy functional containing an edgebased length term and a region-based area term for function as follows: where > 0 and V are constants and the terms ( ) and spf ( ) are defined as follows: respectively. Here, = is the univariate Dirac function and is the Heaviside function defined in (9) and (10), respectively, while ( ) is edge indicator function and spf( ) is locally computed SPF function defined in (11) and (14), respectively. The zero level curve is driven into a smooth curve from a complicated curve to minimize the function ( ) which utilizes edge information in regularization process, while spf ( ) contains the locally computed image intensity information which derives the contour to the weak and blur edges by distinguishing inhomogeneous regions.
The energy ,spf ( ) drives the zero level set toward the object boundaries. The coefficient V of spf ( ) in (4) can be positive or negative, depending on the relative position of the initial contour to the object of interest. For example, if the initial contours are placed outside the object, the coefficient V in the weighted area term should take a positive value, so that the contour can shrink faster. If the initial contours are placed inside the object, the coefficient V should take a negative value to speed up the expansion of the contours. By the calculus of variations [16], the Gateaux derivative (first variation) of the functional ,spf ( ) in (4) can be written as The function that minimizes this functional satisfies the Euler Lagrange equation ,spf / = 0. A classical iterative process for minimizing the function is the gradient flow with artificial time given as

Computational and Mathematical Methods in Medicine
After evolving the level set using (7) and (8) we smooth it by using = 2 * . It not only regularizes the level set function but also eliminates the need of reinitialization, which is computationally very expensive. Moreover, it gives energy leakage free reinitialization.
In the proposed work, the Dirac function and Heaviside function used in (2), (3), and (7) are the smoothed version of the Dirac function and Heaviside function of the entire region. The approximations and as proposed in [6] are and the edge indicator function ( ) is a positive and strict decreasing function defined follows: In outmoded level set methods, it is essential to initialize the level set function as a signed distance function (SDF) 0 . If the initial level set function is expressively different from the SDF, then the reinitialization schemes are unable to reinitialize the function to the SDF. In our formulation, not only is the reinitialization procedure completely eliminated but also the level set function no longer needs to be initialized as an SDF. The initial level set function 0 is defined as where > 0 is a constant and we use = 1. Finally, the principle steps of the algorithm can be summarized as follows.
(iii) Compute local mean values 1 and 2 from (2) and (3), respectively, where 1 is the standard deviation of the truncated Gaussian kernel used to compute 1 and 2 . (iv) Calculate SPF function spf( ) using (14).
(vi) Regularize the level set function by a Gaussian kernel; that is, = 2 * , where 2 is standard deviation of a Gaussian kernel.
(vii) Check whether solution is stationary and if not, go to step (ii), = + 1, and repeat.

Locally Computed SPF Function
The SPF function defined in [17] has values in the range [−1, 1]. It modulates the signs of the pressure force inside and outside the region of interest so that the contour shrinks when outside the object and expands when inside the object. Traditional SPF is formulated using global properties of image; therefore, it works poorly with intensity inhomogeneous images. Here, we introduce a new SPF function based on the local properties of image inside and outside of the contour. This newly formulated SPF function formulates the signs of the signed pressure force function inside and outside the boundary of the region of interest using locally computed mean values. A local fitted image formulation is defined as Using the above defined local fitted image we construct the SPF function as follows: The terms 1 and 2 are defined in (2) and (3), respectively. The SPF function computed using the local properties of image is shown in Figure 3 The sign and value of SPF function ranges in [−1 1] for both local and global SPF functions; the only difference is the construction method used. As mentioned earlier global SPF function uses global mean values inside and outside the contour; that is why it cannot distinguish between the inhomogeneous changes in the intensity. Therefore, it assigns values with the same sign to both inner and outer Computational and Mathematical Methods in Medicine 5 region when the intensity change inside and outside region is not distinctive to the global intensity mean computing function. On the other hand local SPF function uses local mean value inside and outside the contour which helps to distinguish between inhomogeneous changes in the intensity. The local intensity mean values are computed by utilizing the Gaussian kernel. It helps computing local maxima for intensity inhomogeneity region which global model fails to find. It offers different signs for both inside and outside the region although there is inhomogeneous intensity change between inside and outside the region. Figure 4 illustrates the signs and segmentation result comparison of local and global SPF function based active contour models. Figure 4(a) shows the image with initial contour for both local and global cases. Figure 4(b) shows the signs of global SPF function, in which we can see that for the inhomogeneous intensity change in the image the sign assigned to inside and outside the regions are the same, while for distinctive change in intensity by the global mean computing function sign assignment for both inside and outside the region is different. Figure 4

2D Synthetic and Brain MR Image Segmentation Results.
The proposed method is implemented using MATLAB 7.12, in Windows 7 environment on a 2.4 GHz Intel Quad-Core personal computer with 8 GB of RAM. The range of intensities in all images is represented from 0 to 255, while the size in pixels (length × width) of images is variable for all images. In this section we applied the proposed method to synthetic and real images of different modalities and used the parameters which are = 1, V = 22, = 1, = 1.5, 1 = 5, 2 = 1, = 5, and = 1, where V is force term constant which controls the contour evolution speed and time complexity of desired contour. 1 is standard deviation of the truncated Gaussian kernel 1 ( ) with the size 4 + 1 by 4 +1. 2 is the standard deviation of the smoothing Gaussian kernel which is used to regularize the level set. Selection of 1 and 2 may be different for different types of images. If we select small values for 1 and 2 then contour will evolve faster but it will not be accurate; that means small values of 1 and 2 can reduce the time complexity but can decrease the accuracy. For large values of 1 and 2 time complexity will increase but segmentation accuracy will also increase. In case of noisy images selection of 2 should be bigger than normal to smooth the level set curve. Figure 5 shows the segmentation result on a synthetic image with nine different intensities. Figure 5(a) is the initial contour; Figure 5(b) is the segmentation result of synthetic image without noise. Obviously, these objects with different intensities both homogeneous and inhomogeneous are successfully extracted because the proposed method also works well with the intensity inhomogeneity. We then add Gaussian noise to the clean synthetic image. The noisy image that is shown in Figures 5(c) and 5(d) shows the corresponding segmentation result of our method on the noisy image. From the segmentation results obtained from both clean and noisy synthetic image we can see that segmentation result for both clean and noisy image is similar, which implies that the proposed method segments well under the dense Gaussian noise as applied in this case. Figure 5(e) displays the central row intensity profile of the input synthetic image with both clean and noisy data along with the final contour. It shows that irrespective of data the resultant contour followed the edges perfectly. In Figure 5(e) we normalized the intensity scale to [−1 1] in order to visualize data from input image profile and final contour profile at the same time with same peak values. The number of iterations used to evolve contour from initial to final form is 200. Figure 6 shows the importance of 2 in the proposed model using a real brain MR image. Figure 6(a) is the initial contour, Figure 6(b) is the final contour with 2 = 1.0, and Figure 6(c) is the final contour with 2 = 0.5. From Figures 6(b) and 6(c) we can see that if we choose large value of 2 then the proposed method does not segment small details, while the selection of small 2 makes segmentation algorithm more sensitive to noise. By using the small value of 2 we can segment more detailed objects. For the blurry images 2 should be small and for noisy image 2 should be large. The total number of iterations used in the contour evolution process is 150.
Segmentation of brain MR image into disjoint regions based on white matter (WM), gray matter (GM), and cerebral spinal fluid (CSF) is a well-known problem. Due to the geometric complexity of the human brain cortex, manual slice by slice segmentation is quite difficult and time consuming [18]. Numerous methods of image segmentation have been developed to solve such problems [19]. Active contour method is one of those methods which are used in this context. Because of the complex intensity inhomogeneous regions, brain MR images are hard to segment successfully with high accuracy [20]. The proposed method is formulated in the context of intensity inhomogeneity problem. To show the robustness and effectiveness of the proposed algorithm for inhomogeneous images we applied it on 2D real brain MR images with intensity inhomogeneity. Figure 7 shows brain MR image segmentation results using the proposed method.

3D Brain MR Image Segmentation Results.
In this section segmentation results of different brain regions are displayed using 3D brain MR models [21] by applying the proposed method. The range of intensities in all images is represented from 0 to 255, while the size in voxels (length × width × height) of images is (217 × 260 × 362). We have chosen the models of five different regions of head which are involved in 3D brain MR scan of a human test subject [21]. Figure 8 shows the segmentation using five anatomical models of human subject in which Figure 8(a) shows the initial contour, Figure 8 We can remove this artifact by applying smoothing kernel on the data input before using the segmentation algorithm. We have applied the proposed method on 3D MR Dataset to show its application in volume visualization and data exploration.

Comparison with Traditional Active Contour Methods
Using SPF Function in Their Model. Zhang et al. in [7] and Jiang et al. in [8] used SPF function in their proposed method but their methods cannot segment well the images with intensity inhomogeneity. The SPF function they used in their model is grounded on CV region-based active contour method which computes mean of intensity globally that is not sufficient in order to segment the images with intensity inhomogeneity, while the proposed SPF function uses local mean value which provides well segmentation results for images with intensity inhomogeneity but with drawback of high time complexity [22]. In order to show advantages of the proposed method over other active contour methods which utilize traditional global SPF function in their models, we compare their results using a synthetic image with intensity inhomogeneity. The parameters used for this comparison for  The proposed method can segment image with intensity inhomogeneous regions yet it takes additional processing time as compared to the active contour methods using traditional CV based SPF function. Table 1

Comparison with the LBF, LIF, and CV Energy Models.
In this paper we proposed a region-based active contour method using locally computed SPF function. The proposed method works similar to previously developed LBF [11,12] and LIF [13] energy models. To show the effectiveness, accuracy, and robustness of the proposed algorithm we compared the segmentation results with the LBF, LIF, and CV energy models. For the comparison we used both synthetic and real brain MR images. The parameters used for this comparison for the LBF energy model are 1 = 1, 2 = 1, = 1, V = 0.001 × 255 2 , = 1, = 1.5, = 4, = 5, and = 0.1. The parameters used for the LIF energy model are 1 = 5,    Figure 10 shows a comparison between the proposed method, LBF, LIF, and CV energy models using synthetic image. Figure 10(c) shows that the final contour computed by the LIF energy model could not strictly follow the object boundaries. Figure 10(d) shows that the LBF energy model could not accurately segment the small inhomogeneous objects. Figure 10(e) shows that the CV model cannot properly segment intensity inhomogeneous regions; moreover, small objects are also missed in the segmentation process. Figure 10(b) shows that the proposed method provided better segmentation from the entire state-of-the-art active contour methods used in the comparison. Figure 11 shows a comparison between the proposed method, LBF, LIF, and CV energy models using three brain MR images. For the comparison shown in Figure 11, parameters for the proposed method are same as mentioned in Section 4.1. The parameters of CV energy model are same as described earlier in this section, while for the LBF and LIF energy models = 5 and 1 = 5, respectively, and the rest of the parameters are same as mentioned earlier in this section. Figures 11(b), 11(g), and 11(l) show the segmentation results using the proposed algorithm, which, compared to the generated results by other methods, are better in every aspect. There is no over lapping in the region during the segmentation process and contour accurately evolved to the boundary of the object need to be segmented. The proposed method even segmented the sharp details of the objects at the boundary.  Figure 11(n) there are also some signs of contour overlapping during segmentation process. Figures 11(e), 11(j), and 11(o) display the segmentation result using CV energy model. It shows that the CV energy model is unable to segment the intensity inhomogeneous regions in the brain MR images. Table 2 shows the time complexity analysis of the proposed method, LBF, LIF and, CV energy models from Figures  10 and 11. In Figures 10 and 11, the proposed method has less time complexity compared to LBF and LIF energy models. For Figure 10  segment images with intensity inhomogeneity; the proposed method being second in time complexity is still the best among remaining methods. 1 plays an important role in time complexity and segmentation accuracy. If we select same value of 1 = the proposed method has less time complexity and better segmentation results compared to the state-of theart locally computed active contour methods.

Quantitative Analysis.
As discussed earlier white matter (WM) and gray matter (GM) are the main regions of interest in the segmentation of brain MR images. In order to segment WM and GM we split the segmentation result into two regions based on two phases. WM region represents the phase at > 0 and the GM represents the phase at < 0. The WM and GM regions represent the brain region which is the region of interest, while the regions outside the brain, for example, skull, fats, and vacuum, can be taken as unnecessary regions. Therefore, we have used a hand drawn brain mask to extract the WM and GM only and remove the other needless regions outside the brain. Figure 12 shows the computed WM and GM from the proposed method phases and compares it with the given ground truth. It shows how we extracted only GM and WM regions by scaling it to the brain area using hand drawn brain mask. Figure 12(a) shows the initial contour at = 0. After some time = , we obtained the final contour shown in Figure 12(b). Then we divided the final contour into two phases based on the value of ; at > 0 we acquired WM region and at < 0 we obtained GM region as shown in Figures 12(c) and 12(g), respectively. Figures 12(d) and 12(f) show hand drawn brain mask. After obtaining WM and GM regions from two phases of final contour we then multiply these regions with hand drawn brain mask in order to scale them to brain area only and remove other regions outside the brain, for example, skull, fats, and vacuum, as shown in Figures 12(e) and 12(i), respectively. Finally, we compare the scaled WM and GM regions computed using the proposed method with their respective ground truths to visually analyze the segmentation accuracy of the proposed method.
In order to do the quantitative analyses we used the ground truths of 2D slices from the 3D anatomical brain models [21]. The quantitative analysis for the proposed method, LBF, and LIF models is shown in Table 3 using the following expression for the percentage accuracy: We compute the accuracy by using the given ground truth data and the computed segmentation results. In the above expression is the ground truth of the WM or GM region and is the brain mask scaled WM or GM region from > 0 or < 0, respectively, using the proposed method, LBF model, and LIF model.  Figure 12: Visual based segmentation accuracy analysis of WM and GM regions. (a) Initial contour; (b) final contour; (c) WM region of final contour which is phase at > 0; (d) hand drawn brain mask; (e) WM region after brain mask scaling; (f) WM region ground truth; (g) GM region of final contour which is phase at < 0; (h) hand drawn brain mask; (i) GM region after brain mask scaling; (j) GM region ground truth.
accuracy for both WM and GM regions as compared to LBF and LIF models. Selection of standard deviation of truncated Gaussian kernel plays an important role both in the time complexity and segmentation accuracy of the algorithm. If or 1 is big then contour evolves faster (less number of iterations) to its final form with much of segmentation accuracy but its time complexity also increases. If or 1 is small then contour evolves slower (more number of iterations) to its final form but there can be segmentation accuracy problem. 2 is used to regularize the contour and remove the need of initialization. The bigger 2 is the smoother the contour would be at the boundary of the object to be segmented.

Conclusion
In this paper a region-based active contour method is presented which embeds both edge-based and region-based terms in its model. As the proposed model contains both edge-based and regions-based terms, it works better than traditional region-based methods and segments well images with weak and blur edges. A new SPF function is formulated which utilizes image local information and helps to segment intensity inhomogeneous regions. The proposed method is applied to the 2D synthetic and real brain MR images with intensity inhomogeneity to show its robustness and effectiveness. We also applied the proposed method to 3D brain anatomical models to show its application in volume visualization and data exploration.
A comparison is shown with other active contour methods that use traditional global SPF function formulated by CV model, using synthetic image with intensity inhomogeneity. It shows that the proposed method accurately segments images with intensity inhomogeneity unlike previously formulated SPF based active contour methods. However, it has higher time complexity than the active contour methods using traditional global SPF function.
We also compared the proposed method with the previously formulated locally computed active contour methods to show the advantages of the proposed method. The visual comparison shows that the proposed method generates better segmentation results as compared to the state-of-the-art active contour methods for both synthetic and real brain MR images. Moreover, it also has less time complexity compared to the local based active contour methods.
Global region based active contour method is fast as compared to the local active contour method but it does not work well for the images with intensity inhomogeneity. On the other hand local active contour method using local maxima can segment images with intensity inhomogeneity but it has high time complexity and it is sensitive to noise. In future we will address the segmentation problem with an active contour method by using a new SPF function that uses both local and global information. By using that we want to target advantages of both local and global region based models.