Automatic segmentation of brain MR images using an adaptive balloon snake model with fuzzy classification

Abstract

Skull-stripping in magnetic resonance (MR) images is one of the most important preprocessing steps in medical image analysis. We propose a hybrid skull-stripping algorithm based on an adaptive balloon snake (ABS) model. The proposed framework consists of two phases: first, the fuzzy possibilistic c-means (FPCM) is used for pixel clustering, which provides a labeled image associated with a clean and clear brain boundary. At the second stage, a contour is initialized outside the brain surface based on the FPCM result and evolves under the guidance of an adaptive balloon snake model. The model is designed to drive the contour in the inward normal direction to capture the brain boundary. The entire volume is segmented from the center slice toward both ends slice by slice. Our ABS algorithm was applied to numerous brain MR image data sets and compared with several state-of-the-art methods. Four similarity metrics were used to evaluate the performance of the proposed technique. Experimental results indicated that our method produced accurate segmentation results with higher conformity scores. The effectiveness of the ABS algorithm makes it a promising and potential tool in a wide variety of skull-stripping applications and studies.

Appendices

Appendix 1: Cluster analysis

Cluster analysis aims to partition a set of objects into several groups (clusters) by assigning different labels to each individual object. In general, groups are identified according to some specific senses such as connectivity-based, distribution-based, and centroid-based classification. A clustering algorithm can be generally defined as clustering an unlabeled data set $D=\left\{d_1,d_2,\ldots,d_n\right\}$ into c subgroups, where 1 < c < n is the number of the clusters usually predefined by the user. Vectors in the same partition are assigned with an identical label. The c-partition of D can be arrayed as a matrix U _c × n  = [u _ik ] with size c × n. The value of u _ik varies between different clustering algorithms. Moreover, a set of centroids or prototypes $V=\left\{v_1,v_2,\ldots,v_c\right\}$ is further defined for representing each cluster in centroid-based models.

The fuzzy c-means (FCM), proposed by Bezdek [2], is a centroid-based model that incorporates the fuzzy concept in the membership, i.e., the membership of each object d _k is defined as the inversely relative distance of d _k to the centroid V. The membership of the data set can be arrayed as U _FCM = [u _ik ] and is constrained by the following equations:

$$\sum_{i=1}^{c}u_{ik}=1,\quad \forall{k}$$

(17)

and

$$\sum_{k=1}^{n}u_{ik}>0,\quad \forall{i}.$$

(18)

Equation (17) indicates that the membership of an object d _k in each cluster is a normalized value such that the membership of d _k sums to 1 while Eq. (18) implies that there exists at least one object associated with positive membership in each centroid.

Since the membership is the inverse function of the object distance to the centroid, one problem of the FCM is that for points that are naturally identified as outliers but equidistant from two prototypes, the same membership is given to these points. The calculated values of membership thus provide unrepresentative associations between objects and prototypes. Such drawback is resulted from Eq. (17) in that the membership is normalized and unable to present the actual spatial relationship.

Subsequently, Krishnapuram and Keller [15] proposed the possibilistic c-means (PCM), in an attempt to relieve the problem in the FCM. The PCM model replaces the column sum constraint with a looser form as:

$$0<\sum_{i=1}^{c}u_{ik}<c \quad \forall k\,\exists i\ni u_{ik}>0,$$

(19)

where each element u _ik is between 0 and 1. The value of u _ik is interpreted as typicality instead of membership of d _k relative to the centroid v _i . It is also recommended to interpret each row of U as a possibility distribution over D. The PCM model somewhat overcomes the drawbacks of the FCM while it sometimes suffers from coincident cluster problems, which refer to incorrect assignment of an object that naturally belongs to another [1].

Fuzzy possibilistic c-means (FPCM), as the name implies, combines the features of the fuzzy and possibilistic c-means and address the problems of these two models. The FPCM finds out the optimal solution of classification by minimizing the following objective function,

$$\min_{{\bf U},{\bf T},{\bf V}}J_{m,\eta}({\bf U},{\bf T},{\bf V};{\bf X}) =\sum_{i=1}^{c}\sum_{k=1}^{n}(u_{ik}^{m}+t_{ik}^{\eta})||x_k-v_i||^{2},$$

(20)

where U denotes the relative typicality (membership) described in the FCM, T represents the absolute typicality in the PCM, and V is a vector of element v _i representing the center belonging to cluster i. The symbols m and η are weighting exponents with m > 1 and η > 1. The constant c is the number of clusters and k the number of data points. Details of the minimization process and proof are given in [20]. Herein, we briefly describe the necessary conditions for minimizing the objective function J _m,η as follows:

$$u_{ik}=\left[\sum_{j=1}^{c}{\left( \frac{||x_{k}-v_{i}||}{||x_{k}-v_{j}||}\right)}^{\frac{2}{m-2}}\right] ^{-1},\quad \forall i,\,k$$

(21)

and

$$t_{ik}=\left[\sum_{j=1}^{n}{\left(\frac{||x_{k}-v_{i}||}{||x_{j}-v_{i}||}\right)} ^{\frac{2}{\eta-2}}\right]^{-1},\quad \forall i,\,k.$$

(22)

Based on these conditions, vector V is updated using

$$v_i = \frac{\sum_{k=1}^{n}(u_{ik}^{m} +t_{ik}^{\eta})x_{k}}{\sum_{k=1}^{n}(u_{ik}^{m}+t_{ik}^{\eta})},\quad \forall i.$$

(23)

Appendix 2: Balloon snake models

Snakes (also known as parametric active contours), proposed by Kass et al. [14], have been widely applied in image segmentation and object tracking. A snake can be defined as a set of ordered points or snaxels v(s) = [x(s), y(s)], usually generated counter-clockwise. The parameter $s\in [0,1]$ is a normalized arc length starting from the first snaxel. The deformation of each snaxel is governed by both internal and external forces. The internal force is related to the stretching ability (or tension) and smoothness of the curve, which shrinks to a tiny circle when the driving forces are only internal. While the internal force is independent from the image data, the external force is related to salient features such as terminations and edges in images. The balance of internal and external forces drives the snake curve moving toward the object boundary while simultaneously maintaining the tension and stiffness.

Subsequently, Cohen et al. [8] embedded a balloon force into the traditional snake models to solve the problem of limited moving distances. This enhanced version of snakes simulates the action of balloons including deflation and inflation along the normal direction. An overview of the balloon snake model is shown in the following equation:

$$E_{\rm snake}=\int\limits_0^1 E_{\rm int}(v(s))+E_{\rm ext}(v(s))+E_{\rm bal}(v(s)) \hbox{d}s,$$

(24)

where E _int,   E _ext and E _bal represent the internal, external, and balloon energy, respectively.

The internal energy in Eq. (24) is defined as:

$$E_{\text{int}} = \alpha(s)|v_{s}|^2+\beta(s)|v_{ss}|^2,$$

(25)

where α(s) and β(s) are weighting functions, and subscripts are used to indicate derivatives. The first order derivative with respect to s controls the distance between adjacent snaxels. During energy minimization, the second order term makes the contour resist bending. Consequently, the snake contour tends to collapse in the absence of other constraints or forces. The relative strength of tension and stiffness can be adjusted by controlling the values of α(s) and β(s).

Alternatively, the external energy in Eq. (24) is defined as

$$E_{\text{ext}}=-\gamma \mid \nabla G_{\sigma}(x,y) \times I(x,y)\mid^2,$$

(26)

where G _σ denotes the Gaussian filter with σ controlling the spatial extent of the local minima of the convolution kernel, ∇ is the gradient operator, I is the image intensity, and the notation $|\cdot |$ represents norm. The parameter γ is a weighting function for controlling the magnitude of the external energy. Accordingly, regions with salient features have relatively smaller external energy while homogeneous regions are associated with higher external energy.

By embedding a normal force into each snaxel, the balloon energy increases the moving distance using:

$$E_{\text{bal}}= \kappa n(s).$$

(27)

where κ is a weighting function and n(s) is the normal vector that is further resolved into x and y components:

$$n_x(s)=\frac{y_{s+1}-y_{s-1}}{\sqrt{(x_{s+1}-x_{s-1} )^2+(y_{s+1}-y_{s-1})^2}}$$

(28)

and

$$n_y(s)=\frac{-(x_{s+1}-x_{s-1})}{\sqrt{(x_{s+1}-x_{s-1} )^2+(y_{s+1}-y_{s-1})^2}}$$

(29)

where x _s and y _s in pair represent the coordinates of the sth snaxel.

According to the theory of the fundamental Euler-Lagrange differential equation and minimizing E _snake, the snake evolution is achieved when the following equation is satisfied:

$$F_{\rm int}+F_{\rm ext}+F_{\rm bal}=0$$

(30)

where

$$F_{{\rm int}}=\alpha (s)v_{ss}-\beta (s)v_{ssss},$$

(31)

$$F_{{\rm ext}}=-\nabla E_{{\rm ext}},$$

(32)

and

$$F_{{\rm bal}}=\nabla E_{{\rm bal}}.$$

(33)