Brain Extraction Using Active Contour Neighborhood-Based Graph Cuts Model

: The extraction of brain tissue from brain MRI images is an important pre-procedure for the neuroimaging analyses. The brain is bilaterally symmetric both in coronal plane and transverse plane, but is usually asymmetric in sagittal plane. To address the over-smoothness, boundary leakage, local convergence and asymmetry problems in many popular methods, we developed a brain extraction method using an active contour neighborhood-based graph cuts model. The method deﬁned a new asymmetric assignment of edge weights in graph cuts for brain MRI images. The new graph cuts model was performed iteratively in the neighborhood of brain boundary named the active contour neighborhood (ACN), and was e ﬀ ective to eliminate boundary leakage and avoid local convergence. The method was compared with other popular methods on the Internet Brain Segmentation Repository (IBSR) and OASIS data sets. In testing cross IBSR data set (18 scans with 1.5 mm thickness), IBSR data set (20 scans with 3.1 mm thickness) and OASIS data set (77 scans with 1 mm thickness), the mean Dice similarity coe ﬃ cients obtained by the proposed method were 0.957 ± 0.013, 0.960 ± 0.009 and 0.936 ± 0.018 respectively. The result obtained by the proposed method is very similar with manual segmentation and achieved the best mean Dice similarity coe ﬃ cient on IBSR data. Our experiments indicate that the proposed method can provide competitively accurate results and may obtain brain tissues with sharp brain boundary from brain MRI images.

The region-based methods first use thresholding or clustering techniques to divide brain MRI image into several regions by considering that the voxels in the same tissue have similar intensities, and then the brain region can be extracted from these regions by morphological operations or region merging. These thresholding or clustering techniques frequently used in brain extraction are Gaussian mixture model [8], intensity thresholding [9] and watershed algorithm [10]. The parameters in the region-based methods have remarkable effects on the brain extraction and need to be determined properly [12].
The boundary-based methods partition the brain MRI image into internal part (brain) and external part (non-brain) by detecting the boundary between brain and non-brain tissues. Smith [2] proposed Symmetry 2020, 12, 559 2 of 14 a Brain Extraction Tool (BET) by using smoothing and pushing forces to push a tessellated mesh to the brain boundary. Shattuck et al. [6] proposed the Brain Surface Extraction (BSE) method by using a Marr-Hildreth edge detector to separate brain and non-brain tissues. The boundary-based methods [11,12] used active contour model that has made great progress in medical image segmentation during recent years.
The hybrid methods try to use a two steps strategy to improve the extraction result. In the first step, a rough brain region or boundary is obtained to be the initial brain region or boundary in the next step. Next, the brain contour or region is refined to get a more accurate result since it is close to the brain boundary at the first stage. Huang et al. [13] employed the expectation maximization algorithm on a mixture of Gaussian models to determine the initial brain contour for the geodesic active contour evolution. Ségonne et al. [14] proposed a hybrid watershed algorithm (HWA) applying the watershed algorithm to get an initial brain volume for deformable brain mesh. Sadananthan et al. [15] used graph cuts based image segmentation to refine the preliminary binary mask generated by intensity thresholding. Jiang et al. [16] used BET to generate the initial brain boundary, and then refined the brain boundary using a hybrid level set model.
To improve the robustness, some hybrid-based methods warp the brain volume to an atlas using registration techniques before brain extraction, then use the parameter learning techniques such as meta-algorithm, random forest and neural networks to get the proper initial region or parameters to perform brain extraction. Wang et al. [17] warped an atlas to the brain volume to obtain the initial brain mask, and then refined the mask with a deformable-surface similar with BET. Iglesias et al. [18] proposed a robust, learning-based brain extraction system (ROBEX). ROBEX used a trained random forest classifier to detect the brain boundary. Then the contour was refined to obtain the final brain tissues using graph cuts. Eskildsen et al. [19] proposed an atlas-based method with nonlocal segmentation technique, without the need of time-consuming non-rigid registrations. Huang et al. [20] used a trained Locally Linear Representation-based Classification (LLRC) method for brain extraction. By using registration techniques, the atlas-based methods can address the problem of variability in anatomy and modality and produce more robust results. However, the registration in atlas-based methods is usually highly time-consuming and if the registration fails, bad result may be obtained.
In recent years, the deep learning techniques were used for brain extraction. Kleesiek et al. [21] used a 3D convolutional neural network for brain extraction. Saleh et al. [22] applied U-Net for brain extraction. However, the output of these learning-based brain extraction method tends to be over-smoothed if the training data are obtained by automatic brain extraction method under the checking by human experts or by weak bounding box labeling [23]. In other words, if we want the learning-based brain extraction method to output brain tissue with sharp boundary, the training data should be quite precise and be obtained by manual. For example, Hwang et al. [24] recently applied 3D U-Net for brain extraction and achieved high extraction accuracy by using a fine training data.
Although the above methods have greatly improved the accuracy and robustness for brain extraction, these methods can not completely substitute for manual method for the appearances of over-smoothness, leakage through a weak boundary and missing brain tissues caused by local convergence. Unfortunately, it is a conflicting problem to improve all of them. Brain extraction is a compromised problem where a semi-global understanding of the image is required as well as a local understanding [2]. For example, it is effective that we use local region features to obtain sharp brain boundary and eliminate the leakage, but it is easy to lead to local convergence at the edges between the white matter and gray matter if the initial region is far from the true brain boundary. Trying to address these problems, we proposed a new brain extraction method from T1-weighted MRI volume. The main contributions of our work can be summarized as follows: 1. We defined a new asymmetric assignment of edge weights in graph cuts for brain MRI image to obtain brain boundary. 2. We performed the new graph cuts model iteratively in the neighborhood of brain boundary named the active contour neighborhood (ACN), and the new model was effective to eliminate boundary leakage and avoid local convergence.
3. By the asymmetric function of edge weights, we reduced the edge weights in the region between the ACN and the region in the current boundary to obtain a sharp brain boundary.
The remainder of this paper is organized as follows. In Section 2, the details of the designed brain extraction method including the ACN and graph cuts are introduced. In Section 3, experiments and results are presented. Finally, we give a detailed discussion on the analysis of our experiments in Section 4.

