Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
IEEE Trans Biomed Eng. Author manuscript; available in PMC 2011 March 31.
Published in final edited form as:
PMCID: PMC3068615

A Hybrid Geometric–Statistical Deformable Model for Automated 3-D Segmentation in Brain MRI

Albert Huang, Student Member, IEEE,corresponding author Rafeef Abugharbieh, Member, IEEE, Roger Tam, and Alzheimer’s Disease Neuroimaging Initiative


We present a novel 3-D deformable model-based approach for accurate, robust, and automated tissue segmentation of brain MRI data of single as well as multiple magnetic resonance sequences. The main contribution of this study is that we employ an edge-based geodesic active contour for the segmentation task by integrating both image edge geometry and voxel statistical homogeneity into a novel hybrid geometric–statistical feature to regularize contour convergence and extract complex anatomical structures. We validate the accuracy of the segmentation results on simulated brain MRI scans of both single T1-weighted and multiple T1/T2/PD-weighted sequences. We also demonstrate the robustness of the proposed method when applied to clinical brain MRI scans. When compared to a current state-of-the-art region-based level-set segmentation formulation, our white matter and gray matter segmentation resulted in significantly higher accuracy levels with a mean improvement in Dice similarity indexes of 8.55% (p < 0.0001) and 10.18% (p < 0.0001), respectively.

Index Terms: 3-D image segmentation, brain segmentation, deformable models, geodesic active contour

I. Introduction

Magnetic resonance (MR) has become the main modality for brain imagingthat facilitates safe, noninvasive assessment and monitoring of patients with neurodegenerative diseases such as Parkinson’s disease, Alzheimer’s disease (AD), epilepsy, schizophrenia, and multiple sclerosis (MS) [1]–[6]. The ability to diagnose and characterize these diseases in vivo using MR image data promises exciting developments both toward understanding the underlying pathologies, as well as conducting clinical trials of drug treatments. One important biomarker that is often used to assess patients with neurodegenerative disease is brain tissue volume. The typical rate of global brain atrophy in MS patients has been shown to be 0.6%–0.8% annually, which is two to three times the normal atrophy rate [7]. Evidence has shown that white matter (WM) and gray matter (GM) atrophy at different rates, and each correlates differently to disability [8]–[10]; thus, accurate measurement of the WM and GM brain tissues can provide valuable quantitative indicators of disease progression and, potentially, treatment outcomes [7], [11]. Thus, the main goal of this paper is to introduce an automatic algorithm for robust WM, GM, and cerebrospinal fluid (CSF) segmentation to facilitate accurate measurement of brain tissues.

Previously, to measure various tissue volumes in MRI head scans, manual WM and GM segmentations were often performed by skilled experts. Manual segmentation, however, is extremely time-consuming, mostly limited to 2-D slice-based segmentation, and prone to significant intra- and interrater variability [12]. In particular, manual segmentation cannot be practically and efficiently performed in situations where precise measurements on a large number of scans are required, such as in clinical trials. Therefore, a fully automatic, highly accurate, and robust tissue segmentation technique that provides systematic quantitative analysis of tissue volumes in brain MRI is an invaluable tool in many studies of neurodegenerative diseases. A wide variety of methods have been proposed for automating the segmentation process over the past decade that provided either semi- or fully automated frameworks for the segmentation of brain tissues. A review of some of these methods can be found in [13] and [14].

One popular family of brain tissue segmentation methods is based on normalizing the brain scans by registering (or aligning) them to a predefined atlas of brain tissues. One example is the popular statistical parametric mapping (SPM) technique, which relies on voxel-based morphometry (VBM) [15]. A number of extensions to the original SPM technique have been proposed. For example, SPM is utilized to initialize an expectation–maximization (EM) segmentation framework [16], which has been extended to nonrigid registration [17]. Although such atlas-based methods are typically robust to artifacts such as acquisition noise and distortions, concerns and discussion [18]–[20] have ensued regarding the use of templates from one population when analyzing data from another population. For example, morphing patient scans with pathologically enlarged ventricles to match a normal template could potentially distort the surrounding tissues in an unpredictable manner. Such structural differences might lead to systematic biases and misregistration errors that are difficult to quantify [19]. Such concerns introduce yet another level of complications arising from image registration and atlas generation procedures that add to the already nontrivial segmentation problem, especially in the presence of anomalies such as tumors, lesions, and tissue atrophy.

A second family of brain tissue segmentation methods assigns a label for each tissue based on image statistics either by clustering [21] or by modeling the brain tissue intensity distributions as a finite mixture of distributions such as EM [22], maximum a posteriori (MAP) [23], simulated annealing [24], and Gaussian mixture modeling (GMM) [25]. Other approaches incorporate additional regional information, which is lacking from these statistical methods, into their segmentation framework. Such methods extend clustering or EM by integrating with fuzzy connectedness [26], topological constraints [27], Gibbs random field (GRF) [28], and hidden Markov random field (HMRF) [29] in the segmentation task. A common difficulty with many of these methods, particularly the random field approaches, is the requirement for proper parameter settings in a supervised setting.

A third family of brain tissue segmentation methods is based on utilizing geometric information such as deformable models or active contours [30] that delineates region boundaries using a minimization of an energy functional [31], [32]. Deformable models employing level sets [33] provide an effective implicit representation rather than explicit parameterization of the evolving contour. However, a common problem of directly applying the active contour approach in segmenting brain MR images is leakage through weak or noisy edges that are ubiquitous, especially for edge-based deformable models, e.g., geodesic active contour [34], which describe the evolution of propagating curve as a function of image gradient features. Some researchers incorporated image statistics into their deformable models in various segmentation applications by utilizing coupled surface principle [35], [36] and fuzzy logic [37], [38] to achieve better stability. Others employed a region-based model [39] by utilizing regional homogeneity in a curve evolution perspective and a hierarchical implementation on brain pathology images [40]. More recently, tissue segmentation was performed and quantitatively evaluated [41]–[43] by using the multiphase 3-D level set segmentation (M3DLS) algorithm [41]. M3DLS utilizes a multiphase extension [44] of the region-based deformable model [39] based on the Mumford–Shah functional [45] by iteratively deforming two closed curves separating four regions. This minimal partitioning approach assumes piecewise constant or piecewise smooth data and optimizes a sum of terms, including the lengths and areas for the two closed curves, and the sums of square intensity differences from the means for all four separated regions. This minimization is also performed in a level set framework [33] implicitly. Further extending this model to N -phase allows separation of 2N regions but the number of classes to be segmented is limited to two to the power of the number of closed curves defined. Moreover, complexity increases as more level sets are required and the rate of convergence is typically slow [46].

