Level Set Based Hippocampus Segmentation in MR Images with Improved Initialization Using Region Growing

The hippocampus has been known as one of the most important structures referred to as Alzheimer's disease and other neurological disorders. However, segmentation of the hippocampus from MR images is still a challenging task due to its small size, complex shape, low contrast, and discontinuous boundaries. For the accurate and efficient detection of the hippocampus, a new image segmentation method based on adaptive region growing and level set algorithm is proposed. Firstly, adaptive region growing and morphological operations are performed in the target regions and its output is used for the initial contour of level set evolution method. Then, an improved edge-based level set method utilizing global Gaussian distributions with different means and variances is developed to implement the accurate segmentation. Finally, gradient descent method is adopted to get the minimization of the energy equation. As proved by experiment results, the proposed method can ideally extract the contours of the hippocampus that are very close to manual segmentation drawn by specialists.


Introduction
Magnetic Resonance Imaging (MRI) is a noninvasive medical imaging technology which provides anatomical information on disease-related for detection of pathology and treatment planning. Many clinical applications depend on the segmentation results of MRI brain structures that allow us to know our brain anatomy changes during the deterioration of disease. However, accurate and reliable automatic segmentation of medial temporal lobe structures, such as the hippocampus, is still a challenging. Hippocampus, being one of a small temporal lobe brain structure, has been found to be highly related to many neurological pathologies and psychiatric disorders including Alzheimer's disease [1,2], epilepsy [3,4], schizophrenia [5,6], and mild cognitive impairment [7]. Unfortunately, at the spatial resolution of MRI, the structure of hippocampus has relatively low contrast and no distinguishable boundaries, as shown in Figure 1. Traditionally, clinical hippocampus segmentation almost relies on the manual, which is rater-dependent and highly time-consuming. So, how to segment hippocampus automatically in medical image analysis is critical in clinical aspects.
In recent years, many scholars have been engaged in hippocampus segmentation algorithm, and various types of works also have been done. For example, Kim et al. [8] built a surface model of hippocampus and used principal component analysis to find the shape asymmetry of hippocampus between schizophrenic patients and healthy controls. Wang et al. [9] proposed a region growth algorithm based on seed which is simple and effective but could fail due to the indistinct boundary of hippocampus. Tong et al. [10] proposed a new method for nearly automated detection of the hippocampus contours on MRI images based on combined sparse coding techniques and discriminative dictionary learning. Hajiesmaeili et al. [11] utilized contraction wave atom as an efficient method for enhancing the noisy of MR images to improve segmentation accuracy. Kwak et al. [12] presented graph cuts algorithm to segment hippocampus. Similarly, Wolz et al. [13] applied a graph cuts algorithm to simultaneously measure hippocampal atrophy in longitudinal studies. Another hotspot of segmentation approach is level set method. It allows integration of prior information within a single optimization framework and has proved to be powerful tools in image segmentation. One popular example of edgebased methods is the Geometric Active Contour (GAC) model [14], which used the image gradient to stop evolving contours on object boundaries. The Chan-Vese model [15], which based on the Mumford-Shah segmentation framework [16], is one of the most popular region-based models to detect objects whose boundaries are not defined by strong edges. Recently, Li et al. [17] proposed a new edge-based variational level set method without reinitialization. It incorporated a penalty term to preserve the level set function close to a signed distance function. However, these previous works usually present some problems such as oversegmentation and lower speed.
In this paper, we propose a level set framework combined with seeded region growing algorithm for the segmentation of hippocampus. The method has several advantages: (1) the outcome of the region growing approach is provided automatically as the initial contour of level set evolution method; (2) the global Gaussian distribution with different means and variances is integrated into level set framework. By effectively combining region growing and level set algorithm, the segmentation accuracy is improved while the processing time is decreased. The remaining part of the paper is organized as follows. In Section 2, the proposed method and its theoretical background is presented. In Section 3, results and discussion are given. In Section 4, we state our conclusion.

Methods
Level set method is defined as energy minimization process which can be applied for image segmentation. However, traditional level set algorithm has some commonly known weaknesses: (1) the initial contour is difficult to extract; (2) for high intensity inhomogeneity and discontinuous boundaries images, the edge stop function will produce high values and the contour may be attracted to false edges. In this section, an accurate and efficient level set algorithm is presented. Firstly, the evolution of the level set model is initialized using region growing as the prior information. Besides that, we integrated the global Gaussian distributions information into a conventional edge-based level set model, proposed in [17]. The role of the edge-based is to guide segmentation towards apparent boundaries, while the global Gaussian distributions term is to take lead on regions with weak boundaries. The detail of the proposed segmentation approach is described in Sections 2.1 to 2.2. The flowchart of the proposed method is shown in Figure 2.