Data Sets
To measure the extraction accuracy of our method, we used the following three data sets: 1. Data set 1 (IBSR18): 18 normal T1-weighted MR Image data sets with expert segmentations from IBSR. Each volume has around 128 coronal slices, with 256 × 256 pixels per slice and 1.5mm thickness.
2. Data set 2 (IBSR20): 20 normal T1-weighted MR Image data sets from IBSR, with 256 × 256 pixels per coronal slice and 3.1 mm thickness. Obvious intensity inhomogeneity and other significant artifacts present in most of the MR images in this data set. Another challenge of this data set is that the neck and even shoulder areas are included.
3. Data set 3 (OASIS77): 77 T1-weighted MR Image data sets were obtained from OASIS project. Each volume has around 208 coronal slices, with 176 × 176 pixels per slice and 1mm thickness. The brain masks for this set were obtained by an atlas-based brain extraction method automatically, checked by human experts before releasing the data. Although the lack of exactitude makes it unfit for testing the precision of a method, we can use it to test the robustness of a method because it includes scans from a very diverse population with a very wide age range as well as diseased brains [19].
We downloaded IBSR data sets from https://www.nitrc.org/projects/ibsr and OASIS77 data set from http://www.oasis-brains.org/ by agreeing the license agreements to access all the data sets.

