# Unsupervised Neural Techniques Applied to MR Brain Image Segmentation

Authors: A. Ortiz; J. M. Gorriz; J. Ramirez; D. Salas-Gonzalez
Journal: Advances in Artificial Neural Systems
 (2012)

---

## Abstract

The primary goal of brain image segmentation is to partition a given brain image into different regions representing anatomical structures. Magnetic resonance image (MRI) segmentation is especially interesting, since accurate segmentation in white matter, grey matter and cerebrospinal fluid provides a way to identify many brain disorders such as dementia, schizophrenia or Alzheimer’s disease (AD). Then, image segmentation results in a very interesting tool for neuroanatomical analyses. In this paper we show three alternatives to MR brain image segmentation algorithms, with the Self-Organizing Map (SOM) as the core of the algorithms. The procedures devised do not use any a priori knowledge about voxel class assignment, and results in fully-unsupervised methods for MRI segmentation, making it possible to automatically discover different tissue classes. Our algorithm has been tested using the images from the Internet Brain Image Repository (IBSR) outperforming existing methods, providing values for the average overlap metric of 0.7 for the white and grey matter and 0.45 for the cerebrospinal fluid. Furthermore, it also provides good results for high-resolution MR images provided by the Nuclear Medicine Service of the “Virgen de las Nieves” Hospital (Granada, Spain).


---

## Body

## 1. Introduction

Nowadays, magnetic resonance imaging (MRI) systems  provide an excellent spatial resolution as well as a high tissue contrast. Nevertheless, since actual MRI systems can obtain 16-bit depth images corresponding to 65535 gray levels, the human eye is not able to distinguish more than several tens of gray levels. On the other hand, MRI systems provide images as slices which compose the 3D volume. Thus, computer-aided tools are necessary to exploit all the information contained in an MRI. These are becoming a very valuable tool for diagnosing some brain disorders such as Alzheimer’s disease [1–5]. Moreover, modern computers, which contain a large amount of memory and several processing cores, have enough process capabilities for analyzing the MRI in reasonable time.Image segmentation consists in partitioning an image into different regions. In MRI, segmentation consists of partitioning the image into different neuroanatomical structures which corresponds to different tissues. Hence, analyzing the neuroanatomical structures and the distribution of the tissues on the image, brain disorders or anomalies can be figured out. Hence, the importance of having effective tools for grouping and recognizing different anatomical tissues, structures and fluids is growing with the improvement of the medical imaging systems. These tools are usually trained to recognize the three basic tissue classes found on a healthy brain MR image: white matter (WM), gray matter (GM), and cerebrospinal fluid (CSF). All of the nonrecognized tissues or fluids are classified as suspect, to be pathological.The segmentation process can be performed in two ways. The first consists of manual delimitation of the structures present within an image by an expert. The second consists of using an automatic segmentation technique. As commented before, computer image processing techniques allow exploiting all the information contained in an MRI.There are several automatic segmentation techniques. Some of them use the information contained in the image histogram [6–11]. This way, since different contrast areas should correspond with different tissues, the image histogram can be used for partitioning the image. Nevertheless, variations on the contrast of the same tissue are found in an image due to RF noise or shading effects due to magnetic field variations, resulting in tissue misclassification. Other methods use statistical classifiers based on the expectation-maximization (EM) algorithms [12–14], maximum likelihood (ML) estimation [15], or Markov random fields [16, 17]. Other segmentation techniques are based on artificial neural network classifiers [8, 18–21] such as self-organizing maps (SOMs) [18, 19, 21–23].In this paper we present three segmentation alternatives based on SOMs, which provide good results over the internet brain image repository (IBSR) [16] images.


## 2. SOM Algorithm

