I Introduction
In recent years, medical image segmentation has become a standard technique for visualizing structures of the human brain as well as performing various types of volumetric and shape comparisons among these structures. Since the introduction of medical image segmentation, many methods have been implemented for brain structure segmentation from magnetic resonance imaging (MRI). These methods can be categorized into manual, semiautomatic, and fully automatic methods. Manual segmentation is tedious, requires training and much attention to detail, and the results are not reproducible. On the other hand, fully automatic methods require no training and are completely reproducible for the same data, but these methods do not allow for human intervention or manipulation and severely limit the autonomy of the one performing the segmentation. These issues have caused semiautomatic methods to become the preferred type of medical image segmentation [1].
Semiautomatic segmentation can also be divided into several categories, but the two primary classifications for medical image segmentation include regionbased (region growing, region merging) [2] and boundarybased (snake and balloon) [3, 4, 5] techniques. Regionbased methods provide quick segmentation results by assigning membership to voxels according to homogeneity statistics, but the inhomogeneity among MRI voxel intensities can result in inaccurate segmentation (i.e. holes and irregular boundaries) [2].
Boundarybased methods attempt to align an initial deformable boundary with the object boundary by minimizing an energy functional which quantifies the gradient features near the boundary. This technique works well for images with little noise but these methods are also generally unreliable because image noise and low contrast edges between brain structures can result in false or nonexistent boundaries causing under or oversegmentation [3, 4]. For boundarybased methods, defining the initial geometry prior to deformation (i.e. the seed) is also a critical issue that has yet to be resolved. Without a proper method for seed contour initialization, the seed will deform to local rather than global minima in most circumstances due to image noise. [6, 7] provide evidence for these advantages and drawbacks. These unrefined seeding methods increase the time required to complete each segmentation.
The most effective way to overcome the issues with regionbased and boundarybased methods is to utilize both methods in parallel. [8, 9] used regionbased information to drive the explicit deformable models in their techniques, while [10, 11, 12, 13, 14] have addressed these issues by interlacing regionbased and boundarybased methods into a united, iterative segmentation process. The efficacy of these types algorithms exceed that of regionbased or boundarybased methods independently, but the most notable disadvantages of these methods are that they limited to slicebyslice (2D) segmentation and have relatively low efficiency compared with other segmentation algorithms. Also, texture has been used as a criteria to drive a hybrid deformable model as in [15]. However, this method also lacks an efficient seed initialization technique.
Ii Overview
In this paper, we present a hybrid semiautomatic segmentation method in 3D that integrates both regionbased and boundarybased procedures. The goal of this algorithm is to achieve robust, efficient, and reproducible results while avoiding the downfalls of regionbased and boundarybased procedures methods individually. Unlike previous methods, we use a clustering technique to initialize the deformable model and then deform the model with a unique force equation independent of image edges. Separating the two methods allows for more efficient segmentation in comparison with previous hybrid procedures. By initializing the model with a regionbased technique, we are able to avoid most, if not all, local minima in which the model may fall into during deformation. The speed of segmentation is also substantially improved because of this nearboundary seeding method. We utilize an implicit deformable model to alleviate the problem of topology changes, such as selfintersection [16]. The model is driven by the PDE similar to the one derived in [17], which allows the evolution of the model to be less affected by image noise and lowcontrast edges.
Iii Methods
The algorithm consists primarily of two independent phases: seed initialization and implicit deformation. To initialize the seed, a user must first specify a voxel inside the target structure in the MRI. Next, the MRI is clustered (using KMeans) based on voxel intensity to obtain a general outline of the structure. We refine this seed by performing a mixture of mathematical morphology combined with a connected components search. Specifically, the seed obtained from clustering is first eroded several times. Once erosion is complete, a connected components search originates from the point the user has previously specified, and the voxels that are not found during this search (i.e. not connected to the primary cluster) are removed from seed. This ensures that the seed is now a connected structure, like those in the MRI we wish to segment. A matching number of dilation steps are performed on the seed to ensure that the seed s volume has not diminished to an insufficient size. These steps remove pieces of the cluster that may not be a part of the target structure.
To remove artifacts from the seed and complete the segmentation, the seed is deformed based on an implicit level set PDE. Before deformation, the fast sweeping method for distance field initialization is used to create an implicit representation of the seed. The distance field is a discrete scalar function that determines the distance from a voxel to the nearest point on the boundary of the implicit surface. Once a distance field has been created, the field is deformed according to a PDE that does not rely heavily on image edges, but rather takes into account surrounding voxel information. This type of PDE is used to ensure accurate segmentations on structures that have lowcontrast edges, such as the hippocampus. A narrow band algorithm for level set evolution is also employed during the deformation in order to increase the efficiency of the segmentation. Once the equilibrium of the PDE is reached, the segmentation is complete and an implicit representation of the structure is obtained. We use a standard marching cubes algorithm to obtain an explicit mesh representation of the object for later shape and volumetric comparison. Figure 1 diagrams the main phases of the algorithm as well as the intermediate stages within each phase.
Iiia Seed Generation
One of the most important process of the algorithm is the initialization process in which a suitable seed is determined for the deformable model. For the purposes of this paper, the seed is defined as the initial geometry prior to deformation (i.e. ). The seed created during this initial phase provides a coarse representation of the target shape. We are able to create this seed using a simple clustering technique as well as several iterations of connected components and mathematical morphology.
Let be a mapping such that where is defined by the dimensions of the 3D MRI. In order to ensure that is valued in the interval , we normalize the MRI data by setting all voxel intensities below the percentile to , all intensities above the percentile to
, and linearly interpolating the values in between so that all voxel intensities fall between
and . To initiate the algorithm, a user must specify a point inside a target structure ().The next step is to cluster using the kmeans clustering algorithm as found in [18]. This is done by minimizing the error term :
(1) 
where is the number of clusters to be used, each is a subset of and disjoint from one another such that , and is the average intensity of all voxel intensities in the cluster . Either a mean shift algorithm can be applied to to determine , or a user can specify this value.
Clustering separates the domain, into clusters. The most important cluster is the one which contains
since this cluster provides an initial estimate of the target shape. Thus, the other clusters can be disregarded. Assume that
where . Then, let and . We use to denote the voxels inside of the deformable model, and to denote the voxels outside of the model. This notation is important when discussing the future levelset deformation of the model. The next steps will further refine the seed () by ensuring it’s connectedness and removing voxels from the seed which are not strongly affiliated with the target structure.Voxels which are not affiliated with the target structure will in most cases be in or near the boundary of , denoted . To remove these voxels, is eroded with a mathematical morphology operation a certain number of iterations (dependent on the size of ). This essentially removes from at each iteration, and then recalculates after each step of erosion. In some cases, essential voxels of the target structure are removed from , but the deformation stage (see section 3.2) is used to overcome these seeding artifacts.
As with most mathematical morphology operations, should then be dilated the same number of steps that it was eroded. However, in our implementation, this step is proceeded by a connected components step to ensure the connectedness of because there is no guarantee, and little possibility, that is connected. Since most structures in the brain are connected, we make an assumption that our seed should be connected as well to produce more accurate segmentations. To alleviate this possible unconnectedness, a simple connected components algorithm is implemented originating from the point . During the connected components search, only voxels in are available to be searched. is recalculated based on all voxels that are visited during the search, and the outside set is also recalculated with the equation .
Finally, the set is dilated the same number of steps that it was eroded, recalculating and dynamically. These seed creation steps can be done very efficiently, and allow for the development of a seed that is roughly equivalent to the target structure. Figure 2 shows the progression of the seed initialization phase from start to finish.
IiiB Deformation
Some seeds will be coarse and also contain various artifacts or imperfections due to image noise. To increase the accuracy of segmentation, a subsequent deformation phase is applied to the initial seed. Figure 3 illustrates the necessity and benefits of this subsequent deformation phase.
To begin deformation, the seed, is transformed into an implicit distance field first and will deform based on a levelset PDE. A signed distance function, is initialized such that is the signed euclidean distance (positive on the set and negative on the set ) from to the closest voxel in .
We have found that the most efficient method to perform this distance field initialization is to use a 3D fast sweeping method as described in in [19]. The fast sweeping method works by first initializing to be very large at all points in its domain. Then, directly compute for all voxels that are in or adjacent to with the following equation:
(2) 
The rest of ’s domain is computed by propagating an approximation of the actual distance using elements which have previously been computed. To achieve this, the domain is ”swept” a total of eight times (starting from a new corner of the 3D rectangle and ending at the opposite corner for each sweep) using the below method of computation to solve for the values of . In three dimensions, this calculation is done by solving the following equation at each grid point, assuming a uniform grid size of one:
(3) 
Where is the maximum between and 0, is the value of at the point in question, and are the minimum surrounding values in each direction of the point such that . For a detailed explanation of the fast sweeping algorithm in dimensions, refer to [19]. To make a true signed distance field, it must be ensured that is negative when acting on any voxels in and positive for all values in .
is then deformed based on a PDE similar to those described in [17, 20]. The PDEs in these works focus heavily on using information other than edge detection functions to drive the levelset deformation. The idea behind the deformation is to introduce an artificial time variable and to update the levelset function as time elapses. We update by numerically solving the following PDE at each time step:
(4) 
where and are the arithmetic means of the intensities of all voxels in their respective sets, and , , , and are coefficients. More specifically, determines the weight to be associated with the curvature term, is an external force to be applied to , and and are weights for the distance functions and respectively. It is important to note that the sets and are updated at each iteration to match the way they had been previously defined with respect to , i.e. and . Thus, and will need to be recalculated at each iteration as well.
For the purpose of efficiency, we also implement a narrow band algorithm so that only values of the distance field that are within a certain threshold are updated. For example, while updating the distance field, we only solve the aforementioned PDE near the voxel points where . A similar approach and more detailed implementation can be found in [21].
Once the deformation is complete, the signed distance function now provides us with an implicit representation of the target shape; however, an explicit mesh must be extracted from for further shape analysis and shape comparison. For our implementation, this is done with a standard marching cubes algorithm as discussed in [22]. This algorithm takes a three dimensional signed distance function as input and returns a mesh of triangles and vertices that represents a discrete approximation of the zeroth levelset of the distance function.
Iv Results and Validation
This method has been used to segment various brain structures from actual patient MRI as well as test data. We have successfully segmented the corpus callosum, lateral ventricles, hippocampi, and thalami from these data sets with few failures. Figure 4 shows several examples of completed segmentations from several different orientations and perspectives.
To test the validity and accuracy of the segmentations, we have compared segmented structures using our algorithm to manual (ground truth) segmentations of the structure from the same MRI. Only real patient MRI were used in these experiments to verify that our technique is applicable in clinical research settings. For the comparisons, we calculate volumetric comparison statistics (dice similarity and overlap coefficients) between the two segmentations. The range of these statistics lie between 0 and 1, with 1 indicating a perfect agreement between our segmentation and the manual segmentation. We tested the corpus callosum, left and right ventricles, and left and right hippocampi segmentations using our algorithm against corresponding manual segmentation of these structures performed by a trained expert to obtain the dice similarity and overlap coefficients. Overall, our method was able to closely reproduce the accuracy of the manual segmentations. Table I provides a summary of the results.
Structure  N  