Graph Cuts
In our previous work [16] we used BET to pre-process T1-weighted MRI scans to generate a robust initial contour and obtained the active contour neighborhood (ACN) [25] by dilating the contour with a small size. In this paper, we used a modified graph cuts to refine the brain contour in ACN. Thus, the ACN and brain contour can be obtained iteratively. Since the initial contour obtained by BET is close to the real brain boundary, the real brain boundary should be the global minimum within the ACN obtained by the initial contour, and the redefined edge weights are more powerful than GCUT [15] for brain extraction.
The graph cuts approaches [26][27][28][29] model the image as a weighted undirected graph G = (V, E, W), where V, E and W are sets of vertices, edges and edge weights respectively. Let A = (A 1 ,...,A p ,..., A |P| ) be a binary vector whose components A p specify assignments to pixels p in data set P. Each A p can be either "obj" or "bkg" (abbreviations of "object" and "background"). In the graph described by Boykov and Jolly [26], an energy function was defined as below: where R(A) is the regional term and B(A) is the boundary term. R p ("obj") and R p ("bkg") reflect how the intensity of the pixel p fits the histogram of the object (O) and background (B) models respectively. The boundary term B(A) was described by: A well-known Min-Cut/Max-Flow algorithm can be used to minimize the energy function, and the cut will separate the image into deferent regions. Usually, the graph G = (V, E, W) has two special vertices (terminals), namely, the source S and the sink T . Let {p, q} be the two vertices in the graph and N be the neighbors of them. The edge weights between all the vertices in the graph are described in Table 1. Table 1. The edge weights in the graph.

Edge
Weight for

Active Contour Neighborhood-Based Graph Cuts Model
To address the over smoothness, boundary leakage and asymmetry problem, we used a piecewise function to define the edge weights in graph cuts, which was more powerful than GCUT. To avoid local convergence, we performed graph cuts in ACN.

Description of ACN Model (ACNM)
The elements of the pipeline for the ACN model are shown in Figure 1. The first step is to obtain a rough boundary by a modified BET method for 2D images [30] from a round initial contour. Although the rough boundary is too smooth, it is close to the brain tissue and is still fit well for providing a good initial contour for brain extraction, due to the good robustness of BET. Then, the ACN is obtained by dilating the rough boundary. In the ACN, a graph cuts model is defined and performed to obtain a new boundary. Finally, a more accurate and sharper brain boundary is obtained by iteratively replacing the new boundary and the ACN. The maximum iteration number in Figure 1 was 10 in this paper.
Steps of the proposed brain extraction method.

Edge weights assignment in ACNM
In the defined graph cuts model for brain extraction, the edge weights between all the vertices in the graph are described in Table 2. The brain tissues and non-brain tissues are considered as the Object () and Background () respectively in our work, thus the ACN in Table 2 exactly corresponds to the area , Table 1 (see in Figure 2). R("brain") and R("nonbrain") are the likelihood items reflecting how the vertex is fit for brain or non-brain tissues respectively. For ACN, we firstly classified the ACN into "brain" and "nonbrain" regions using Fuzzy C-Means clustering method (FCM) [31,32], then obtained the their likelihood items from Equation (3) and Equation (4).
Recalling that  should include most of the non-brain tissues and inversely in , we simply set the edge weights between the terminal ( and ) and vertex in  and  as the minimal or maximal value of R("brain") and R("nonbrain") in ACN respectively. Table 2. The edge weight in the graph used in ACN.

Edge
Weight for Steps of the proposed brain extraction method.

Edge Weights Assignment in ACNM
In the defined graph cuts model for brain extraction, the edge weights between all the vertices in the graph are described in Table 2. The brain tissues and non-brain tissues are considered as the Object (O) and Background (B) respectively in our work, thus the ACN in Table 2 exactly corresponds to the area p ∈ P, p O ∪ B in Table 1 (see in Figure 2). R("brain") and R("nonbrain") are the likelihood items reflecting how the vertex is fit for brain or non-brain tissues respectively. For ACN, we firstly classified the ACN into "brain" and "nonbrain" regions using Fuzzy C-Means clustering method (FCM) [31,32], then obtained the their likelihood items from Equation (3) and Equation (4). Recalling that B should include most of the non-brain tissues and inversely in O, we simply set the edge weights between the terminal (S and T ) and vertex in O and B as the minimal or maximal value of R("brain") and R("nonbrain") in ACN respectively. Table 2. The edge weight in the graph used in ACN.