In this paper, we propose a MR brain tissue segmentation approach that integrates both geometric and statistical image features into an edge-based deformable model formulation to achieve accurate segmentation results. By utilizing this novel hybrid image feature, we present one solution to the challenging problem of stabilizing the active contour. Similar existing work used a topology preservation principle enforced at non-simple points in a geometric deformable model [47], or a curve shortening prior for smoothness in a level set framework to minimize leakage [40], [48]. Here, we do not explicitly apply any smoothness and topological constraints (e.g., topology preservation at nonsimple points) to the geometric deformable models but rather rely on the proposed hybrid feature to regularize the level sets. Other approaches used prior knowledge such as a distance penalizing term in the level set function between two boundary classes [35], a fuzzy decision system on contour distance to an anatomical target or atlas [37], or a dissimilarity measure between the contour and a shape model in the energy term [49]. Here, we demonstrate our proposed approach in segmenting complex anatomical structures such as WM, GM, and CSF without a priori knowledge. Hence, the proposed approach is truly automated and data-driven in both statistical and geometric sense. Furthermore, we compare the segmentation performance of our proposed edge-based level set method to the region-based M3DLS approach [41] on real clinical MRI scans. We demonstrate the improved WM, GM, and CSF segmentation accuracy and robustness when using the proposed method.

The paper is organized as follows: In Section II, we introduce our novel hybrid geometric–statistical feature implemented on the edge-based geodesic active contour formulation. Our modified deformable model is then used to design a new automated 3-D brain tissue segmentation algorithm for both single and multiple MR sequence data. In Section III, we present quantitative and qualitative results and analysis obtained on simulated and real clinical MRI images, as well as comparisons to results reported by using M3DLS. We then conclude in Section IV.

II. Methods

In this section, we present our 3-D brain segmentation method that integrates both geometric and statistical features in an edge-based geodesic active contour framework. We describe the proposed model in its general form. We then present a segmentation framework for brain MRI for both single and multiple MR sequence data.

A. Edge-Based Deformable Model

In this study, we utilize the geodesic active contour model [34] rather than the region-based formulation [39] due to its computation soundness and extendibility. The geodesic model delineates region boundaries by describing the evolution of a curve or surface C from an initial position C0 as finding the minima of the Riemannian curve distance


where g is a general feature function, |[nabla]I| is the gradient norm of intensity I, and q is the parameterization of the curve C. The right side of the equation describes the parameterized curve C(q) such that the Euclidean length of C can be represented as L(C) = [contour integral operator]|C′ε;(q)|dq = [contour integral operator]ds, where ds = |C′(q)|dq is the Euclidean arc length or the Euclidean metric. Note that this geodesic formulation of the active contour relies on g, the speed and halting feature for the evolving surface in 3-D applications derived based on the geometric gradient feature of an image. Generally, g is chosen as a positive-valued function of the intensity gradient as in (2), where Î is a smoothed version of I and ρ = 1 or 2 [34]. Other similar monotonically decreasing functions, such as the sigmoid function (3) with parameters α (width of intensity window) and β (center of intensity window) are also often utilized [50].


The value of this feature function determines the propagation of the surface by searching for the minimal Riemannian distance. An ideal edge would ultimately have a feature value of zero at all the pixel points along this boundary. However, propagation relying solely on edge feature is typically sensitive to noisy and weak edges that are frequently observed in medical images. In particular, with the presence of complex anatomical structures, it is often impossible to automatically and accurately derive the desired geometric edge term to prevent contour leakage into the surrounding regions. Consequently, achieving accurate segmentation results with edge-based geodesic active contour requires either user intervention or careful adjustment of parameters such that the ideal boundary is minimal. This process is subjective and ideal parameters are often difficult to derive for a fully automated segmentation framework.

B. Proposed Hybrid Geometric–Statistical Feature

We propose to transform the feature function g in the traditional geodesic active contour formulation into a hybrid feature function by incorporating geometric image features with voxel statistics to help automate and regularize the evolving contours. The minimization of the active contour is thus represented by (4)


where for gray-scale intensity MR images, P (I|Φ) represents the probability distribution function of a mixture model (5) from which voxel statistics are drawn, assuming that all voxels are identically and independently distributed and the image is to be described with K class labels.


where P (k) represents the prior probability of the class label k and P (Ik) is the conditional density function of the kth class given Φ, the parameter set of the distribution. We employ Gaussian distributions as


which require a parameter set Φ = {μ, σ}, where μ and σ are the mean and standard deviation. This parameter estimation problem for GMM is solved by applying the EM algorithm [51] to the image intensity histogram.

The design of g in (4) utilizes both a geometric term and a statistical term. Geometrically, the presence of strong image gradients indicates significant structural content. As a result, the contour propagation speed slows to a halt. On the other hand, a lack of edge features often indicates the presence of a homogeneous region. Statistically, high voxel probability indicates a high likelihood of the voxel belonging to the class of interest, warranting a fast contour propagation. If the voxel likelihood is reduced, the contour propagation is slowed down accordingly. The contribution of voxel likelihood to the contour propagation exhibits an inverse behavior to that of image gradients. Since both geometric and statistical features are essential to the contour stability, they can be combined into a single hybrid feature function by modeling the aforementioned behavior as


where the first term is the traditional geometric feature as in (3) and the second term models the inverse behavior of voxel likelihood to image gradients using an inverse sigmoid function with magnitudes between −1 and 1. Complementarily, these two components in the new hybrid feature help regularize the evolving contour in both the geometric and statistical sense. The minimization of (4) is then achieved by computing the Euler–Lagrange equation [34]


where κ is the Euclidean curvature, N is the inward unit normal, and C0 is the initial curve or surface, and performing steepest gradient search (9), to deform C toward a minima


where {ψ, c, ε} are the free parameters introduced to govern curvature, propagation, and advection strengths, respectively. With the designed hybrid feature, the algorithm uses only the propagation term. Other terms are shown here for completeness. This curve evolution equation is then embedded in a level set function u and solved for the steady state solution:


The numerical implementation is based on the curve evolution algorithm via level sets [33], which utilizes an upwind piecewise continuous approximation scheme to provide a numerically stable solution in the presence of singularities.

As summarized in Table I, the rationale behind using this new hybrid feature is to enable handling of situations where the image gradient is high (small sigmoid(|[nabla]I|) value) and the posterior probability of voxel is low, in which the voxel is considered to be a significant feature but lies outside of the desired region. We, therefore, aim to steer the contour slowly away from this voxel by assigning a small negative feature value. On the other hand, if the posterior probability is high, this indicates a significant feature within the desired region; therefore, a small positive feature value is assigned. In contrast, if the image gradient is low (large sigmoid(|[nabla]I|) value) and the posterior probability is low, the voxel is considered to be a weak gradient feature that lies outside of the desired region, warranting a large negative feature value such that contour can be quickly steered away from that region. If the posterior probability is high, a homogeneous area in the desired region is indicated, and is rewarded with a large positive feature value for fast contour expansion. In summary, the proposed hybrid feature provides an adaptive active contour propagation based on local information reflecting both geometry and statistical homogeneity.

Effects of Regularizing Contour Propagation Using Geometric and Statistical Features

C. Segmentation of Brain MRI

Based on the proposed active contour model, we develop a fully automated 3-D brain tissue segmentation algorithm for MR images. The segmentation procedures (Fig. 1) are a generalization and extension of earlier work [52]. We first present the proposed algorithm for T1-weighted (T1w) MRI scans, which are most often used for brain tissue segmentation due to the generally high WM and GM contrast and the reduced effects of WM lesions in patients with neurodegenerative diseases. We later extend the proposed method to simultaneously incorporate additional MR sequence data, such as T2-weighted (T2w) and PD-weighted (PDw) images, in addition to T1w.

Fig. 1
Block diagram of the proposed MRI brain segmentation algorithm. Asterisk (*) denotes the proposed novel hybrid geometric–statistical image feature described in Section II-B.

To segment the brain tissues, we first estimate the GMM parameters such that each mixture distribution represents one single class. Based on these estimated distributions, the normalized posterior probability of each voxel is calculated. We derive the hybrid geometric–statistical feature as described above by combining both the voxel statistics and the image gradient information. To initialize the active contour, we first introduce a voxel threshold ts on the posterior probability to further regularize this contour initialization by effectively removing voxels outside of the desired region boundaries due to noise or partial volume artifacts. Based on the thresholded masks, we form the skeleton representation using standard 2-D morphological sequential thinning (3 × 3 kernel) and sequential pruning (3 × 3 kernel) [53] iteratively until no further changes occur. This is done slice-by-slice in 2-D as morphological operations are performed on discrete numbers of voxels, and to do this in 3-D on anisotropic scans would require resampling of the acquired data, which would potentially introduce additional artifacts. The hybrid feature and the contour initialization are determined individually for each tissue to be segmented, and each contour is then propagated independently. After the contours converge, minor overlaps and gaps occurred between the individually evolving contour boundaries are re-assigned to a single label by comparing the individual z-scores, the difference between voxel intensity to the sample mean normalized by the sample standard deviation, of all tissue classes. A brief description of the algorithm is described as follows:

  • Step 1: Given raw MRI scan U.
  • Step 2: Preprocess with intensity inhomogeneity correction, noise filtering, and brain extraction to obtain V.
  • Step 3: Calculate gradient norm |[nabla]I| of pixel intensity I in V.
  • Step 4: For k = 1 to K, tissue class labels:
    1. determine the GMM distribution parameter set by using EM: Φk = {μk, σk};
    2. for each I, calculate normalized probability P (Ik).
  • Step 5: For k = 1 to K, tissue class labels:
    1. derive hybrid geometric–statistical feature: gk (I) = g(|[nabla]I|, P (Ik));
    2. derive initial contour C0,k by thresholding P (Ik) at tS = 0.1 and skeletonizing;
    3. propagate curve Ck from C0,k on gk (I) until convergence to obtain class label image Lk (I).
  • Step 6: For each I, if Lk (I) = 1 for more than one k or Lk (I) = null for all k, assign Lk (I) = 1 for k with the highest z-score.

In step 5, the performance gain by setting ts > 0 is demonstrated in Table II, and we show that for 0.1 ≤ ts ≤ 0.6, the final segmentation accuracy in a simulated brain volume is shown not to be sensitive to the parameter value tested but is necessary for removing extreme outliers before skeletonizing is performed. The parameter is set to ts = 0.1 for all subsequent experiments, both simulated and clinical.

Parameter Sensitivity of a Simulated Brain Volume (9% Noise, 40% Inhomogeneity)

D. Extension to Multiple MR Sequence Data

To further demonstrate the flexibility of the proposed segmentation approach, we extend our method so that information from multiple MR sequences with different contrast properties can be incorporated when the data is available. Assuming registered images, we first replace the geometric feature component in the proposed hybrid active contour feature (7) with the multidimensional vector gradient norm derived from all available data sequences. To derive the statistical active contour feature term, the GMM parameters and the voxel statistics are individually estimated for each contrast modality m given M input modalities. Since our main goal is to improve the ability to differentiate between various tissue labels, we derive an pairwise intensity contrast term Ci,jm as a ratio of absolute intensity difference between the EM estimated means μ of labels i and j over the observed intensity range of [IMaxIMin] (11) and then normalize (12).


This intensity contrast factor is used as an optimized weighting factor to linearly combine individual posterior probability, P(IΦkm) as in (13), for label k in modality m to form a new posterior probability. This is done by using all relevant pairwise contrast terms, defined between label k and all others. This ensures that the MR sequence with greater contrast receives a greater weight with respect to other lower contrast MR modalities used. Equation (13) shows how the final posterior probability is computed as a function of the means and variances of each tissue and modality. In the case where tissues at different modalities exhibit equal variances, more emphasis is placed on the scan that possesses the greater difference between the expected values of the tissues of interests, thereby facilitating separation. In the other case where the differences in the sample means are equal, the differences in variance are taken into account inherently in the individual posterior probability terms. The resulting probability replaces the posterior probability derived from single contrast input in the segmentation procedure.


E. Data Preprocessing