Corpus  16  0.79  0.05  0.79  0.03 
callosum  
Lateral ventricle  8  0.81  0.04  0.73  0.05 
(left)  
Lateral ventricle  8  0.78  0.06  0.81  0.03 
(right)  
Hippocampus  8  0.69  0.08  0.66  0.06 
(left)  
Hippocampus  8  0.70  0.07  0.65  0.05 
(right) 
is the standard deviation for these calculations. Similarly,
is the arithmetic mean of the overlap coefficient calculation for each structure, and is the standard deviation for these calculations. For both the dice similarity and overlap coefficient means, each value ranges from 0 to 1, with larger values indicating better results.We have also compared the efficiency of our segmentation method to a boundarybased method using a generic seed (a sphere centered at the userspecified point). The only difference between these two methods was the way in which the seed was created. The results are shown in Table II, where Method A uses with our seeding method for segmentation, and Method B uses the aforementioned generic seeding method for segmentation. The results show a large increase in efficiency when using the seeding technique described in this paper, with a speedup of nearly 300 to 500 percent based on the structure. Combined with the accuracy our algorithm attains, this efficiency comes at little to no cost to the overall segmentation results.
Structure  Method A  Method B 

Corpus callosum  33s  102s 
Lateral ventricle (left)  54s  263s 
Lateral ventricle (right)  58s  240s 
Hippocampus (left)  47s  142s 
Hippocampus (right)  46s  148s 
V Conclusion and Future Work
Differing from previous hybrid segmentation methods, this paper presents a technique for seeding that can improve the accuracy and efficiency of current deformable model semiautomatic segmentation methods. This initialization method is novel and can substantially increase the segmentation time of modern deformable geometry segmentation, both explicit and implicit. We have also introduced a new method for combining regionbased and boundarybased segmentation methods, which is significantly different than the previous hybrid methods. This semiautomatic segmentation algorithm maintains accuracy while significantly improving the speed of modern boundarybased segmentation; the validation results attest to this. The method works with a number of structures in the brain and is not limited to a specific structure. Also, the algorithm is designed to segment 3D structures outright, as opposed to performing 2D slicebyslice segmentation and stitching the resulting contours together to obtain a 3D shape.
While our algorithm has been able provide consistent results among the data sets we have tested, there have been a few failed segmentations. This was primarily caused by excessive noise in the MRI, or due to low image resolution. In the future, we hope to address these issues by perfecting the seeding and deformation phases, or perhaps adding several MRI preprocessing stages to our algorithm. We would also like to increase the robustness of our algorithm to allow for segmentation of several other complex structures, such as the cerebellum, by incorporating a texture parameter into the PDE in which our deformable model is governed by, such as in [15].
Acknowledgment
This work is supported in part by a NIH predoctoral training grant for Clinical Biodetectives, the Department of Defense Autism Concept Award, and the NARSAD Foundation Young Investigator Award.
References
 [1] Y. J. Hu, M. D. Grossberg, and G. S. Mageras, “Semiautomatic medical image segmentation with adaptive local statistics in conditional random fields framework,” in 30th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 2008, pp. 3099–3102.

