A Fast, Semi-Automatic Brain Structure Segmentation Algorithm for Magnetic Resonance Imaging

TL;DR

This paper introduces a hybrid semi-automatic 3D brain structure segmentation algorithm for MRI that combines region-based and boundary-based methods, offering improved efficiency and accuracy over previous techniques.

Contribution

The proposed method uniquely separates region-based and boundary-based steps, enabling more efficient segmentation and robust seed initialization for brain MRI analysis.

Findings

Achieves high accuracy in brain structure segmentation.

Demonstrates increased efficiency compared to existing methods.

Effective seed initialization reduces local minima issues.

Abstract

Medical image segmentation has become an essential technique in clinical and research-oriented applications. Because manual segmentation methods are tedious, and fully automatic segmentation lacks the flexibility of human intervention or correction, semi-automatic methods have become the preferred type of medical image segmentation. We present a hybrid, semi-automatic segmentation method in 3D that integrates both region-based and boundary-based procedures. Our method differs from previous hybrid methods in that we perform region-based and boundary-based approaches separately, which allows for more efficient segmentation. A region-based technique is used to generate an initial seed contour that roughly represents the boundary of a target brain structure, alleviating the local minima problem in the subsequent model deformation phase. The contour is deformed under a unique force equation independent of image edges. Experiments on MRI data show that this method can achieve high accuracy and efficiency primarily due to the unique seed initialization technique.

Tables

Table 1. TABLE I: Comparisons of the semi-automatic segmentation algorithm presented in this paper to manual (ground truth) segmentations. N denotes the number of segmentations used for each structure. μ D subscript 𝜇 𝐷 \mu_{D} is the arithmetic mean of the dice similarity calculation for each structure, and σ D subscript 𝜎 𝐷 \sigma_{D} is the standard deviation for these calculations. Similarly, μ O subscript 𝜇 𝑂 \mu_{O} is the arithmetic mean of the overlap coefficient calculation for each structure, and σ O subscript 𝜎 𝑂 \sigma_{O} 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.

Structure	N	${\mu }_D$	${\sigma }_D$	${\mu }_O$	${\sigma }_O$
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)

Table 2. TABLE II: Comparisons of the efficiency of segmentation using differing seeding methods. Method A was seeded with the technique described in this paper; Method B was seeded by initializing a sphere centered at a user-specified point. All other aspects of the segmentation remained constant. The amount of time each segmentation took is measured in seconds.

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

Equations

$$\mathcal{E}=\sum _{j=1}^k\sum _{\mathbf{x}\in C_j}|I(\mathbf{x})-I({\tilde{\mathbf{x}}}_j)|$$

$$\phi(\bf{x})=\left\{\begin{array}[]{l}\displaystyle{-\min_{\bf{y}\in{\mathcal{\partial\mathcal{I}}}}}\|\bf{x}-\bf{y}\|\;\mathrm{if}\,\bf{x}\in{\mathcal{I}}\\ \displaystyle{\min_{\bf{y}\in{\mathcal{\partial\mathcal{I}}}}}\|\bf{x}-\bf{y}\|\;\mathrm{if}\,\bf{x}\in{\mathcal{O}}\end{array}\right.$$

$$[(x-a_1{)}^{+}{]}^2+[(x-a_2{)}^{+}{]}^2+[(x-a_3{)}^{+}{]}^2=1$$

$$\frac{\partial \varphi }{\partial t}=\alpha \nabla \cdot \left(\frac{\nabla \varphi }{\Vert \nabla \varphi \Vert }\right)-\beta -{\gamma }_1(I-\tilde{\mathcal{I}}{)}^2+{\gamma }_2(I-\tilde{\mathcal{O}}{)}^2$$

Full text

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, semi-automatic, 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 semi-automatic methods to become the preferred type of medical image segmentation [1].

Semi-automatic segmentation can also be divided into several categories, but the two primary classifications for medical image segmentation include region-based (region growing, region merging) [2] and boundary-based (snake and balloon) [3, 4, 5] techniques. Region-based 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].