SOM is an unsupervised classifier proposed by Kohonen and it has been used for a large number of applications regarding classification or modelling [24]. The self-organizing process is based on the distance (usually the Euclidean distance) computation among each training sample and all the units on the map as a part of a competitive learning process. On the other hand, several issues such as topological map, number of units on the map, initialization of weights, and the training process on the map are decisive for the classification quality. Regarding the topology, a 2D hexagonal grid was selected since it fitted better in the feature space as shown in the experiments.The SOM algorithm can be summarized as follows. LetX⊂ℝd be the data manifold. In each iteration, the winning unit is computed according to(1)Uω(t)=argmin⁡i⁡{‖x(t)-ωi(t)‖},
where x(t), x∈X, is the input vector at time t and ωi(t) is the prototype vector associated with the unit i. The unit closer to the input vector Uω(t) is referred to as winning unitand the associated prototype is updated. To complete the adaptive learning process on the SOM, the prototypes of the units in the neighborhood of the winning unit are also updated according to:(2)ωi(t+1)=ωi(t)+α(t)hUi(t)(x(t)-ωi(t)),
where α(t) is the exponential decay learning factor and hUi(t) is the neighborhood function associated with the unit i. Both, the learning factor and the neighborhood function decay with time; thus the prototypes adaptation becomes slower as the neighborhood of the unit i contains less number of units:(3)hUi(t)=e(-‖rU-ri‖2/2σ(t)2).Equation (3) shows the neighbourhood function, where ri  represents the position on the output space and ∥rU-ri∥ is the distance between the winning unit and the unit i on the output space. The neighbourhood is defined by a Gaussian function which shrinks in each iteration as shown in (4). In this competitive process, the winning unit is named the best matching unit (BMU). On the other hand, σ(t) controls the reduction of the Gaussian neighborhood in each iteration. τ1 is a time constant which depends on the number of iterations and the map radius and computed as τ1 = number_of_iterations/map_radius:(4)σ(t)=σ0e(-t/τ1).The quality of the trained map can be computed by the means of two measures. These two measures are the quantization error(te), which determines the average distance between each data vector and its best matching unit (BMU) and the topological error (qe), which measures the proportion of all data vectors for which first and second BMUs are not adjacent units. Both, the quantization error and the topological error are defined by the following:(5)te=1N∑i=1Nu(x⃗i),(6)qe=∑i=1N‖x⃗i-b⃗x⃗i‖.In (5), N is the total number of data vectors, and u(x⃗i) is 1 if the first and the second BMU for x⃗i are nonadjacent and 0 otherwise. In (6) the quantization error is defined where x⃗i is the ith data vector on the input space and bω⃗i is the weight (prototype) associated with the best matching unit for the data vector x⃗i. Therefore, lower values of te and qe imply a better topology preservation, which is equivalent to a better clustering result. That is to say, the lower the values on the quantization error (qe) and the topological error (te), the better the goodness of the SOM [25, 26]. In this paper, SOM toolbox [27] has been used to implement SOM.


## 3. MR Image Segmentation with SOM

In this section we present two image segmentation algorithms based on unsupervised SOM. The first uses the histogram to segment the whole volume (i.e., classify all the voxels on the volumetric image). The second extracts a set of features from each image slice and uses an SOM to classify the feature vectors into clusters using the devised entropy gradient clustering method. Thus, Figure1 shows the block diagram of the presented segmentation algorithms.Figure 1
Block diagram of the segmentation methods.


### 3.1. Image Preprocessing

Once the MR image has been acquired, a preprocessing is performed in order to remove noise and to homogenize the image background. The brain extraction for undesired structures removal (i.e., skull and scalp) can be done at this stage. There are several algorithms for this purpose such as brain surface extractor (BSE), brain extraction tool (BET) [8], Minneapolis consensus strip (McStrip), or hybrid watershed algorithm (HWA) [2]. Since IBSR 1.0 images have these undesired structures already removed, brain extraction is not required. Nevertheless, images provided by IBSR 2.0 are distributed without the scalp/skull already removed. In these images, the brain has been extracted in the preprocessing stage using BET.


### 3.2. Segmentation Using the Volume Image Histogram (HFS-SOM)