Edge
Weight for Symmetry 2020, 12, x FOR PEER REVIEW 5 of 14 Figure 1.
Steps of the proposed brain extraction method.

Edge weights assignment in ACNM
In the defined graph cuts model for brain extraction, the edge weights between all the vertices in the graph are described in Table 2. The brain tissues and non-brain tissues are considered as the Object () and Background () respectively in our work, thus the ACN in Table 2 exactly corresponds to the area , Table 1 (see in Figure 2). R("brain") and R("nonbrain") are the likelihood items reflecting how the vertex is fit for brain or non-brain tissues respectively. For ACN, we firstly classified the ACN into "brain" and "nonbrain" regions using Fuzzy C-Means clustering method (FCM) [31,32], then obtained the their likelihood items from Equation (3) and Equation (4).
Recalling that  should include most of the non-brain tissues and inversely in , we simply set the edge weights between the terminal ( and ) and vertex in  and  as the minimal or maximal value of R("brain") and R("nonbrain") in ACN respectively.  When setting the edge weights between the vertex and its neighbors (N) B{p,q}, we used the following equations: B p, q = D(p, q) · G(p, q) · R(p, q) where I max = max(t m , I(0), I(1), · · · , I(d)) Compared with the B{p,q} used in GCUT [15], we redefined the D(p,q) and G(p,q) in Equations (9)-(12) and added a new regional item R(p,q) in Equations (13)- (15).
In (10) d pb is the minimal distance between p and the current brain boundary, and W is the width of the image. This assignment increases the weight of edges locating deep within the foreground region, making cuts here less likely. k D is a parameter to control the smoothness. If over-smoothness occurs, a small k D can be selected to reduce the edge weights in the region between the ACN and R b (the region in the current boundary), and then the current brain boundary will move inside to touch the true brain boundary. The new G(p,q) in Equation (11) was defined by a piecewise function (12), which depends on the grayscale value of p, the mean grayscale value µ of the ACN and the maximal value I m of the ACN. The plot corresponding to Equation (12) is shown in Figure 3.  The regional term R(p,q) in Equations (13)−(15) is used to eliminate boundary leakage trough the edge between the brain and eyeball. R(p,q) depends on the regional maximal value max I which was firstly used in BET and then used by Zhuang et al. [11] and Liu et al. [12]. max I is searched along a line pointing inward from the current vertex within a distance d. tm is the median intensity of the brain tissue on each slice, which is approximated from the pixels within the initial brain region. max I increases the edge weights on the eyeball to avoid the brain boundary leaking through the eyeball tissue. In Figure 3, µ = 0.5, I m = 0.9, k 0 = 6 and k 1 = 0.55k 0 , they control the range of the grayscale values that could be considered inside the brain tissues. Generally, the true brain boundary should locate rather on the vertices with grayscale values close to µ than on the vertices with lower or higher grayscale values. So, wherever I(p) equals to µ, we set I G (p) = 1 to make the voxels whose grayscale values are close to µ touched by the active contour. Wherever I(p) µ, we set I G (p) > 1 and the active contour will move away from the voxels whose grayscale values are far from µ. Considering the true Symmetry 2020, 12, 559 7 of 14 brain boundary more likely locates on the lower intensity side of the vertices whose grayscale values are close to µ, we set k 1 = 0.55k 0 . Thus, Equation (12) is a piecewise function, which is also well fit for the asymmetric distribution of intensity in the brain images.

Evaluation Metrics
The regional term R(p,q) in Equations (13)-(15) is used to eliminate boundary leakage trough the edge between the brain and eyeball. R(p,q) depends on the regional maximal value I max which was firstly used in BET and then used by Zhuang et al. [11] and Liu et al. [12]. I max is searched along a line pointing inward from the current vertex within a distance d. t m is the median intensity of the brain tissue on each slice, which is approximated from the pixels within the initial brain region. I max increases the edge weights on the eyeball to avoid the brain boundary leaking through the eyeball tissue.