Boundary-based 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 non-existent boundaries causing under- or over-segmentation  [3, 4]. For boundary-based 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 region-based and boundary-based methods is to utilize both methods in parallel. [8, 9] used region-based information to drive the explicit deformable models in their techniques, while [10, 11, 12, 13, 14] have addressed these issues by interlacing region-based and boundary-based methods into a united, iterative segmentation process. The efficacy of these types algorithms exceed that of region-based or boundary-based methods independently, but the most notable disadvantages of these methods are that they limited to slice-by-slice (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 semi-automatic segmentation method in 3D that integrates both region-based and boundary-based procedures. The goal of this algorithm is to achieve robust, efficient, and reproducible results while avoiding the downfalls of region-based and boundary-based 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 region-based 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 near-boundary seeding method. We utilize an implicit deformable model to alleviate the problem of topology changes, such as self-intersection [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 low-contrast 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 K-Means) 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 low-contrast 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.

III-A 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. $t=0$). 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 $I$ be a mapping such that $I:\mathrm{\Omega}\rightarrow[0,1]$ where $\mathrm{\Omega}$ is defined by the dimensions of the 3D MRI. In order to ensure that $I$ is valued in the interval $[0,1]$, we normalize the MRI data by setting all voxel intensities below the $2^{nd}$ percentile to [math], all intensities above the $98^{th}$ percentile to $1$, and linearly interpolating the values in between so that all voxel intensities fall between [math] and $1$. To initiate the algorithm, a user must specify a point inside a target structure ($\mathbf{x}_{0}\in\mathrm{\Omega}$).

The next step is to cluster $\mathrm{\Omega}$ using the k-means clustering algorithm as found in  [18]. This is done by minimizing the error term $\mathcal{E}$:

[TABLE]

where $k$ is the number of clusters to be used, each $C_{j}$ is a subset of $\mathrm{\Omega}$ and disjoint from one another such that $\bigcup_{j=1}^{k}C_{j}=\mathrm{\Omega}$, and $\mathbf{\tilde{x}}_{j}$ is the average intensity of all voxel intensities in the cluster $C_{j}$. Either a mean shift algorithm can be applied to $I$ to determine $k$, or a user can specify this value.

Clustering separates the domain, $\mathrm{\Omega}$ into $k$ clusters. The most important cluster is the one which contains $\mathbf{x}_{0}$ since this cluster provides an initial estimate of the target shape. Thus, the other clusters can be disregarded. Assume that $\mathbf{x}_{0}\in{C_{m}}$ where $m\in{\{1,\ldots,k\}}$. Then, let ${\mathcal{I}}=C_{m}$ and ${\mathcal{O}}=\mathrm{\Omega}\setminus C_{m}$. We use $\mathcal{I}$ to denote the voxels inside of the deformable model, and $\mathcal{O}$ to denote the voxels outside of the model. This notation is important when discussing the future level-set deformation of the model. The next steps will further refine the seed ($\mathcal{I}$) 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 ${\mathcal{I}}$, denoted $\partial{\mathcal{I}}$. To remove these voxels, ${\mathcal{I}}$ is eroded with a mathematical morphology operation a certain number of iterations (dependent on the size of $\mathrm{\Omega}$). This essentially removes $\partial{\mathcal{I}}$ from ${\mathcal{I}}$ at each iteration, and then recalculates ${\mathcal{I}}$ after each step of erosion. In some cases, essential voxels of the target structure are removed from $\mathcal{I}$, but the deformation stage (see section 3.2) is used to overcome these seeding artifacts.

As with most mathematical morphology operations, $\mathcal{I}$ 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 $\mathcal{I}$ because there is no guarantee, and little possibility, that ${\mathcal{I}}$ 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 $\mathbf{x}_{0}$. During the connected components search, only voxels in $\mathcal{I}$ are available to be searched. $\mathcal{I}$ is recalculated based on all voxels that are visited during the search, and the outside set is also recalculated with the equation ${\mathcal{O}}=\mathrm{\Omega}\setminus\mathcal{I}$.