The first step after preprocessing the image consists in computing the volume image histogram which describes the probability of occurrence of voxel intensities in the volume image and provides information regarding different tissues. A common approach to avoid processing the large number of voxels present on MR images consists in modelling the intensity values as a finite number of prototypes, which deals to improve the computational effectiveness. After computing the histogram, the bin 0 is removed since it contains all the background voxels. Thus, only information corresponding to the brain is stored.Figure2 shows the rendered brain surface from the IBSR volume 12, and its histogram.Rendered brain surface extracted from IBSR_12 volume (a) and computed histogram (b).
(a)(b)Histogram data including the intensity occurrence probabilities(pi) and the relative position (bin number), bi,are used to compose the feature vectors F⃗=(pi,bi),pi∈ℝ,bi∈ℤ, to be classified by the SOM.On a trained SOM, the output layer is composed by a reduced number of prototypes (the number of units on the output layer) modelling the input data manifold. In addition, the most similar prototypes are closely located in the output map at the time the most dissimilar are located apart. Nevertheless, since all the units have an associated prototype, it is necessary to cluster the SOMs in order to define the borders between clusters. In other words, each prototype is grouped so that it belongs to a cluster. Thus, thek-means algorithm is used to cluster the SOMs, grouping the prototypes into a number of different classes, and the DBI [28], which gives lower values for better clustering results, is computed for different k values to provide a measurement of the clustering validity.The clusters on the SOM group the units so that they belong to a specific class. As each of these units will be the BMU of a specific set of voxels, the clusters define different voxel classes. This way, each voxel is labeled as belonging to a class (i.e., segment).


### 3.3. MR Image Segmentation with SOM and the Entropy-Gradient Algorithm (EGS-SOM)

The method described in this section is also based on SOM for voxel classification, but histogram information from the image volume is replaced by computing a set of features, selecting the most discriminant ones. After that, SOM clustering is performed by the EGS-SOM method described here in after, which allows us to obtain higher-resolution images providing good segmentation results as shown in the experiments.


#### 3.3.1. Feature Extraction and Selection

In this stage some significant features from the MR image are extracted to be subjected to classification. As commented before, we perform the image processing slice by slice on each plane. Thus, the feature extraction is carried out by using an overlapping and sliding window of7×7 pixels on each slice of a specific plane.In the feature extraction process, window size plays an important role since smaller windows are not able to capture the second-order features, that is, texture information. The use of higher window sizes results in loosing resolution. Therefore the7×7 size gives a good trade-off between complexity and performance.In this paper we use first- and second-order statistical features [29]. The first-order features we extract from the image are intensity, mean, and variance. The intensity is referred to the gray level of the center pixel on the window. The mean and variance are calculated taking into account the gray level present on the window. On the other hand, we additionally use second-order features such as textural features. Haralick et al. [30] proposed the use of 14 features for image classification, computed using the gray level coocurrence matrix (GLCM) method. The set of second-order features we have used are energy, entropy, contrast, angular second moment (ASM), sum average, autocorrelation, correlation, inverse difference moment, maximum probability, cluster prominence, cluster shade, dissimilarity, and second-order variance as well as moment invariants [31].In order to select the most discriminant features, a genetic algorithm is used to minimize the topological and the quantization error on the SOM through the fitness function shown in (7)(7)FQT=(0.5qe+0.5te).The feature selection process is summarized in Figure3.Figure 3
Feature selection process with genetic algorithm.The stop criterion is reached when the performance of the proposed solutions does not improve the performance significantly (1%) or the maximum number of generations is reached (500).Once the dimension of the feature space has been reduced, we use the vectors of this space for training a SOM. The topology of the map and the number of units on the map are decisive for the SOM quality. In that sense, we use a hexagonal grid since it allows better fitting the prototypes to the feature space vectors. Each BMU on the SOM has an associated pixel on the image. This association is made through a matrix computed during the feature extraction phase which stores the coordinates of the central pixel on each window. This allows associating a feature vector to an image pixel.Nevertheless, these clusters roughly define the different areas (segments) on the image, and a further fine-tuning phase is required. This fine-tuning phase is accomplished by the entropy-gradient method.Entropy-Gradient Method
The procedure devised consists of using the feature vectors associated with each BMU to compute a similarity measurement among the vectors belonging with each BMU and the vectors associated to each other BMU. Next, the BMUs are sorted in ascending order of the contrast. Finally, the feature vectors of each BMU are included on a cluster. For each map unit, we compute the accumulated entropy:(8)Hmi=∑n=1NpHn,
where i is the map unit index and  Np the number of pixels belonging to the map unit in the classification process. This means that the unit i has a number of Np-associated pixels. Since the output layer on the SOM is a two-dimensional space, we calculate the entropy-gradient vector from each map unit (8) and move to the opposite direction for clustering.