Comparison to Other Methods
Because BET, BSE, GCUT and ROBEX are all free and publicly available, a comparison to them was performed on the chosen data sets. The software of our work was programmed with matlab2016 and we tested it on a computer with 8GB RAM and Intel i7-4790 CPU. The software cost about 2 min to process one MRI volume on this computer. In this software, only two parameters k D and k 0 in Equation (10) and Equation (12) were needed to be set by the user for different volumes. k D is positively associated with the smoothness degree of brain boundary. k 0 is negatively associated with the average intensity of brain tissue. In testing, we chose k D = 0.80 and k 0 = 5.5 for IBSR18, k D = 0.83 and k 0 = 1.2 for IBSR20, k D = 1 and k 0 = 6 for OASIS77 to get the best results. The extremely low k 0 used in IBSR20 was to deal with the extremely high value and inhomogeneity of intensities in brain tissue. Before performing our method, the modified BET was done on the coronal middle slice once to obtain the initial contour. The initial contours of the other slices in the 3D volume were obtained by the resultant brain boundary of previous slice in a slice by slice manner, which was described in detail in our previous work [16]. The parameters of the BET were almost the same for all the data sets, which are quite easy to be set according to our previous work [30] for the robustness of BET. When estimating the centroid and the equivalent radius of the brain in BET, extra biases were added to them for IBSR20 because the images in IBSR20 include a lot of tissues in neck. For GCUT and ROBEX, we used the default parameters setting in their software. For BET and BSE, we used the different parameters on different volumes in the three data sets to obtain the best results by processing each volume for several times.
Tables 3-5 display the means and standard deviations of the metrics from each method on the three data sets. The P * -values listed in Tables 3-5 are the adjusted P-values of paired t-tests with the Bonferroni-Holm correction to show the statistical difference between the compared methods and the proposed method. Most of the P*-values are below 0.05 except those of FN RATE and DS for BSE on IBSR18 and IBSR20 data sets, which means the performance of the proposed method was different with the compared methods. Figures 4 and 5 display sample outputs from each method on both IBSR data sets and mark false negative and false positive voxels with different colors to show the extraction errors. To clearly and correctly show the visual extraction errors on all test images for each method, we first warped all the results and the ground truth to the atlas [19], thus every 3D brain extracted  Figures 6 and 7. The brighter the hot color in Figures 6 and 7 is, the bigger the error on the corresponding area is. Because the ground truth of OASIS77 (see in Figure 8) is not as precise as those of IBSR18 and IBSR20 data sets, we do not show the errors for OASIS77 in Figures 4-7, but directly display the sample outputs from each method on both IBSR data sets and OASIS77 data set in Figure 8 to show the differences among them.
BET shows good results on IBSR18 and OASIS77 data sets. However, it performs poorly on IBSR20 data set, because of the included neck and shoulder areas featured with the obvious intensity inhomogeneity. BET then shows both high FP Rate and high FN Rate in Table 4 caused by a bad estimation of brain center. The obvious intensity inhomogeneity also makes it difficult for BET to obtain a good result with the same parameter for all slices of the MRI volume, leading to low DS and JS.
BSE performs well on most of the MRI images in IBSR18, IBSR20 and OASIS77 data sets, however it brought very bad results on some of MRI images whatever we changed the parameters (Some of the tissues in metencephalon and cerebellum were missed by BSE in Figures 5 and 8 due to the obvious intensity inhomogeneity in the image). GCUT shows poor results on both IBSR data sets with the highest FP Rate . A very high FP Rate shows that it tends to maintain non-brain tissue by GCUT. ROBEX performs equally on all data sets, as it is robust across different data. However, the DS and JS for ROBEX are not high in Tables 3 and 4. By the proposed ACNM method, the average DS is 0.957 for IBSR18, 0.960 for IBSR20, and 0.936 for OASIS77. The proposed ACNM method shows the highest extraction accuracy on both IBSR data sets. Both FP Rate and FN Rate are lower than 5%, indicating that a balance of FN Rate and FP Rate is achieved by the proposed method on both IBSR data sets. From the coefficients in Table 5, it seems ACNM does not perform as well as ROBEX and GCUT on OASIS77 data set. It is important to note that the ground truth of OASIS77 data set is not as precise as IBSR data set and has the tendency to over cover the brain, which is similar to the result from ROBEX and GCUT. So, it is reasonable that ROBEX and GCUT obtain better evaluation on OASIS77 than ACNM does. In Figure 8, it is clear that ACNM does the best on OASIS77. Table 3. Comparison of our method with BET, BSE, GCUT and ROBEX using IBSR18 Data Set.