We employed the nonparametric nonuniform intensity normalization (N3) algorithm [54] for intensity inhomogeneity correction using default parameter settings (width of deconvolution kernel = 0.15, number of iterations = 50, sampling factor = 4, characteristic distance = 200 mm, stopping criteria = 0.001). All scans were then noise filtered using an edge-preserving Perona–Malik anisotropic diffusion filter [55] (number of iterations = 4, conductance = 3.0) to further enhance the image signal-to-noise ratio. Brain masks were generated using the provided ground truths consistent with other published methods to facilitate direct comparisons, or alternatively, many methods are also available for this brain extraction task [56]–[58]. For the clinical datasets where the ground truths are not available, the brain masks were generated with the brain extraction tool (BET) [56] using default parameter settings.

III. Results and Discussion

We applied our proposed segmentation to both simulated and real clinical MRI scans, and demonstrated in the following sections: 1) the accuracy of the proposed segmentation method on simulated T1w brain MRIs; 2) the segmentation improvement on multiple MRI sequences; 3) the accuracy of the proposed method on real clinical MRI scans of normal adults; and 4) the qualitative performance of the proposed methods on clinical MRI scans of MS and AD patients.

We first validated our proposed method on 18 simulated T1w BrainWeb [59] MRI images (with 0%/1%/3%/5%/7%/9% noise, 0%/20%/40% inhomogeneity, 181 × 217 × 181 dimension, 1 × 1 × 1 mm3 spacing). We also performed multisequence segmentation based on six T1w/T2w/PDw MRI triplets (with 0%/1%/3%/5%/7%/9% noise, 0% inhomogeneity). Second, 18 real high-resolution clinical T1w MRI scans from the Internet Brain Segmentation Repository (IBSR) [60] (coronally acquired, 256 × 128 × 256 dimension, 0.837 × 0.837 mm2 to 1 × 1 mm2 in-plane spacing, 1.5 mm slice thickness) were also segmented. For both datasets, the “ground truth” is known for comparisons. For the BrainWeb dataset, the ground truth is the phantom atlas used to generate the simulated scans, whereas for the IBSR dataset, the ground truth is the provided expert-guided manual segmentation label for each of the clinical scans. Lastly, from the MS MRI Research Group (MS/MRI), real clinical 1.5 T spoiled gradient (SPGR) MRI scans (axially acquired, 256 × 256 × 120–160 dimension, 0.937 × 0.937 × 1.50 mm3 spacing) were taken at multiple sites. Real clinical 1.5 T magnetization prepared rapid gradient echo (MP-RAGE) MRI scans (sagittally acquired, 256 × 256 × 166 dimension, 0.937 × 0.937 × 1.20 mm3 spacing) were also obtained from the AD Neuroimaging Initiative (ADNI) of the LONI image data archive (IDA) [61] initiated by the National Institute on Aging (NIA), the National Institute of Biomedical Imaging and Bioengineering (NIBIB), the Food and Drug Administration (FDA), private pharmaceutical companies, and nonprofit organizations. These clinical scans were segmented and qualitatively evaluated. For all experiments, the parameters for the level set evolutions are set at {ψ = 0.0, c = 1.0, ε = 0.0} to reinforce propagation and effectively remove boundary attraction and smoothness regularizations. The convergence criteria for the gradient descent optimization is defined as either less than 0.5% root-mean-square change in the level set function or, if not achieved, a maximum of 100 iterations reached.

The average surface distance between the ground truth and the computed segmentation was computed for each test scan by approximate nearest neighbor searching [62]. In addition, the Dice similarity index [63] was also chosen for the quantitative evaluation of the 3-D brain segmentation results to facilitate direct comparisons to other published results


We denote the true positives, true negatives, false positives, and false negatives as T+, T, F+, and F, respectively, between the known ground truth and the segmentation results. We compared our segmentation results with those of the M3DLS method [41].

A. Segmentation Validation Using Simulated Brain MRI

We first validated a three-class (WM, GM, and CSF) segmentation using the proposed method on the simulated T1w brain MRI data. Segmentation was performed using the traditional geometric feature only, the statistical feature only, and the proposed hybrid feature on all 18 datasets with varying noise and intensity inhomogeneity levels. For the edge-only level set evolution, the parameter set {ψ = 2.0, c = 1.0, ε = 4.5} was used to enforce a stronger smoothness constrain; otherwise, contours leaking through weak edges were often observed. For the statistical feature term only, the parameter set {ψ = 0.0, c = 1.0, ε = 0.0} was used, same as the hybrid approach. Qualitative results in Fig. 2 demonstrated very good resemblance between the provided phantom label and the segmentation results based on the hybrid feature. Table III shows that, on average, over all noise and inhomogeneity levels, the proposed method achieved consistently accurate segmentation results for both WM and GM with similarity indexes of 93.71% (σ = 2.11%, average distance = 0.27mm) and 92.59% (σ = 2.40%, average distance = 0.30mm). Segmentation of structures such as CSF by using the proposed approach also achieved considerable (>70%) similarity of 77.75% (σ = 6.15%, average distance = 2.24 mm). The CSF results were not as stable as WM and GM mainly due to the much smaller structural volume, leading to increased sensitivity to estimation errors in the active contour initialization and feature derivation. Nonetheless, the overall segmentation results were on par if not better than other previously published results [16], [27], [64].