## 3.1. Image Preprocessing

Once the MR image has been acquired, a preprocessing is performed in order to remove noise and to homogenize the image background. The brain extraction for undesired structures removal (i.e., skull and scalp) can be done at this stage. There are several algorithms for this purpose such as brain surface extractor (BSE), brain extraction tool (BET) [8], Minneapolis consensus strip (McStrip), or hybrid watershed algorithm (HWA) [2]. Since IBSR 1.0 images have these undesired structures already removed, brain extraction is not required. Nevertheless, images provided by IBSR 2.0 are distributed without the scalp/skull already removed. In these images, the brain has been extracted in the preprocessing stage using BET.


## 3.2. Segmentation Using the Volume Image Histogram (HFS-SOM)

The first step after preprocessing the image consists in computing the volume image histogram which describes the probability of occurrence of voxel intensities in the volume image and provides information regarding different tissues. A common approach to avoid processing the large number of voxels present on MR images consists in modelling the intensity values as a finite number of prototypes, which deals to improve the computational effectiveness. After computing the histogram, the bin 0 is removed since it contains all the background voxels. Thus, only information corresponding to the brain is stored.Figure2 shows the rendered brain surface from the IBSR volume 12, and its histogram.Rendered brain surface extracted from IBSR_12 volume (a) and computed histogram (b).
(a)(b)Histogram data including the intensity occurrence probabilities(pi) and the relative position (bin number), bi,are used to compose the feature vectors F⃗=(pi,bi),pi∈ℝ,bi∈ℤ, to be classified by the SOM.On a trained SOM, the output layer is composed by a reduced number of prototypes (the number of units on the output layer) modelling the input data manifold. In addition, the most similar prototypes are closely located in the output map at the time the most dissimilar are located apart. Nevertheless, since all the units have an associated prototype, it is necessary to cluster the SOMs in order to define the borders between clusters. In other words, each prototype is grouped so that it belongs to a cluster. Thus, thek-means algorithm is used to cluster the SOMs, grouping the prototypes into a number of different classes, and the DBI [28], which gives lower values for better clustering results, is computed for different k values to provide a measurement of the clustering validity.The clusters on the SOM group the units so that they belong to a specific class. As each of these units will be the BMU of a specific set of voxels, the clusters define different voxel classes. This way, each voxel is labeled as belonging to a class (i.e., segment).


## 3.3. MR Image Segmentation with SOM and the Entropy-Gradient Algorithm (EGS-SOM)

The method described in this section is also based on SOM for voxel classification, but histogram information from the image volume is replaced by computing a set of features, selecting the most discriminant ones. After that, SOM clustering is performed by the EGS-SOM method described here in after, which allows us to obtain higher-resolution images providing good segmentation results as shown in the experiments.


### 3.3.1. Feature Extraction and Selection