Finally, the set ${\mathcal{I}}$ is dilated the same number of steps that it was eroded, recalculating ${\mathcal{I}}$ and ${\mathcal{O}}$ 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.

III-B 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, $\mathcal{I}$ is transformed into an implicit distance field first and will deform based on a level-set PDE. A signed distance function, $\phi:\mathrm{\Omega}\rightarrow\mathbb{R}$ is initialized such that $\phi(\mathbf{x})$ is the signed euclidean distance (positive on the set $\mathcal{O}$ and negative on the set $\mathcal{I}$) from $\mathbf{x}$ to the closest voxel in $\partial\mathcal{I}$.

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 $\phi$ to be very large at all points in its domain. Then, directly compute $\phi$ for all voxels that are in or adjacent to $\partial\mathcal{I}$ with the following equation:

[TABLE]

The rest of $\phi$’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 $\phi$. In three dimensions, this calculation is done by solving the following equation at each grid point, assuming a uniform grid size of one:

[TABLE]

Where $(m)^{+}$ is the maximum between $m$ and 0, $x$ is the value of $\phi$ at the point in question, and $a_{1},a_{2},a_{3}$ are the minimum surrounding values in each direction of the point such that $a_{1}\leq a_{2}\leq a_{3}$. For a detailed explanation of the fast sweeping algorithm in $n$ dimensions, refer to [19]. To make $\phi$ a true signed distance field, it must be ensured that $\phi$ is negative when acting on any voxels in $\mathcal{I}$ and positive for all values in $\mathcal{O}$.

$\phi$ 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 level-set deformation. The idea behind the deformation is to introduce an artificial time variable $t$ and to update the level-set function $\phi$ as time elapses. We update $\phi$ by numerically solving the following PDE at each time step:

[TABLE]

where $\tilde{\mathcal{I}}$ and $\tilde{\mathcal{O}}$ are the arithmetic means of the intensities of all voxels in their respective sets, and $\alpha$, $\beta$, $\gamma_{1}$, and $\gamma_{2}$ are coefficients. More specifically, $\alpha$ determines the weight to be associated with the curvature term, $\beta$ is an external force to be applied to $\phi$, and $\gamma_{1}$ and $\gamma_{2}$ are weights for the distance functions $(I-\tilde{\mathcal{I}})^{2}$ and $(I-\tilde{\mathcal{O}})^{2}$ respectively. It is important to note that the sets $\mathcal{I}$ and $\mathcal{O}$ are updated at each iteration to match the way they had been previously defined with respect to $\phi$, i.e. $\mathcal{I}=\{\bf{x}\in\mathrm{\Omega}:\phi({\bf x})\leq\rm{0}\}$ and $\mathcal{O}=\{\bf{x}\in\mathrm{\Omega}:\phi({\bf x})>\rm{0}\}$. Thus, $\tilde{\mathcal{I}}$ and $\tilde{\mathcal{O}}$ 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 $\phi=0$. A similar approach and more detailed implementation can be found in [21].

Once the deformation is complete, the signed distance function $\phi$ now provides us with an implicit representation of the target shape; however, an explicit mesh must be extracted from $\phi$ 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 level-set 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.

We have also compared the efficiency of our segmentation method to a boundary-based method using a generic seed (a sphere centered at the user-specified 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 speed-up 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.

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 semi-automatic 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 region-based and boundary-based segmentation methods, which is significantly different than the previous hybrid methods. This semi-automatic segmentation algorithm maintains accuracy while significantly improving the speed of modern boundary-based 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 slice-by-slice 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 pre-doctoral training grant for Clinical Biodetectives, the Department of Defense Autism Concept Award, and the NARSAD Foundation Young Investigator Award.