The segmentation of the brain is a task commonly developed in neuroimaging laboratories. The difficulty and importance of the skull stripping problem has led to a wide range of proposals being developed to tackle it. Some techniques reported in the literature to solve this issue are, for example, surfacing models,1 deformable models,12.–3 watershed,4,5 morphology,6 atlas-based methods,7 hybrid techniques,8,9 fuzzy regions of interest,10 histogram analysis,11 active contours,12,13 multiresolution approach,14 multiatlas propagation and segmentation (MAPS),15,16 topological constraints,17 and others. Some revision papers concerning brain segmentation can be found in Refs. 1819.20.–21.
The problem that arises when having many viable techniques is to choose the ones that have the best performance for a particular visualization task. In Refs. 2223.–24, the authors selected the popular skull-stripping algorithms reported in the literature and carried out a comparison among them. These algorithms include Brain Extraction Tool (BET),3 3dIntracranial,25 Hybrid Watershed Algorithm,8 Brain Surface Extractor (BSE),26 and Statistical Parametric Mapping v.2 (SPM2).27 The two common and popular methods mentioned in Refs. 22–24 are BET and BSE. According to the results presented in Ref. 23, the authors found that BET and BSE produce similar brain extractions if adequate parameters are used in those algorithms. Interesting information about BET is reported in Ref. 28, where the authors found that the BET algorithm’s performance is improved after the removal of the neck slices. Due to the popularity of BET, we will compare our results to that algorithm. Two important characteristics of BET are: it is fast and it generates approximated segmentations.
In this paper, a morphological transformation that disconnects chained components is proposed and applied to segment brain from magnetic resonance imaging (MRI) T1. The operator is built as a composition between the viscous opening29 and the lower leveling30 and it is implemented in MATLAB R2010a on a 2.5 GHz Intel Core i5 processor with 2 GB RAM memory. To illustrate the performance of our proposals, two brain MRI datasets of 20 and 18 normal subjects,31 obtained from the Internet Brain Segmentation Repository (IBSR) and developed by the Centre for Morphometric Analysis (CMA) at the Massachusetts General Hospital ( http://www.nitrc.org), were processed.
In order to introduce our proposals, Sec. 2 provides a background on some morphological transformations such as opening and closing by reconstruction,32 viscous opening, and lower leveling. Other approaches on viscous transformations can be found in Refs. 3334.35.–36; however, these transformations work differently because they consider structuring elements that change dynamically, while our proposals work with a geodesic approach.29
In Sec. 3, a new transformation is built through the composition of the viscous opening and the leveling. Because the composed transformation uses several parameters, a simplification of it is introduced for facilitating its application. Such a reduction results in approximated segmentations and less time is utilized during its execution. It is noteworthy to mention that the two operators employ determinate size parameters deduced from a granulometric analysis.37
The experimental results are presented in Sec. 4. In Sec. 4.1, an explanation is given about the parameters involved in the proposed transformations and the performance of each one is illustrated with several pictures. In Sec. 4.2, the results obtained with the morphological transformations are compared using the mean values of the Jaccard38 and Dice39 indices with respect to those obtained from: (i) the BET algorithm and (ii) the results reported in Refs. 11, 12, and 40, which utilize the same databases. In Sec. 4.3, the advantages and disadvantages of our method for brain MR image extraction are presented. Section 5 contains our conclusions.
Background on Some Morphological Transformations
Opening and Closing by Reconstruction
In mathematical morphology (MM), the basic transformations are the erosion and dilation , where represents the three-dimensional (3-D) structuring element which has its origin in the center. Figure 1 illustrates the shape of the structuring element used in this paper. denotes the transposed set of with respect its origin, , is a size parameter, is the input image, and is a point on the definition domain.
The next equations represent the morphological erosion and dilation :41
In addition, the opening (closing) by reconstruction has the characteristic of modifying the regional maxima (minima) without affecting the remaining components to a large extent. These operators use the geodesic transformations.32,42 The geodesic dilation and erosion are expressed as with , and considering , respectively. When the function is equal to the morphological erosion or dilation, the opening or closing by reconstruction is obtained. Formally, the next expressions represent them
Viscous Opening and Viscous Difference
The viscous opening and closing defined in Ref. 29 allow one to deal with overlapped or chained components. These transformations are denoted by
Equation (2) uses three operators, the morphological erosion , the opening by reconstruction , and the morphological dilation . The morphological erosion allows one to discover and disconnect the -components (all components where the structuring element can go from one place to another by a continuous path made of squares and whose centers move along this path). Then, the opening by reconstruction removes all regions less than around the -components. Finally, because the viscous opening is defined on the lattice of dilatation, the must be obtained on . An example of Eq. (2) is given in Fig. 2. The original image is exhibited in Fig. 2(a). Figure 2(b) shows the morphological erosion . Notice that there are several components around the brain. The image in Fig. 2(c) corresponds to the transformation . In this image, several components have been eliminated by the process of opening by reconstruction. Figure 2(d) displays the result of the transformation .
Viscous openings permit the sieving of the image through the viscous difference.29 This is defined in
According to the explanation given above, the erosion discovers the -components and the difference with sieving the image, whereas is necessary to obtain the viscous component. The viscous difference gives the information of all discovered disconnected components of a certain size when is increased.
The lower leveling transformation is presented as follows:304) is iterated until stability is reached with the purpose of reconstructing the marker at the interior of the original mask , i.e.,
On the other hand, the selection of the marker is very important. In Ref. 30, for example, the following marker was used to segment the brain:
The parameter helps to control the reconstruction of the marker into . An example to illustrate the performance of Eq. (5) is given in the next section.
Segmentation of Brain in MRI T1
Equations (2), (5), and (6) were applied in Refs. 29 and 30 to separate the skull and the brain on slices of an MRI T1 for the two-dimensional (2-D) case. These transformations allow disconnecting overlapped components, because they can control the reconstruction process. In 2-D, the structuring element moves and uniquely touches one image. However, for the 3-D case, neighbors within the structuring element are taken from three adjacent brain slices. Due to this, several regions are connected through the shape of the 3-D structuring element among the different brain sections, resulting in, as a consequence, a major connectivity. The increase in connectivity originates that Eqs. (2) and (5) for the 3-D case do not show the same performance as in the 2-D case, and the component of the brain cannot be separated by applying such operators once.
This situation is illustrated in Fig. 3 where Eq. (5) has been applied considering the marker obtained by Eq. (6). The original volume appears in Fig. 3(a). Figure 3(b) displays a portion of the brain obtained from Eq. (6) with . The set of images in Figs. 3(c)–3(f) illustrates the control in the reconstruction process using different slopes . However, the extracted brain contains additional components since this condition is inadequate. The next section provides a solution to this problem.
Composition of Morphological Connected Transformations
As previously stated, viscous opening and lower leveling allow separating the chained components, and it is natural to think of combining both operators to get one transformation capable of having better control on the reconstruction process. Following this idea, an option is to use the viscous opening as a marker of the lower leveling to eliminate a great portion of the skull; posteriorly, the resulting image is again processed with a similar filter to eliminate the remaining regions around the brain. Such a procedure represents a sequential application of the combined transformations considering different parameters in order to have increasing control in the reconstruction process. The purpose is to eliminate the skull softly in two steps. The following equation permits the disconnection of the chained components and comes from the combination of Eqs. (2) and (5):
Nevertheless, Eq. (7) produces an unsatisfactory performance. Figure 4 shows an example of the transformation . Some slices of the original volume can be seen in Fig. 4(a) where the segmentation results in the creation of holes on the brain, as is illustrated in Fig. 4(b). The viscous opening causes this behavior, since all components not supporting the morphological erosion of size will merge with the background because the 3-D structuring element produces stronger changes as the size of the structuring element is increased. One way to get better segmentations consists of computing the operator in Eq. (8), where is the input image, represents the mean filter of size , and expresses a threshold between sections. The mean filter partially closes the holes, and permits the selection of certain regions of interest, and helps to obtain a portion of the original image
To apply Eq. (9), the reasoning below needs to be considered.
Parameters , , and
The quantification of the time when Eq. (5) is applied on a volume of 60 slides—using as a marker and the lower leveling with , 3, 6, 9, 12—is presented in Fig. 5. This figure displays several slices belonging to the output volumes obtained for different values.
ParameterFig. 6(a) corresponds to the application of Eq. (10) taking . This graph has three important intervals. The interval where shows the elimination of an important part of the skull. Hence, in order to detect the brain, will take values greater than 6, i.e., .
To apply Eq. (11), the next parameters will be considered: , , and (from the previous analysis for ), in order to detect the brain component.
Simplification of Eq. (9)
The fact of sequentially applying two transformations along with a mask to obtain Eq. (9) brings as a consequence the use of a large number of parameters. This problem is considered and a simplification is proposed as follows:
According to Eq. (12), the marker is obtained from the viscous opening, and it is propagated by the lower leveling transformation . Equation (12) presents the following benefits when compared with respect to Eq. (9): (1) the use of fewer parameters and (2) a reduction of the execution time. The performances of Eqs. (9) and (12) are illustrated in Sec. 4.
For the purpose of measuring the performance of our proposed method, the following MRI databases taken from the IBSR and developed by the CMA at the Massachusetts General Hospital ( http://www.nitrc.org)31 are utilized: (i) 20 simulated T1W MRI images (denoted as IBSR1) and (ii) 18 real T1W MRI images (denoted as IBSR2), with a slice thickness of 1.5 mm.
Tables 1Table 2Table 3–4 contain the parameters used in Eqs. (9) and (12) to segment the volumes belonging to the IBSR1 and IBSR2 datasets. In the volumes of IBSR1, the neck was cropped to obtain similar images to those of IBSR2. Differences among the volumes with respect to the intensity, size, and connectivity, originate the parameters’ variation. The guidelines for the parameter selection are given below:
Analysis for Eq. (9):
i. Parameters , , , and take their values into the interval [8, 14]. Furthermore, those volumes complying with and fulfill the order relation .
ii. The parameter takes its values within the interval [3, 12]. The mean filter will partially close several holes produced by the viscous opening as those presented in Fig. 4(b). Into MM, the closing transformation fills holes; however, from Fig. 4(b), the background and the holes represent the same region. This means that larger sizes of the structuring element will close the holes; nevertheless, this practice will increase the execution time of Eq. (9).
iii. The parameter represents a threshold. This varies in the interval [1, 120]. The application of a threshold obeys two things: (a) it eliminates some of the undesirable dark components such as dura matter, skin, and fat and (b) it obtains an appropriated marker.
iv. The lower leveling defined in Eq. (4) and utilized in Eq. (5) uses the parameters and . The size of the morphological dilation keeps its value during the processing with the purpose of detecting the different structures of the brain closer to the input image, whereas the slope varies in the interval [7, 120]. When parameter increases, a finer control is obtained, i.e., a smooth transition is generated in each iteration. Figure 7 shows an example of Eq. (9), considering the information of Table 1.
Analysis for Eq. (12):
v. The intervals defined previously are valid for Eq. (12). However, notice that the viscous opening and the lower leveling are applied once. For this situation, the viscous opening must get an appropriated marker containing the brain, and the lower leveling will reconstruct this marker inside the original volume. The transition (dura mater) between the skull and the brain avoids the lower leveling reconstructing the skull completely. Figure 8 displays an example of Eq. (12) by considering the information of Table 2. The input volume used to exemplify Eq. (12) was the same as that used in Fig. 7.
Parameters corresponding to Eq. (9). The processed dataset was IBSR1.
Parameters corresponding to Eq. (12). The processed dataset was IBSR1.
Parameters corresponding to Eq. (9). The processed dataset was IBSR2.
Parameters corresponding to Eq. (12). The processed dataset was IBSR2.
Figure 9 illustrates the resulting segmentation corresponding to the volume IBSR1_16 under the application of Eqs. (9), (12), and by BET (default parameters) implemented in the MRIcro software.43 Figure 9(a) shows the original slices. Figure 9(b) displays the respective manual segmentations. Figure 9(c) presents the application of Eq. (9) with the parameters given in Table 1. Figure 9(d) presents the application of Eq. (12) with the parameters given in Table 2. Figure 9(e) shows the brain extraction using the BET algorithm. The parameters selected for the set of 20 brains are intensity and . In order to compare the segmentations, the Jaccard and the Dice coefficients are computed. Table 5 contains the indices corresponding to BET and Eqs. (12) and (9) for the IBSR1.
Jaccard and Dice indexes corresponding to IBSR1 dataset and segmented with BET, Eqs. (9) and (12) considering the parameters presented in Tables 1 and 2.
|Volume||BET||Eq. (12)||Eq. (9)|
A similar procedure is applied to the IBSR2 database. Figure 10 presents the segmentations corresponding to the IBSR2_04 volume. Figure 10(a) shows the original slices. Figure 10(b) displays the respective manual segmentations. Figure 10(c) presents the application of Eq. (9) with the parameters given in Table 3. Figure 9(d) presents the application of Eq. (12) with the parameters given in Table 4. Figure 9(e) shows the brain extraction using the BET algorithm with the default parameters. Table 6 contains the Jaccard and Dice indices corresponding to BET and Eqs. (12) and (9) for the IBSR2.
Jaccard and Dice indexes corresponding to IBSR2 dataset and segmented with BET, Eqs. (9) and (12) considering the parameters presented in Tables 3 and 4.
|Volume||BET||Eq. (12)||Eq. (9)|
Jaccard and Dice mean indexes corresponding to IBSR1 and IBSR2 datasets and some reported in the literature.
|Method||Jaccard mean||Dice mean|
|SMHASS (Ref. 11) IBSR1||0.904||0.950|
|SMHASS (Ref. 11) IBSR2||0.905||0.950|
|ACNM One (Ref. 12) IBSR1||0.890||0.940|
|ACNM One (Ref. 12) IBSR2||0.900||0.950|
|Equation (9) IBSR1||0.935||0.966|
|Equation (9) IBSR2||0.924||0.963|
|Equation (12) IBSR1||0.866||0.938|
|Equation (12) IBSR2||0.919||0.957|
|Method in Ref. 40 using IBSR1||0.923||0.960|
Some commentaries on the segmented volumes are presented as follows:
i. The time performance of our operators is slow compared to the BET algorithm. The table in Fig. 11(a) presents several times measured during the execution of Eq. (9) and considering an increasing numbers of slices. Its corresponding graph is shown in Fig. 9(b). A similar behavior is presented using Eq. (12), but within half the time.
In Ref. 24, the measured time to separate brain components varies between 40 s (BET algorithm) to 35 min (SPM2) when considering a complete volume. In our case, Eq. (9) takes 16 min and Eq. (12) took 8 min.
iii. Although the time spent to segment a volume is high compared to the BET, the segmentations obtained with our proposal [Eq. (9)] are better according to the Jaccard and Dice mean indices.
Two morphological transformations were proposed to extract brain in MRIs T1. The first operator [Eq. (9)] presents a better performance than the second one [Eq. (12)] according to the computed Jaccard and Dice mean indices.
The idea to segment the brain consisted of smoothly propagating a marker given by the viscous opening into the original volume. For this, adequate parameters must be obtained from a granulometric analysis. The sequential application of such transformations results in a new morphological operator [Eq. (9)] capable of better controlling the reconstruction process. Nevertheless, the new transformation employs several parameters; due to this, a second morphological transformation was obtained from a simplification of the first one [Eq. (12)].
Our proposals were tested using the two brain databases obtained from the IBSR home page. In total, 38 volumes of MR images of the brain were processed. The segmentations were compared through two popular indices with respect to the segmentations obtained from the BET algorithm, with manual segmentations obtained from the IBSR website and with respect to the values of the indices reported in the current literature. When the mean values of the Jaccard and Dice indices are compared, our proposal outperforms the other methodologies. This means that our segmentations are closer to the manual segmentations obtained from the IBSR website. However, the time spent to segment a volume with 160 slices, along with the number of parameters utilized in Eq. (9), is higher compared to the time and the parameters utilized by the BET algorithm.
The main problem of the BET algorithm is that several regions are not detected; the beginning of Fig. 10(e) clearly illustrates this situation. Due to this, the Jaccard and Dice indices fall considerably.
Finally, as future work, the proposal presented in this paper will be improved by the implementation of fast algorithms and/or parallel implementation on graphics processing units using compute unified device architectures technology, so the performance can improve for real-time applications.
The author Iván R. Terol-Villalobos would like to thank Diego Rodrigo and Darío T.G. for their great encouragement. This work was funded by the government agency CONACyT México under the Grant 133697.
Jorge Domingo Mendiola-Santibañez received his PhD degree from the Universidad Autónoma de Querétaro (UAQ), México. Currently, he is a professor/researcher at the Universidad Autónoma de Querétaro. His research interests include morphological image processing and computer vision.
Martín Gallegos-Duarte is an MD and a PhD student at the Universidad Autonoma de Querétaro. He is head of the Strabismus Service at the Institute for the Attention of Congenital Diseases and Ophthalmology-Pediatric Service in the Mexican Institute of Ophthalmology in the state of Queretaro, Mexico.
Miguel Octavio Arias-Estrada is a researcher in computer science at National Institute of Astrophysics, Optics and Electronics, Puebla, Mexico, with a PhD degree in electrical engineering (computer vision) from Laval University (Canada) and BEng and MEng degrees in electronic engineering from University of Guanajuato (Mexico). Currently, he is a researcher at INAOE (Puebla, México). His interests are computer vision, FPGA and GPU algorithm acceleration for three-dimensional machine vision.
Israel Marcos Santillán-Méndez received the BS degree in engineering from the Instituto Tecnológico de Estudios Superiores de Monterrey and his MS degree and PhD degree in engineering from Facultad de Ingeniería de la Universidad Autónoma de Querétaro (México). His research interests include models of biological sensory and perceptual systems and mathematical morphology.
Juvenal Rodríguez-Reséndiz received his MS degree in automation control from University of Querétaro and PhD degree at the same institution. Since 2004, he has been part of the Mechatronics Department at the UAQ. He is the head of the Automation Department. His research interest includes signal processing and motion control. He serves as vice president of IEEE in Queretaro State.
Iván Ramón Terol-Villalobos received his BSc degree from Instituto Politécnico Nacional (I.P.N. México), his MSc degree from Centro de Investigación y Estudios Avanzados del I.P.N. (México). He received his PhD degree from the Centre de Morphologie Mathématique, Ecole des Mines de Paris (France). Currently, he is a researcher at CIDETEQ (Querétaro, México). His main current research interests include morphological image processing, morphological probabilistic models, and computer vision.