In this stage some significant features from the MR image are extracted to be subjected to classification. As commented before, we perform the image processing slice by slice on each plane. Thus, the feature extraction is carried out by using an overlapping and sliding window of7×7 pixels on each slice of a specific plane.In the feature extraction process, window size plays an important role since smaller windows are not able to capture the second-order features, that is, texture information. The use of higher window sizes results in loosing resolution. Therefore the7×7 size gives a good trade-off between complexity and performance.In this paper we use first- and second-order statistical features [29]. The first-order features we extract from the image are intensity, mean, and variance. The intensity is referred to the gray level of the center pixel on the window. The mean and variance are calculated taking into account the gray level present on the window. On the other hand, we additionally use second-order features such as textural features. Haralick et al. [30] proposed the use of 14 features for image classification, computed using the gray level coocurrence matrix (GLCM) method. The set of second-order features we have used are energy, entropy, contrast, angular second moment (ASM), sum average, autocorrelation, correlation, inverse difference moment, maximum probability, cluster prominence, cluster shade, dissimilarity, and second-order variance as well as moment invariants [31].In order to select the most discriminant features, a genetic algorithm is used to minimize the topological and the quantization error on the SOM through the fitness function shown in (7)(7)FQT=(0.5qe+0.5te).The feature selection process is summarized in Figure3.Figure 3
Feature selection process with genetic algorithm.The stop criterion is reached when the performance of the proposed solutions does not improve the performance significantly (1%) or the maximum number of generations is reached (500).Once the dimension of the feature space has been reduced, we use the vectors of this space for training a SOM. The topology of the map and the number of units on the map are decisive for the SOM quality. In that sense, we use a hexagonal grid since it allows better fitting the prototypes to the feature space vectors. Each BMU on the SOM has an associated pixel on the image. This association is made through a matrix computed during the feature extraction phase which stores the coordinates of the central pixel on each window. This allows associating a feature vector to an image pixel.Nevertheless, these clusters roughly define the different areas (segments) on the image, and a further fine-tuning phase is required. This fine-tuning phase is accomplished by the entropy-gradient method.Entropy-Gradient Method
The procedure devised consists of using the feature vectors associated with each BMU to compute a similarity measurement among the vectors belonging with each BMU and the vectors associated to each other BMU. Next, the BMUs are sorted in ascending order of the contrast. Finally, the feature vectors of each BMU are included on a cluster. For each map unit, we compute the accumulated entropy:(8)Hmi=∑n=1NpHn,
where i is the map unit index and  Np the number of pixels belonging to the map unit in the classification process. This means that the unit i has a number of Np-associated pixels. Since the output layer on the SOM is a two-dimensional space, we calculate the entropy-gradient vector from each map unit (8) and move to the opposite direction for clustering.


## 3.3.1. Feature Extraction and Selection

In this stage some significant features from the MR image are extracted to be subjected to classification. As commented before, we perform the image processing slice by slice on each plane. Thus, the feature extraction is carried out by using an overlapping and sliding window of7×7 pixels on each slice of a specific plane.In the feature extraction process, window size plays an important role since smaller windows are not able to capture the second-order features, that is, texture information. The use of higher window sizes results in loosing resolution. Therefore the7×7 size gives a good trade-off between complexity and performance.In this paper we use first- and second-order statistical features [29]. The first-order features we extract from the image are intensity, mean, and variance. The intensity is referred to the gray level of the center pixel on the window. The mean and variance are calculated taking into account the gray level present on the window. On the other hand, we additionally use second-order features such as textural features. Haralick et al. [30] proposed the use of 14 features for image classification, computed using the gray level coocurrence matrix (GLCM) method. The set of second-order features we have used are energy, entropy, contrast, angular second moment (ASM), sum average, autocorrelation, correlation, inverse difference moment, maximum probability, cluster prominence, cluster shade, dissimilarity, and second-order variance as well as moment invariants [31].In order to select the most discriminant features, a genetic algorithm is used to minimize the topological and the quantization error on the SOM through the fitness function shown in (7)(7)FQT=(0.5qe+0.5te).The feature selection process is summarized in Figure3.Figure 3
Feature selection process with genetic algorithm.The stop criterion is reached when the performance of the proposed solutions does not improve the performance significantly (1%) or the maximum number of generations is reached (500).Once the dimension of the feature space has been reduced, we use the vectors of this space for training a SOM. The topology of the map and the number of units on the map are decisive for the SOM quality. In that sense, we use a hexagonal grid since it allows better fitting the prototypes to the feature space vectors. Each BMU on the SOM has an associated pixel on the image. This association is made through a matrix computed during the feature extraction phase which stores the coordinates of the central pixel on each window. This allows associating a feature vector to an image pixel.Nevertheless, these clusters roughly define the different areas (segments) on the image, and a further fine-tuning phase is required. This fine-tuning phase is accomplished by the entropy-gradient method.Entropy-Gradient Method
The procedure devised consists of using the feature vectors associated with each BMU to compute a similarity measurement among the vectors belonging with each BMU and the vectors associated to each other BMU. Next, the BMUs are sorted in ascending order of the contrast. Finally, the feature vectors of each BMU are included on a cluster. For each map unit, we compute the accumulated entropy:(8)Hmi=∑n=1NpHn,
where i is the map unit index and  Np the number of pixels belonging to the map unit in the classification process. This means that the unit i has a number of Np-associated pixels. Since the output layer on the SOM is a two-dimensional space, we calculate the entropy-gradient vector from each map unit (8) and move to the opposite direction for clustering.