Method DS Mean(SD) JS Mean(SD) FP RATE (%)Mean(SD) FN RATE (%)Mean(SD)
BET P * -value           We also compared the proposed method with the deep learning-based method (CNN) proposed by Kleesiek et al. [21] in literature due to no available soft for it. In [21], only the average DS, SE and SP for all the images in IBSR and OASIS data sets were listed, so we listed them in Table 6. The DS, SE and SP of ACNM are very close to those of CNN, so the proposed method almost reached the deep learning-based method without using any learning and registration techniques. We also compared the proposed method with the deep learning-based method (CNN) proposed by Kleesiek et al. [21] in literature due to no available soft for it. In [21], only the average DS, SE and SP for all the images in IBSR and OASIS data sets were listed, so we listed them in Table 6. The DS, SE and SP of ACNM are very close to those of CNN, so the proposed method almost reached the deep learning-based method without using any learning and registration techniques.

Discussion
Among these methods evaluated in this paper, BET, BSE and GCUT are very sensitive to parameters. If the parameters are well set, BSE can obtain a highly accurate result. However, BSE sometimes failed to obtain a good result no matter how the parameters were set. BET performed with high robustness if the brain center was well predicted. The remain of a lot of neck tissues in IBSR20 made BET predict wrong center of the brain and obtain very bad result for IBSR20. Although our proposed method used BET in the initial step, it achieved much better results than BET due to a better way to calculate the brain center. GCUT did not obtain any good results on both IBSR data sets using fixed parameters suggested in the literature. ROBEX did not achieve good result on IBSR18 with quite different ranges of intensity between MRI subjects. The main reason may be that the graph cuts models in the refine step of ROBEX and GCUT are sensitive to parameters and usually lead to over smoothness. The proposed method used BET to initialize the brain contour which is very close to the brain boundary and ensures that the active contour evolves to the true brain boundary by performing a modified graph cuts model in the neighborhood of the contour (ACN). The modified graph cuts model used an asymmetric edge weighting function which is more powerful than that used in GCUT and can output brain tissue with sharp boundary.
Although the proposed method has the above advantages, there are some disadvantages of the current implementation. First, the computational cost is higher than the other method. Second, two parameters are required to be manually set, if the parameters are set poorly, the extraction result might not be very desirable. Third, our proposed method was implemented as a 2D algorithm and could produce some rough artifacts on the brain surface. On the contrary, the deep learning technique-based brain extraction methods [21,22,24] usually run very fast and stably by using the GPU accelerating and learning techniques. However, the deep learning techniques also face a lot of challenges such as few-shot learning and capacity for transfer and must be supplemented by other techniques if we are to reach artificial general intelligence [33]. Furthermore, the U-Net based brain extraction methods [22,24] only used the original network, the network should be improved for the generalization capabilities across multiple datasets. Recently, Rundo et al. [34] incorporated Squeeze and Excitation blocks into U-Net (named USE-Net) for prostate zonal segmentation of multi-institutional MRI data sets. USE-Net provided excellent cross-dataset generalization when testing was performed on samples of the datasets used during training. As inspired by this, future work is to combine the improved U-Net technique with the ACN method, as the improved U-Net can provide more robust initial brain contour than BET across multiple data sets and predict more appropriate parameters for ACN.