[2]
S. C. Zhu, T. S. Lee, and A. L. Yuille, “Region competition: Unifying snakes,
region growing,and bayes/mdl for multiband image segmentation,” in
Proc. Intl. Conf. on Computer Vision
, 1995, pp. 416–423.  [3] M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: Active contour models,” Int l Journal of Computer Vision, vol. 1, pp. 321–331, 1987.
 [4] L. H. Staib and J. S. Duncan, “Boundary finding with parametrically deformable models,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 14, no. 11, pp. 1061–1075, 1992.
 [5] T. Mcinerney and D. Terzopoulos, “Topology adaptive deformable surfaces for medical image volume segmentation,” IEEE Transactions on Medical Imaging, vol. 18, pp. 840–850, 1999.
 [6] L. D. Cohen and I. Cohen, “Finiteelement methods for active contour models and balloons for 2d and 3d images,” IEEE Trans. on Pattern Analysis and Machine Intelligence, vol. 15, pp. 1131–1147, 1993.
 [7] D. Metaxas, PhysicsBased Deformable Models. Kluwer Academic Publishers, 1996.
 [8] A. Chakraborty, M. Worring, and J. S. Duncan, “On multifeature integration for deformableboundary finding,” in Proc. Intl. Conf. on Computer Vision, 1995, pp. 846–851.
 [9] R. Ronfard, “Region based strategies for active contour models.” International Journal of Computer Vision, vol. 13, no. 2, pp. 229–251, October 1994. [Online]. Available: http://perception.inrialpes.fr/Publications/1994/Ron94