## 4. Results and Discussion

In this section we show the segmentation results obtained using real MR brain images from two different sources. One of these sources is the IBSR database [32] in two versions, IBSR and IBSR 2.0.Figures4(a) and 4(b) show the segmentation results for the IBSR volume 100_23 using the HFS-SOM algorithm and the EGS-SOM algorithm, respectively. In these images, WM, GM, and CSF are shown for slices 120, 130, 140, 150, 160, and 170 on the axial plane. Expert segmentation from IBSR database is shown in Figure 4(c).Segmentation of the IBSR volume 100_23 using the HFS-SOM algorithm (a) and the EGS-SOM algorithm (b). Ground Truth is shown in (c). Slices 120, 130, 140, 150, 160, and 170 on the axial plane are shown on each column. First column corresponds to WM, second column to GM, and third column to CSF.
(a)(b)(c)Figure5(a) shows the segmentation results for the IBSR 2.0 volume 12 using the fast volume segmentation algorithm. In this figure, each row corresponds to a tissue and each image column corresponds to a different slice. In the same way, Figure 5(b) shows the same slices of Figure 5(b) but the segmentation is performed using the EGS-SOM algorithm. Figure 5(c) shows the segmentation performed by expert radiologists provided by the IBSR database (ground truth).Segmentation of the IBSR 2.0 volume 12 using the HFS-SOM algorithm (a) and the EGS-SOM algorithm (b). Ground Truth is shown in (c). Slices 110, 120, 130, 140, 150, and 160 on the axial plane are shown on each column. First column corresponds to WM, second column to GM, and third column to CSF.
(a)(b)(c)Visual comparison between automatic segmentation and the ground truth points up that the EGS-SOM method outperforms the fast volume segmentation method.This fact is also stated in Figure6 where Tanimoto’s index is shown for different segmentation algorithms, where SSOM corresponds to our entropy-gradient algorithm, BMAP is biased map [33], AMAP is adaptative map [33], MAP is maximum a posteriori probability [34], MLC is maximum likelihood [35], FUZZY is fuzzy k-means [36] and TSKMEANS is tree-structured k-means [36]. The performance of the presented segmentation techniques has been evaluated by computing the average overlap rate through Tanimoto’s index, as it has been widely used by other authors to compare the segmentation performance of their proposals [13, 16, 17, 21, 26, 37–41]. Tanimoto’s index can be defined as(9)T(S1,S2)=|S1∩S2||S1∪S2|,
where S1 is the segmentation set and  S2 is the ground truth.Figure 6
Average overlap metric comparison for different segmentation methods.


## 5. Conclusions

In this paper we presented fully unsupervised segmentation methods for MR images based on hybrid artificial intelligence techniques for improving the feature extraction process and self-organizing maps for pixel classification. The use of a genetic algorithm provides a way for training the Self-Organizing map used as a classifier in the most efficient way. This is because the dimension of the training samples (feature vectors) has been reduced in order to be enough discriminant but not redundant. As a result, the number of units (neurons) on the map is also optimized as well as the classification process. Thus, we take advantage of the competitive learning model of the SOM which groups the pixels into clusters. This competitive process discovers similarities among the pixels, resulting in an unsupervised way to segment the image. Moreover, the clusters’ borders are redefined by using the entropy-gradient method presented on this paper. The whole process allows figuring out the segments present on the image without using any a priori information.The results shown in Section4 have been compared with the segmentations provided by the IBSR database that outperform the results obtained by other algorithms such as k-means or fuzzy k-means. The number of segments or different tissues found in an MR image is figured out automatically making possible to find out tissues which could be identified with a pathology.



---

*Source: 457590-2012-06-07.xml*