Fig. 2
Qualitative segmentation performance of two simulated T1w brain images showing the provided phantom label, raw images and the segmentation results obtained by using the proposed hybrid feature. We show three slices for both the best case (0% noise, 0% ...
Quantitative Evaluation of Segmentation Results Obtained by Using the Geometric Feature, the Statistical Feature, and the Proposed Hybrid Approach on 18 Simulated T1w BrainWeb Scans

To statistically evaluate the differences of segmentation results between the proposed hybrid approach, and the contours based on geometric-only and statistical-only features, we calculated the p-values (p < 0.05 indicates a statistically significant difference in the group means). Compare to results from using only the traditional geometric feature, the proposed hybrid approach achieved significantly higher similarity indexes and reduced surface distance across all scans. On average, the proposed method achieved increased similarity indexes of 5.36% (p = 0.0002) in WM, 7.23% (p < 0.0001) in GM, and 9.30% (p < 0.0001) in CSF segmentation results with reduced surface distance of 3.99 mm (p < 0.0001) in WM, 0.27 mm (p < 0.0001) in GM, and 4.86 mm (p < 0.0001) in CSF. Compare to results from using only the statistical feature, the proposed hybrid approach only achieved slightly better similarity indexes. However, we observed that, on average, the proposed method was able to significantly reduce the average surface distance by 4.44 mm (p < 0.0001) in WM and 3.79 mm (p < 0.0001) in CSF. These results showed that using the geometric term alone was highly sensitive to image artifacts and require contour regularization, and using the statistical term alone caused the contours to deviate from the true edges. Only when incorporating both features we can evolve the contour in the appropriate statistical space while maintaining high geometric relevance at the same time. Processing a single BrainWeb volume takes approximately 55 min (dual 3.20 GHz Xeon PC with 3.25 GB memory) to 75 min (3.60 GHz Pentium4 PC with 2 GB memory), comparable to the processing time required by other conventional techniques.

B. Segmentation Improvement Using Multiple MR Sequences

We next performed a three-class (WM, GM, and CSF) segmentation using the proposed method on the simulated T1w/T2w/PDw brain MRI data. Segmentation was performed on six datasets with varying noise levels and 0% intensity inhomogeneity. Qualitative results in Fig. 3 illustrated very good resemblance between the provided phantom label and the segmentation results. Table IV demonstrates the quantitative segmentation accuracies. On average, over all noise levels, the proposed method achieved consistently accurate segmentation results for WM, GM, and CSF with similarity indexes of 96.24% (σ = 1.13%, average distance = 0.15 mm), 94.14% (σ = 1.24%, average distance = 0.24 mm) and 81.57% (σ = 2.82%, average distance = 1.56 mm), respectively. When compared to the experiment on single simulated T1w brain images, segmentation using multiple MR sequence data provided an average improvement in similarity indexes of 1.29% (p = 0.0479), 0.44% (p = 0.4627), and 3.55% (p = 0.1403) for WM, GM, and CSF, respectively. The segmentation improvements are not statistically significant, which is not surprising given that the synthetic T1w brain images by themselves already have the strong image contrast required to distinguish between the majority of WM and GM tissues. Additional MR sequences such as T2w and PDw in this case, helped improve the overall robustness by achieving a much better balance between the WM, GM, and CSF estimation as observed by the T+ and T performance. Table III (column 1) shows that with only T1w scans as inputs, the differences between the T+ and T were 7.77%, 8.67%, and 24.09%, respectively, for WM, GM, and CSF, whereas with T1w, T2w, and PDw inputs (Table IV), these differences were much reduced to 2.84%, 0.24%, and 6.95%.

Fig. 3
Qualitative segmentation performance of a multiple simulated MR sequence (T1w, T2w, PDw) brain images showing the provided phantom label, raw images (0% noise, 0% inhomogeneity), and the segmentation results obtained by using the proposed approach. White, ...
Quantitative Evaluation of Segmentation Results Obtained by Using the Proposed Method on 6 Simulated T1w/T2w/PDw Scans

C. Segmentation Comparisons Using Clinical MR Scans of Normal Adults

We applied the proposed method to segment 18 clinical IBSR brain images. The images were segmented using a three-class (WM, GM, and CSF) segmentation. The qualitative results are shown in Fig. 4. The tissue labels were then postprocessed due to a known limitation of the provided manual segmentation labels. It has been reported previously [41] that the expert-guided manual segmentation label contains much of the cortical CSF being mislabeled as GM. We have confirmed this with our own observations. As observed from Fig. 4, we note that the original segmentation results matched closely with what can be visually observed from the raw images. However, this observation did not correspond well to the provided expert-guided manual segmentation label due to the existing limitation. If quantitative evaluation was to be performed on this result, both GM and CSF would produce much lower accuracies than what were actually present.

Fig. 4
Qualitative segmentation performance of a real clinical T1w brain images (IBSR #08) showing the raw images, the expert-guided manual segmentation label, and the segmentation results obtained by using the traditional edge feature and the proposed hybrid ...

To enable a valid comparison by retaining only the ventricular CSF for comparison, we employed the following postprocessing scheme:

  • Step 1: Re-assign segmented CSF labels close to (≤5 mm) any background labels as GM.
  • Step 2: Define a rectangular region of interest (30% of the mask width in the superior–inferior orientation and 50% for both left–right and anterior–posterior) around the brain mask centroid. Apply 2-D morphological thinning (3 × 3 kernel) [65] on CSF pixels outside the region of interest in coronal slices, and singly links are removed.
  • Step 3: Identify the largest 3-D connected component as the desired ventricular CSF. Apply 3-D morphological dilation (3 × 3 × 3 kernel) to this final mask, and re-assign all original CSF labels outside as GM.

This postprocessing scheme allowed us to generate results comparable to the provided manual label. We show in Table V that the results obtained by using the proposed method achieved similarity indexes of 87.55% (σ = 2.92%, average distance = 0.54 mm), 93.18 (σ = 0.92%, average distance = 0.38 mm), and 77.39% (σ = 11.15%, average distance = 3.12 mm) for WM, GM, and CSF, respectively, averaged across the 18 images tested. To statistically compare the segmentation performance of our hybrid active contour approach against the region-based active contour M3DLS algorithm, a one-sample t-test was performed to calculate the p-values (standard deviations were not reported by M3DLS [41]). Both methods are automatic and based on 3-D level set implementation with the key difference in the active contour formulation and the optimization functional. The proposed method achieved significantly higher accuracies with an average improvement of 8.55% (p < 0.0001) in WM and 10.18% (p < 0.0001) in GM segmentation similarities. CSF results were not reported with the M3DLS method; however, our proposed approach achieved considerable (>70%) similarity. Furthermore, the segmentation results achieved by M3DLS were conservative with low true positives. This can be mainly attributed to the fact that the optimization functional in M3DLS modeled each class as piecewise constant by taking only the intensity square differences from the means. However, variations within each segmented class are not considered, such that if, intrinsically, WM contains significantly larger intensity variance than GM, voxels belonging to WM would potentially be mislabeled as GM if the voxel intensities are closer to the estimated GM mean. This lead to undersegmentation in WM (undersegmentation in GM is attributed to the known limitation in the manual labels), and M3DLS is required to postprocess the segmentation results by incorporating an additional morphological dilation step [42]. On the other hand, by integrating both geometric and statistical features into an edge-based deformable model, the proposed hybrid approach captures both the image edge geometry and the voxel statistical homogeneity, and achieved higher segmentation accuracies with less additional computation complexity.

Quantitative Evaluation of Segmentation Results on 18 Real Clinical T1w IBSR Scans

D. Segmentation Performance Using Clinical MR Scans of MS and AD Patients

Lastly, we applied the proposed method to segment clinical MRI brain scans of MS and AD patients. The images were segmented using a three-class (normal appearing WM, GM grouped with diseased WM, and CSF) segmentation. The qualitative results are shown in Fig. 5, demonstrating that the proposed approach appears stable on clinical scans. Fig. 5(a) shows a typical clinical MRI scan of a MS patient, whereas Fig. 5(b) illustrates an example scan of a MS patient with enlarged lateral ventricles. Both scans were segmented without any gross misclassifications. In addition, we demonstrated segmentation results from an MRI scan where the brain extraction step inaccurately included part of the eyes in the brain mask in Fig. 5(c). This error in the brain mask did not seem to cause visible problems for the algorithm in the adjacent brain tissue. Furthermore, Fig. 5(d) and (e) shows segmentation results of scans from an AD patient one year apart. We again noticed visibly consistent segmentation in the tissue labels in the presence of tissue atrophy/ventricle enlargement over time. In our T1w test scans, the intensity difference between GM and diseased WM is subtle, and separating these two class types is likely not possible without additional MRI sequences that are more sensitive to WM pathology, such as PDw or T2w, or relying on prior probability maps such as those derived from a training set. We have left these experiments for future work. In its current form, the proposed method can potentially be used for the assessment of disease severity by providing stable and consistent segmentations of CSF and normal appearing WM.

Fig. 5
Qualitative segmentation performance of real clinical T1w brain images (UBC MS/MRI and LONI IDA) showing the raw images and the segmentation results obtained by using the proposed approach. We show two slices each for the three-class segmentation. White, ...

IV. Conclusion

We proposed a 3-D brain MR segmentation method based on deformable models and demonstrated accurate and stable brain tissue segmentation on single as well as multiple MR sequence scans. The main contribution of our work is that we employed a geodesic active contour formulation by integrating both image geometry and voxel statistics into a hybrid geometric–statistical feature, which acts as a stabilizing regularizing function for the extraction of complex anatomical features such as WM, GM, and CSF. We validated our technique first by using both single and multiple simulated brain MRI sequence data. Improved segmentation accuracy and robustness were shown in results from the proposed hybrid approach against those using individual geometric or statistical features only. Furthermore, on real clinical MRI datasets, we also demonstrated improved accuracy over a state-of-the-art approach, the region-based M3DLS. We also demonstrated consistent and robust results when segmenting MRI scans of both MS and AD patients.

Issues identified for possible future work include enhancing the statistical distribution estimation process by using complex intensity distribution estimation methods such as nonparametric and partial volume models, and extending additional segmentation classes, hierarchy or feature cues for segmentation of anomalies such as WM lesions or tumors.


This work was supported in part by the Natural Sciences and Engineering Research Council (NSERC) Parts of the data used in the preparation of this article were obtained from Alzheimer’s Disease Neuroimaging Initiative (ADNI) database. As such, investigators within ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report.

The authors thank the UBC MS/MRI Research Group for data collection and sharing. In addition, ADNI (Principal Investigator: M. Weiner; National Institutes of Health (NIH) Grant U01 AG024904; a complete listing of investigators is available at ADNI was funded by the National Institute on Aging, National Institute of Biomedical Imaging and Bioengineering, Pfizer, Inc., Wyeth Research, Bristol-Myers Squibb, Eli Lilly and Company, GlaxoSmithKline, Merck & Company, Inc., AstraZeneca AB, Novartis Pharmaceuticals Corporation, Alzheimer’s Association, Eisai Global Clinical Development, Elan Corporation plc, Forest Laboratory, and ISOA, with participation from U.S. Food and Drug Administration. ADNI industry partnerships are coordinated through NIH with the Veterans Health Research Institute (NCIRE) as the grantee organization. ADNI study was coordinated by Alzheimer’s Disease Cooperative Study at University of California, San Diego. ADNI data are disseminated by the LONI at University of California, Los Angeles (



Albert Huang (S’05) received the S.M. degree in engineering science from Harvard University, Cambridge, MA, in 2002, and the B.A.Sc. degree in electrical and computer engineering in 2000 from the University of British Columbia (UBC), Vancouver, BC, Canada, where he is currently working toward the Ph.D. degree at the Department of Electrical and Computer Engineering.

From 2000 to 2004, he was an Algorithm Developer at the Harvard Communications/Molecular and Cellular Biology Laboratories, and Sony Taiwan Limited. His current research interests include signal and image processing.

An external file that holds a picture, illustration, etc.
Object name is nihms238467b1.gif

Rafeef Abugharbieh (M’03) received the Bachelor’s degree in electrical engineering from Jordan University, Amman, Jordan, in 1995, the Master’s degree (with distinction) in digital communication systems from Chalmers University, Gothenburg, Sweden, in 1997, and the Doctoral degree from the Department of Signals and Systems, Chalmers University, in 2001.

In 2004, she joined as an Assistant Professor in the Department of Electrical and Computer Engineering, University of British Columbia (UBC), Vancouver, BC, Canada., where she cofounded the Biomedical Signal and Image Computing Laboratory (BiSICL), in 2005, and is currently the Co-Director. Her current research interests include image processing and analysis particularly in medical imaging applications.

Dr. Abugharbieh has been a member of the IEEE Engineering in Medicine and Biology Society (EMBS) since 2003. She is also an Associate Founding Member of the IEEE EMBS Vancouver Section (2004) and an Associate Editor for Image and Vision Computing.

An external file that holds a picture, illustration, etc.
Object name is nihms238467b2.gif

Roger Tam received the Ph.D. degree in computer science from the University of British Columbia (UBC), Vancouver, BC, Canada, in 2004.

He is currently an Assistant Professor in the Department of Radiology, UBC, and is also a member of the UBC MS/MRI Research Group, a pioneering research laboratory in the field of MRI analysis for the study of multiple sclerosis. His current research interests include medical image processing, computational shape modeling, and scientific visualization.

An external file that holds a picture, illustration, etc.
Object name is nihms238467b3.gif

Contributor Information

Albert Huang, Department of Electrical and Computer Engineering, University of British Columbia (UBC), Vancouver, BC V6T 1Z4, Canada.

Rafeef Abugharbieh, Department of Electrical and Computer Engineering, University of British Columbia (UBC), Vancouver, BC V6T 1Z4, Canada.

Roger Tam, Department of Radiology, University of British Columbia (UBC), Vancouver, BC V6T 1Z4, Canada.


1. Beyer MK, Janvin CC, Larsen JP, Aarsland D. An MRI study of patients with Parkinson’s disease with mild cognitive impairment and dementia using voxel based morphometry. J Neurol Neurosurg Psychiatry. 2007 Mar;78(3):254–259. [PMC free article] [PubMed]
2. Grossman M, McMillan C, Moore P, Ding L, Glosser G, Work M, Gee J. What’s in a name: Voxel-based morphometry analysis of MRI and naming difficulty in Alzheimer’s disease, frontotemporal dementia and corticobasal degeneration. Brain. 2004;127(3):628–649. [PubMed]
3. Grant PE. Structural MR imaging. Epilepsia. 2004;45(s4):4–16. [PubMed]
4. Wisco JJ, Kuperberg G, Manoach D, Quinn BT, Busa E, Fischl B, Heckers S, Sorensen AG. Abnormal cortical folding patterns within Broca’s area in schizophrenia: Evidence from structural MRI. Schizophrenia Res. 2007 Aug;94(1–3):317–327. [PMC free article] [PubMed]
5. Paty DW, Li D, Zhao GJ. MRI in Multiple Sclerosis-Implications for Diagnosis and Treatment. 2. UBC MS/MRI Research Lab; Vancouver, BC, Canada, Rep. for Ares-Serono SA: 1999.
6. Miller DH. Biomarkers and surrogate outcomes in neurodegenerative disease: Lessons from multiple sclerosis. NeuroRx. 2004;1:284–294. [PubMed]
7. Traboulsee A, Zhao G, Li DKB. Neuroimaging in multiple sclerosis. Neurol Clin. 2005;23:131–148. [PubMed]
8. Chard DT, Griffin CM, Parker GJM, Kapoor R, Thompson AJ, Miller DH. Brain atrophy in clinically early relapsing-remitting multiple sclerosis. Brain. 2002 Feb;125(2):327–337. [PubMed]
9. De Stefano N, Matthews PM, Filippi M, Agosta F, De Luca M, Bartolozzi ML, Guidi L, Ghezzi A, Montanari E, Cifelli A, Federico A, Smith SM. Evidence of early cortical atrophy in MS. Neurology. 2003;60:1157–1162. [PubMed]
10. Sastre-Garriga J, Ingle GT, Chard DT, Ramio-Torrenta L, Miller DH, Thompson AJ. Grey and white matter atrophy in early clinical stages of primary progressive multiple sclerosis. NeuroImage. 2004 May;22(1):353–359. [PubMed]
11. Anderson VM, Fox NC, Miller DH. Magnetic resonance imaging measures of brain atrophy in multiple sclerosis. J Magn Reson Imag. 2006;23:605–618. [PubMed]
12. Warfield SK, Zou KH, Wells WM. Simultaneous truth and performance level estimation (STAPLE): An algorithm for the validation of image segmentation. IEEE Trans Med Imag. 2004 Jul;23(7):903–921. [PMC free article] [PubMed]
13. Pham DL, Xu C, Prince JL. Current methods in medical image segmentation. Annu Rev Biomed Eng. 2000 Aug;2:315–337. [PubMed]
14. Liew AWC, Hong Y. Current methods in automatic tissue segmentation of 3D magnetic resonance brain images. Current Med Imag Rev. 2006 Feb;2(1):91–103.
15. Ashburner J, Friston KJ. Voxel-based morphometry—The methods. NeuroImage. 2000;11(6):805–821. [PubMed]
16. Van Leemput K, Maes F, Vandermeulen D, Suetens P. Automated model-based tissue classification of MR images of the brain. IEEE Trans Med Imag. 1999 Oct;18(10):897–908. [PubMed]
17. Pohl KM, Wells WM, Guimond A, Kasai K, Shenton ME, Kikinis R, Grimson WEL, Warfield SK. Incorporating non-rigid registration into expectation maximization algorithm to segment MR images. Proc MICCAI. 2002:564–571.
18. Bookstein FL. Voxel-based morphometry should not be used with imperfectly registered images. NeuroImage. 2001;14:1454–1462. [PubMed]
19. Crum WR, Griffin LD, Hill DLG, Hawkes DJ. Zen and the art of medical image registration: Correspondence, homology, and quality. NeuroImage. 2003;20:1425–1437. [PubMed]
20. Ashburner J, Friston KJ. Why voxel-based morphometry should be used. NeuroImage. 2001;14:1238–1243. [PubMed]
21. Pham DL, Prince JL. Adaptive fuzzy segmentation of magnetic resonance images. IEEE Trans Med Imag. 1999 Sep;18(9):737–752. [PubMed]
22. Wells MW, Grimson WEL, Jolesz FA. Adaptive segmentation of MRI data. IEEE Trans Med Imag. 1996 Aug;15(4):429–442. [PubMed]
23. Rajapakse J, Kruggel F. Segmentation of MR images with intensity inhomogeneities. Image Vis Comput. 1998;16(3):165–180.
24. Burkhardt J. Diploma thesis. Tech. Univ. Munchen; München, Germany: 2003. A Markov chain Monte Carlo algorithm for the segmentation of T1-MR-images of the head.
25. Ruf A, Greenspan H, Goldberger J. Tissue classification of noisy MR brain images using constrained GMM. Proc MICCAI. 2005:790–797. [PubMed]
26. Zhou Y, Bai J. Atlas-based fuzzy connectedness segmentation and intensity nonuniformity correction applied to brain MRI. IEEE Trans Biomed Eng. 2007 Jan;54(1):122–129. [PubMed]
27. Bazin PL, Pham DL. Topology preserving tissue classification of magnetic resonance brain images. IEEE Trans Med Imag. 2007 Apr;26(4):487–496. [PubMed]
28. Liao L, Lin T. MR brain image segmentation based on kernelized fuzzy clustering using fuzzy Gibbs random field model. Proc IEEE ICME. 2007:529–535.
29. Zhang Y, Brady M, Smith S. Segmentation of brain MR images through a hidden Markov random field model and the expectation–maximization algorithm. IEEE Trans Med Imag. 2001 Jan;20(1):45–57. [PubMed]
30. Terzopoulos D, Witkin A, Kass M. Constraints on deformable models: Recovering 3D shape and nonrigid motion. Artif Intell. 1988;36(1):91–123.
31. Suri JS. Leaking prevention in fast level sets using fuzzy models: An application in MR brain. Proc IEEE EMBS. 2000:220–225.
32. Yushkevich PA, Piven J, Hazlett HC, Smith RG, Ho S, Gee JC, Gerig G. User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability. NeuroImage. 2006 Jul;31(3):1116–1128. [PubMed]
33. Osher S, Sethian J. Fronts propagating with curvature-dependent speed: algorithms based on the Hamilton–Jacobi formulation. J Comput Phys. 1988;79:12–49.
34. Caselles V, Kimmel R, Sapiro G. Geodesic active contours. Int J Comput Vis. 1997;22(1):61–79.
35. Zeng X, Staib LH, Schultz RT, Duncan JS. Segmentation and measurement of the cortex from 3-D MR images using coupled-surfaces propagation. IEEE Trans Med Imag. 1999 Oct;18(10):927–937. [PubMed]
36. Goldenberg R, Kimmel R, Rivlin E, Rudsky M. Cortex segmentation—A fast variational geometric approaches. Proc. IEEE Workshop Variational Level SetMethods Comput. Vis; Vancouver, BC, Canada. Jul. 2001; pp. 127–133.
37. Ciofolo C, Barillot C, Hellier P. Combining fuzzy logic and level set methods for 3D MRI brain segmentation. Proc IEEE ISBI. 2004:161–164.
38. Ciofolo C, Barillot C. Brain segmentation with competitive level sets and fuzzy control. Inform Process Med Imag. 2005;3565:333–344. [PubMed]
39. Chan TF, Vese LA. Active contours without edges. IEEE Trans Image Process. 2001 Feb;10(2):266–277. [PubMed]
40. Tsai A, Yezzi A, Jr, Willsky AS. Curve evolution implementation of the Mumford–Shah functional for image segmentation, denoising, interpolation, and magnification. IEEE Trans Image Process. 2001 Aug;10(8):1169–1186. [PubMed]
41. Angelini E, Song T, Mensh BD, Laine A. Segmentation and quantitative evaluation of brain MRI data with a multi-phase three-dimensional implicit deformable model. Proc SPIE (Med Imag) 2004:526–537.
42. Angelini E, Song T, Laine A. Homogeneity measures for multiphase level set segmentation of brain MRI. Proc IEEE ISBI. 2006:746–749.
43. Angelini E, Song T, Mensh BD, Laine A. Brain MRI segmentation with multiphase minimal partitioning: A comparative study. Int J Biomed Imag. 2007;2007:10526-1–10526-15. [PMC free article] [PubMed]
44. Vese LA, Chan TF. A multiphase level set framework for image segmentation using the Mumford and Shah model. Int J Comput Vis. 2002;50(3):271–293.
45. Mumford D, Shah J. Optimal approximation by piecewise smooth functions and associated variational problems. Commun Pure Appl Math. 1989;42:577–685.
46. Jeon M, Alexander M, Pizzi N, Pedrycz W. Unsupervised hierarchical multi-scale image segmentation: level set, wavelet and additive splitting operator. Pattern Recognit Lett. 2005;26:1461–1469.
47. Han X, Xu C, Prince JL. A topology preserving level set method for geometric deformable models. IEEE Trans Pattern Anal Mach Intell. 2003 Jun;25(6):755–768.
48. Pohl KM, Kikinis R, Wells WM. Active mean fields: Solving the mean field approximation in the level set framework. Proc IPMI. 2007:26–37. [PMC free article] [PubMed]
49. Rousson M, Paragios N, Deriche R. Implicit active shape models for 3D segmentation in MR imaging. Proc MICCAI. 2004:209–216.
50. Ibanez L, Schroeder W, Ng L, Cates J. The ITK Software Guide. 2. New York: Kitware, Inc; 2005.
51. Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm. J R Stat Soc Ser B (Methodological) 1977;39(1):1–38.
52. Huang A, Abugharbieh R, Tam R, Traboulsee A. Automatic MRI brain tissue segmentation using a hybrid statistical and geometric model. Proc IEEE ISBI. 2006:394–397.
53. Gonzales RC, Woods RE. Digital Image Processing. 2. Englewood Cliffs, NJ: Prentice-Hall; 2002. p. 549.
54. Sled JG, Zijdenbos AP, Evans AC. A non-parametric method for automatic correction of intensity non-uniformity in MRI data. IEEE Trans Med Imag. 1998 Feb;17(1):87–97. [PubMed]
55. Perona P, Malik J. Scale-space and edge detection using anisotropic diffusion. IEEE Trans Pattern Anal Mach Intell. 1990 Jul;12(7):629–639.
56. Smith SM. Fast robust automated brain extraction. Human Brain Mapping. 2002 Nov;17(3):143–155. [PubMed]
57. Shan ZY, Yue GH, Liu JZ. Automated histogram-based brain segmentation in t1-weighted three-dimensional magnetic resonance head images. NeuroImage. 2002;17:1587–1598. [PubMed]
58. Huang A, Abugharbieh R, Tam R, Traboulsee A. MRI brain extraction with combined expectation maximization and geodesic active contours. Proc IEEE ISSPIT. 2006:107–111.
59. Cocosco CA, Kollokian V, Kwan RKS, Evans AC. BrainWeb: Online interface to a 3D MRI simulated brain database. NeuroImage. 1997;5(4):425.
60. Worth A. Internet brain segmentation repository. 2008. Jan 20, [Online]. Available:
61. Weiner MW. Alzheimer’s Disease Neuroimaging Initiative. VA Medical Center and University of California; San Francisco: 2008. Nov, [Online]. Available:
62. Arya S, Mount DM, Netanyahu NS, Silverman R, Wu AY. An optimal algorithm for approximate nearest neighbor searching. J ACM. 1998;45:891–923.
63. Dice LR. Measures of the amount of ecologic association between species. Ecology. 1945;26:297–302.
64. Shattuck DW, Sandor-Leahy SR, Schaper KA, Rottenberg DA, Leahy RM. Magnetic resonance image tissue classification using a partial volume model. NeuroImage. 2001;13(5):856–876. [PubMed]
65. Lam L, Lee SW, Suen CY. Thinning methodologies—A comprehensive survey. IEEE Trans Pattern Anal Mach Intell. 1992 Sep;14(9):869–885.