[10]
A. Chakraborty and J. S. Duncan,
Integration of boundary finding and regionbased segmentation using game theory
. Kluwer, 1995, pp. 189–201.  [11] T. N. Jones and D. N. Metaxas, “Segmentation using deformable models with affinitybased localization,” CVRMed, 1997.

[12]
——, “Image segmentation based on the integration of pixel affinity and
deformable models,” in
Proc. Conf. Computer Vision and Pattern Recognition
, 1998, pp. 330–337.  [13] ——, “Automated 3d segmentation using deformable models and fuzzy affinity,” in IPMI, 1997, pp. 113–126.
 [14] T. Chen and D. N. Metaxas, “Gibbs prior models, marching cubes, and deformable models: A hybrid framework for 3d medical image segmentation,” in MICCAI (2), 2003, pp. 703–710.
 [15] X. Huang, D. Metaxas, and T. Chen, “Metamorphs: Deformable shape and texture models,” IEEE Computer Society Conference on Computer Vision and Pattern Recognition, vol. 1, pp. 496–503, 2004.
 [16] S. Osher and R. Fedkiw, “Level set methods: An overview and some recent results,” J. Comput. Phys., vol. 169, pp. 463–502, 2001.
 [17] T. A. Chan and L. A. Vese, “Active contour without edges,” IEEE Transactions of Image Processing, vol. 10, no. 2, pp. 266–277, 2001.
 [18] K. Alsabti, S. Ranka, and V. Singh, An Efficient KMeans Clustering Algorithm, ser. Lecture Notes in Computer Science: Methodologies for Knowledge Discovery and Data Mining. Springer, 1999, vol. 1574, pp. 355–360.
 [19] H. K. Zhao, “A fast sweeping method for eikonal equations,” Mathematics of Computation, vol. 74, pp. 603–627, 2005.
 [20] F. Gibou and R. Fedkiw, “A fast hybrid kmeans level set algorithm for segmentation,” in 4th Annual Hawaii International Conference on Statistics and Mathematics, 2005, pp. 281–291.
 [21] A. E. Lefohn, J. M. Kniss, C. D. Hansen, and R. T. Whitaker, “A streaming narrowband algorithm: Interactive computation and visualization of level sets,” IEEE Transactions on Visualization and Computer Graphics, vol. 10, no. 4, pp. 422–433, 2004.
 [22] W. E. Lorensen and H. E. Cline, “Marching cubes: A high resolution 3d surface construction algorithm,” ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, pp. 163–169, 1987.
Comments
There are no comments yet.