Region Growing.
Traditional level set algorithms can handle the segmentation problem accurately by using some contour or image dynamics such as curvature, intensities, or derivatives. But all approaches need the initialization step so that the automatic operation of the system is difficult to achieve. Consequently, some researchers have used shapebased methods to segment the hippocampus and proposed local solutions strictly depending on the shape information of the interest. In this section, the outcome of the region growing approach is provided automatically as the initial contour to the distance regularized level set method. The contour based Li's method [17] is thus aided by region growing approach and improves the segmentation performance compared to both the region growing and Li's models.
The region growing [18][19][20] is a simple region-based image segmentation method that attempts to segment images into many small regions based on predefined seed points, grow rule, and stop condition. As the name is implied, region growing starts from the seed points to the adjacent pixels according to the similarity between the pixel and the marked region and then determines whether the pixel neighbors should be added to the regions. Obviously, selecting seed points and the growth rules is the key of regional growth. However, this method may encounter difficulties with images in which regions both are noisy and have blurred or indistinct boundaries. Therefore, in this study, we proposed a novel and adaptive region growing method for semiautomatically generating the preliminary boundary of hippocampus in MRI. Detailed process will be introduced as follows.
Step 1. Select initial seed points: in MRI images, the hippocampus structure has similar gray values with its adjacent anatomical structures and it is only a small region. In order to reduce the computation and improve the segment reliability, we only process a local area, including the whole hippocampus and little other region. Firstly, an initial seed point is manually selected on the original image of hippocampus, as shown in Figure 3(a). Then, we obtain slice image, with the size of 45 * 45 pixels, of which center is the initial seed point position, as shown in Figure 3(b). After shearing, the majority of the cerebrospinal fluid, white matter, and gray matter can be filtered out, so that it can reduce the interference of the hippocampus segmentation. From the gray histogram shown in Figure 4, we can see that the slice images can better reflect the characteristics of the hippocampus and improve the segmentation accuracy.
Step 2. Region Growing: given the seed, regions are grown from it by appending those neighboring pixels whose properties are most similar to the region. According to the special position of hippocampus in MR image, we choose the standard deviation as the growth rule. It can be described as follows: where ( , ) is the slice image, seed is the mean intensity of current seed points, Global is the global standard deviation of the slice image, and is the weight of Global . This criterion allows including all connected pixel blocks, obtained by the quad-tree decomposition, of which intensity does not differ from the seed point intensity more than a user-defined threshold value. Based on this observation, Global can be defined as where × is the size of slice image and is the mean intensity of slice image. If any pixel belonging to the seed point's 8-connected neighborhood satisfies the growth criterion, we add it to the seed point set. The experiment results of our region growing algorithm are presented in Figure 3(c).
Step 3. Morphological operations: morphological operations are a collection of nonlinear operations related to the shape or morphology of features in an image. Firstly, the perimeter of the homogeneous region obtaining by region growing is extracted, as shown in Figure 3(d). However, the result is not the hippocampus region we desire to extract exactly. Therefore, level set method is introduced to modify the result which is produced by region growing. In order to prevent falling into local minima, a convex polygon is gained by fitting the perimeter and used to initialize the level set algorithm, as shown in Figure 3(e).

Level Set
where , > 0, ] are vary constants and is a level set function. The ( ) is the penalizing term, ( ) is the length of the interface, ( ) is the area of the subregion enclosed by the zero level curve, the function (⋅) and (⋅) are the Heaviside and Dirac function, and is an edge indicator function, defined by where is the Gaussian kernel having the standard deviation .
Due to the penalizing term ( ), the involving level set function driven by (3) can automatically close to a signed distance function; thus reinitialization procedure is eliminated completely. However, this method above has at least three intrinsic limitations in applications. First, it is highly sensitive to the contour initializations; therefore we have to choose suitable initial contours to correctly detect object of interest in given images. Second, it is prone to getting "trapped" by extraneous edges caused by image noise or concavity caused by high intensity variations near these points. Third, it only can detect objects with edges defined by gradient.

The Proposed Level Set Algorithm.
As discussed in the above section, Li's method is highly sensitive to strong noise Computational and Mathematical Methods in Medicine 5 and cannot extract edge without gradient or cognitive edges. For the accurate and efficient detection of the hippocampus, we use the region growing as the prior information and integrate the global Gaussian distributions information into a conventional edge-based level set model. In image segmentation, the role of global information plays a leading on regions with weak boundaries. Particularly it is crucial to reduce the sensitivity to the initialization of the contour. From statistical point of view, in the proposed model, the fitting energy is described by a combination of Li's model and global Gaussian distributions with different means and variances, respectively. The global Gaussian distributions fitting energy is defined as follows: where Ω 1 = in( ), Ω 2 = out( ), and ( ( )) is where and are global intensity means and standard deviation, respectively.
Let Ω ⊂ 2 be a two-dimensional image space and let : Ω → be a given gray image. We assume that the image domain can be partitioned into two regions. These two regions can be represented as the regions outside and inside the zero level ; that is Ω 1 = { > 0} and Ω 2 = { < 0}. Using the Heavside function , the global Gaussian distribution fitting energy function can be rewritten as follows: In this study, we define the energy function as follows: ( ( )) = ( ( )) + ( ( )) + ] ( ( )) where is weighting constants of GDF .
From (9), we obtain which minimize the energy functional ( ) for fixed . Minimization of the energy functional (8) with respect to can be achieved by solving the gradient descent flow equation:

Algorithms and Numerical Approximation.
In this subsection, we briefly present the numerical algorithms and procedures to solve the evolution (12). Due to the regularization term ( ), (12) in the continuous domain can be discretized by means of simple finite difference rather than complex upwind difference scheme. We recall first the usual notation. Let ℎ be the space step, let Δ be the time step, and let ( , ) = ( ℎ, ℎ) be the grid points. Let , = ( Δ , , ) be an approximation of the level set function ( , , ), with ≥ 0, 0 = 0 , where 0 is the initial level set function. The central differences of spatial partial derivatives are expressed in the following notations: All the spatial partial derivatives / and / are approximated by the central difference, and the temporal partial derivative / is discretized as the forward difference. Set = 0 and start with initial level set function 0 , ; knowing , , we first compute 1 and 2 according to (13) and (14). Then, the numerical approximation to (12) is given by the following discretization: where the corresponding curvature , can be discretized using a second-order central difference scheme: with Next, a direct implementation of the proposed level set algorithm is presented in Algorithm 1.

Results and Discussion
To assess the performance of the proposed method, quantitative and visual experiments have been carried out on three groups of the whole brain MRI images (the image size is 157 * 189 and 2D axial brain images with T1-weighted sequence of size 0.78 * 0.78 * 2.00 mm). For each image set, 17 to 19 frames which correctly resemble the shape of hippocampus were used to make the model. These images are from West China Hospital, Sichuan University, China. Li's model and the method of this paper are applied to the hippocampus brain segmentation, while the experiment results were compared. Two methods set the same parameters and initial conditions as follows: the level set function can be simply initialized as a binary value − 0 inside the region and a positive value 0 outside the region. Unless otherwise specified, we use the follow-  Figure 5 offers the segmentation results of five hippocampus images, where the red curve surrounded is the segmentation area. Our method requires a seed point located inside the targeted structure by the user, as shown in Figure 5(a). The segmentation results by Li's model and our method are shown in the second row and the third row, respectively. The last row is the results of artificial segmentation by specialist. The experiment results with the same seed points show that when the images are badly under the conditions of intensity inhomogeneity, low contrast, and discontinuous boundaries, only using the edge information may fail to discriminate the intensity between an object and background, leading to inaccurate segmentation. However, our method both uses the edge information and region information of images, making our method have a very strong discriminative capability for the object and background. So, the segmentation results of this paper are more accurate, which can distinguish regions similar intensity means but different variances.

Segmentation of Hippocampus Images.
Furthermore, in order to show the superiority of our model, we utilize the dice similarity coefficient (DSC) [27] to evaluate the performances on three different group images. If where (⋅) indicates the number of pixels in the enclosed region and Ω is the image domain. The value of DSC ranges from 0 to 1, with a higher value representing a more accurate segmentation result. A comparison plot of DSC for the three groups of the whole brain MRI images which correctly resemble the shape of hippocampus is provided in Figure 6. By quantitative comparison, we can see from the overlap ratios that there is a very good correspondence between our method and manually extracted boundaries and this demonstrates the significant advantage of our model in terms of segmentation accuracy.
In addition, Figure 7 shows the results of our method on other brain tissue images. The first to last rows show temporal lobe, amygdale, and midbrain, respectively. It can be observed that the good result can be obtained by our method.

Robustness to the Initial Seed Point.
In order to evaluate the robustness of our method, we apply it to the hippocampus images in Figure 5 with different initial seed points, as shown in Figure 8. In Figure 8  Image with different initial seed points Figure 9: Segmentation accuracy of our method for images with different initial seed points.
tial seed points and the corresponding segmentation results are shown. In spite of huge difference of these seed points, the corresponding results are nearly the same, and the object boundaries are accurately obtained. To further illustrate the robustness to initialization of our method, the segmentation accuracy is quantitatively verified with DSC values in Figure 9. It can be seen that the -axis represents five original images of hippocampus with four different initial seed points, and the -axis represents the DSC value of each initialization. From the results we can see that the segmentation results obtained by each initial seed point are almost the same. Thus, experiment proved that the proposed model has a good robustness to the initial seed point.

Conclusion
In summary, we have introduced an approach to extract the pathway of hippocampus based on level set algorithm with an improved initialization using region growing, which allows us to optimize the location of hippocampus in a global sense. This model is derived from Li's model by incorporating the global Gaussian distributions with different means and variances into level set framework, avoiding the problems of boundary leakage low contrast and discontinuous boundaries of hippocampus images. Meanwhile, to solve the problem of initialization sensitivity, we have discussed and analyzed the adaptive region growing and morphological operations which are used to gain a rough segmentation result of hippocampus to initialize the level set function. Compared with Li's model, our method provides more accurate and efficient segmentation results that are closer to the manual segmentation results obtained by a specialist. In future work, we will further improve the segmentation efficiency by GPU acceleration and study the adaptive change of the shape constraint's effect according to the quality of the hippocampus images to acquire better segmentation results.