Abstract
We propose a novel segmentation method based on regional and nonlocal information to overcome the impact of image intensity inhomogeneities and noise in human brain magnetic resonance images. With the consideration of the spatial distribution of different tissues in brain images, our method does not need preestimation or precorrection procedures for intensity inhomogeneities and noise. A nonlocal information based Gaussian mixture model (NGMM) is proposed to reduce the effect of noise. To reduce the effect of intensity inhomogeneity, the multigrid nonlocal Gaussian mixture model (MNGMM) is proposed to segment brain MR images in each nonoverlapping multigrid generated by using a new multigrid generation method. Therefore the proposed model can simultaneously overcome the impact of noise and intensity inhomogeneity and automatically classify 2D and 3D MR data into tissues of white matter, gray matter, and cerebral spinal fluid. To maintain the statistical reliability and spatial continuity of the segmentation, a fusion strategy is adopted to integrate the clustering results from different grid. The experiments on synthetic and clinical brain MR images demonstrate the superior performance of the proposed model comparing with several state-of-the-art algorithms.
1. Introduction
Magnetic resonance imaging (MRI) is a helpful method for diagnosis of brain diseases, such as Alzheimer’s disease and schizophrenia. Accurate tissues segmentation, including gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF), plays an important role in clinical practice and hence has attracted extensive research attention.
Many methods have been proposed for MR image segmentation in the past several decades. These approaches can be classified in terms of different criteria. For example, edge based methods [1, 2], region based methods [3, 4], and clustering based methods [5–7]. Unfortunately, most segmentation methods are hindered by various imaging artifacts such as noise and intensity inhomogeneities.
Intensity inhomogeneity, also known as bias field, arises from the imperfections of the image acquisition process and changes the absolute intensity for a given tissue class in different locations, which usually makes the intensity distribution within a particular tissue class flatter. Most traditional intensity based methods cannot obtain satisfactory results due to the impact of intensity inhomogeneity.
The observed MRI signal is the product of the true signal generated by the underlying anatomy and a spatially varying field factor and an additive noise :where , , , and are the observed intensity, true intensity, bias field, and noise at the th voxel, respectively. is the total number of pixels in the MR image. Many techniques [6, 8, 9] often ignore the noise and take the logarithmic transform on both sides of (1):
Many methods have been proposed to correct or estimate intensity inhomogeneities. Collewet et al. proposed a method based on measuring the coil sensitivity functions [8]. Based on the observation that the bias field is smooth, another group of methods overcome the impact of intensity inhomogeneity without estimating the bias field [9, 10]. However, most of them may lose edge information [11].
In this paper, we first propose an improved nonlocal Gaussian mixture model by introducing the nonlocal information into GMM model to reduce the effect of the noise. Then, a new multigrid generation method is presented, and we simplify the NGMM into a local version to eliminate the effect of the bias field and the intensity variation of intratissues. Finally, we propose a fusion strategy to integrate the results from different local regions. The experiments on both synthetic and clinical brain MR images show that our method can obtain more accurate results.
2. Nonlocal Gaussian Mixture Model
Gaussian mixture model (GMM) has been widely used in many applications due to its excellent approximation properties. Suppose a MR image has a mixture of components, the mixture density function of pixel can be written as where is the mixture weights and is the parameter vector. is a standard Gaussian distribution of the th component and contains the parameters of the Gaussian distribution. and are the mean and variance, respectively. Then the entire distribution can be written as
The problem is how to find the best parameters :
Equation (5) can be calculated by using the expectation-maximization (EM) method [7]. In the E-Step, the algorithm calculates the expected value of the th weight:
In the M-Step, the parameters of the th Gaussian distribution can be calculated:
Based on the initialization, and are calculated iteratively until the stop criteria are reached. Finally, the pixel can be classified into the th class when . From (3), we can find that the GMM only considers the intensity distribution, which makes the method sensitive to the intensity inhomogeneity and noise.
In order to reduce the effect of the intensity inhomogeneity, Wells et al. [6] proposed a method to simultaneously estimate the bias field and segment the image into different classes. However, the method only addressed the bias field without analyzing the inhomogeneities in inner tissues. In order to ensure the smoothness of the bias field, a low-pass filter is utilized to convolve the bias field, which makes the estimated bias field inaccurate. Furthermore, the method is sensitive to the noise.
In order to reduce the effect of the noise and preserve more detailed information, we improve the Gaussian mixture model by using nonlocal information. The nonlocal information has been widely used for denoising purposes [12, 13]. Following the idea of nonlocal means method [12], we use the nonlocal information to adapt which can be defined aswhere is the weight function based on the similarity between the neighbor patch of each neighbor pixel to that of center pixel and satisfies the conditions and . For pixel and its neighbor pixel , the weight function is defined aswhere and are the neighbor patches around pixels and with width . Broadly speaking, if the neighbor patches of two pixels and are similar, it is more probable that these pixels belong to the same tissue, and the corresponding values of weight function would be higher. Conversely, if these two pixels are quite different, the value of the weight function should be small. acts as a filtering parameter to control the decay of the exponential function. is the Euclidean distance and is the standard deviation of the Gaussian kernel. Due to the fast decay of the exponential kernel, large distances between estimated patches lead to nearly zero weights. Essentially, the weight function aims to take advantage of the redundancy present in natural structures. Therefore, by using nonlocal information, the nonlocal Gaussian mixture model can reduce the effect of noise and preserve details of the edges.
3. Multigrid Nonlocal Gaussian Mixture Method (MGMM)
Without considering the intensity inhomogeneity, the nonlocal Gaussian mixture model can only reduce the effect of noise but cannot obtain satisfactory results for the image containing severe intensity inhomogeneities. Zhu and Jiang [14] proposed a multicontext fuzzy clustering method (MCFC) to reduce the effect of intensity inhomogeneity by using fuzzy clustering method on each nonoverlapping regions and a fusion strategy to integrate the clustering outcomes form different regions. However, this method is sensitive to noise and it only uses traditional multigrid generation method, which makes the method inaccurate. Following the idea of MCFC, we improve the GMM by using the nonoverlapping multigrid. This idea is based on the following four assumptions:(1)The brain image has been skull-stripped. In this paper, we use the cut based method [15].(2)Bias field is smooth and slowly varying.(3)Within each grid, the number of clusters must equal and there are considerable numbers of pixels in each tissue class.(4)Within a grid, all pixels of the same tissue have similar true intensities.
The brain image without skull only has cerebrospinal fluid, gray matter, and white matter. Then, we set with assumption (3). The bias field is smooth and slowly varying, which makes it probable that the bias field values in small local region be regarded as constant. Then the multigrid segmentation method is less sensitive to the bias field. However, each pixel is processed only in its single local grid, which makes the method very sensitive to the size of the grid and unable to preserve the statistical reliability and spatial continuity of segmentation results. This can be illustrated in Figure 1.
(a)
(b)
(c)
(d)
(e)
Figure 1 shows the segmentation results of GMM and MGMM applied on a synthetic 3 T MR image, which were created by using the MRI simulator (Brain Imaging Center at the Montreal Neurological Institute, McGill University, Montreal). There are many advantages for using these synthetic images rather than real image data volumes for validating segmentation methods. The simulator can provide full three-dimensional data volumes which have been simulated using three sequences (T1-, T2-, and PD-weighted) and a variety of slice thicknesses, noise levels, and intensity inhomogeneity levels, providing the ground truth of the image data. In this paper, the parameters of the simulated data sets are as follows: Phantom: normal; slice thickness: 1 mm; scan technique: SFLASH; TR = 18 msec; flip angle = 30 degrees; TE = 10 msec. The dimension of the image data sets is 181 × 217 × 181. The parameters of Figure 1(a) are as follows: noise level 0% and INU (RF) level 100%. Figure 1(a) shows the initial image, which is uniformly divided into 4 × 4 nonoverlapping grids. Figure 1(b) shows the results of GMM. Due to the effect of the bias field, some voxels of WM and CSF have been misclassified into GM. Figures 1(c)–1(e) are segmentation results of MGMM when the image is divided into 3 × 3, 4 × 4, and 5 × 5 grids, respectively.
It can be seen from the grids (2, 2) and (3, 2) in Figure 1(c) or from the grids (3, 2) and (4, 2) in Figure 1(d) that the variation of intensity distributions of neighbor patches would, more or less, lead to discontinuities across the grids’ boundaries. Furthermore, when the bias field is severe, grids with large size can satisfy assumption (3); however, this would make the method sensitive to the bias field and does not satisfy assumption (4). It is hard to hold assumption (3) in some grids when the size of the grids is small. For example, there are no WM pixels in the grid (3, 1) of Figure 1(e), which makes some pixels of GM misclassified into WM. It can also be found in grid (3, 2) of Figure 1(e) that tissues of WM and GM predominate the grid, which yields a deviant intensity distribution far from that typical of the brain. The GMM misclassifies some GM pixels with relatively low intensity into the class of CSF.
3.1. New Multigrid Generation Method
In this paper we present a new method to generate the multigrid. Firstly, the boundary of the brain needs to be found, because there are a large number of pixels belonging to the background in brain MR images, which usually affect the accuracy of segmentation methods. Secondly, the brain region is divided into small nonoverlapping grids. The generated nonoverlapping grids may not satisfy assumption (3). Figure 2 shows the generated multigrid on brain MR images. In this paper, we set . It can be seen from Figure 2(a) that the grids (1, 1) and (6, 3) only contain some CSF pixels and the grids (1, 6) and (6, 6) have no brain tissues. Furthermore, the grid (6, 1) has no pixels of the WM. The NGMM cannot obtain accurate results based on these grids. In order to deal with this problem, the small grids need to be combined. The combine process includes 6 steps as follows.
(a)
(b)
(c)
(d)
(e)
(f)
Step 1. Define a matrix with the same size of the multigrid and a variable to count the serial number of the final patches.
Step 2. Classify every grid into classes.
Step 3. Find the worst grid and set .
Step 4. Find the preferred grid to be combined into the patch and set . Repeat this process times.
Step 5. Set and go to Step until all grids have been labeled except those nonbrain grids.
Step 6. Combine single patches into the preferred patches nearby.
In Step , matrix is defined to sign the multigrid and initialized to zero for . The variable is used to count the serial numbers of the final patches and initialized as 1. In Step , every grid is classified into classes by using Fuzzy Clustering Means (FCM) method, which can classify grids efficiently. We assume that assumption (3) holds in all grids and all patches are classified into 3 classes except those background grids. The histograms of the grids and (3, 3) are shown in Figures 2(b) and 2(c), respectively. There are only GM and some CSF pixels in the grid , which make the distribution of the pixels’ intensities more compact than that of (3, 3).
In Step , we first calculate an inner distance for every grid. The inner distance is defined as , where , , and are the intensity centers of CSF, GM, and WM, respectively. The worst grid with the smallest inner distance can be found easily and is set as . In Step , the grids , which are next to the patch need to be found. Then, every grid is combined into the patch independently and the inner distance of the corresponding combined patch is calculated. The grid with the largest inner distance is regarded as the preferred grid to be combined into patch . This process should be repeated for times, where satisfies . is the size of the grids. is the amount of tissue pixels in the input image. is the control parameter. It is hard to hold assumption () with large . In this paper, we set , then . Furthermore the width and height of every patch are no more than 1/2 of the brain region’s width and height.
After Step , all grids have been labeled; however, there may be some grids that have not been combined with any neighbor patch. These grids are set as single grids. It is hard to hold assumption (3) in these grids. In order to deal with this problem, we use similar strategy shown in Step to combine the single grids into best neighbor patch. Firstly, we find one of the single grids. Secondly, we analyze the neighbor patches, which are next to the single grid, to find the patch with largest inner distance as the preferred patch. The single grid is combined into the preferred patch. The process is repeated until all single grids have been combined into their preferred patch.
Figures 2(d)–2(f) show the results of the new multigrid generation method on three synthetic brain MR images, with the following parameters: noise levels are 0%, 0%, and 5%, respectively, and INU (RF) levels are 0%, 100%, and 40%, respectively. It can be seen that our method can generate good multigrid even when the image has strong noise or bias field.
3.2. Multigrid Nonlocal Gaussian Mixture Model (MNGMM)
Assumption (3) can be retained by using the improved multigrid. However, due to the effect of the bias field, the variation of intensity distributions of neighbor patches would, more or less, lead to boundaries discontinuous across edges of some grids. Despite this problem, the assumptions for simplifying the image model are still correct in most cases, due to the complicated and convoluted structure of human brain as mentioned above. In other words, most grids would yield good judgments based on the expected values . A very natural idea is that the adverse impact of some grids could be overcome by all the good clustering results from other neighbor patches. Such a consideration leads to the development of a novel method called the multigrid nonlocal Gaussian mixture model (MNGMM) that can take advantage of the adaptation of local clustering but also keep the classifications spatially continuous and statistically reliable.
MNGMM includes two basic parts: clustering and information fusion. The rationale of MNGMM can be described as follows. First, all grids need to be classified by using NGMM to obtain the corresponding expected values . For every pixel , it is easy to find the corresponding grid and obtain its expected values in the grid. The distribution information of the neighbor grids next to is used to calculate , where is the distribution parameter of and is the total number of the neighbor grids. In the information fusion stage, all these expected values are integrated with a strategy to adapt the current expected values. Consequently, for each pixel , MNGMM can be summarized as follows.
Step 1. One has
Step 2. One hasHere, we name as and as ; the strategy is
In our method, only pixels on edge regions of the grid may be misclassified into other class. The intensities of the pixels are continuous across the edges of the neighbor grids, which makes it probable that the expected values of the misclassified pixels be weighted averaged by using the information of the neighbor grids. As a result, the weighted averaging values of different grids make the final result more reliable. Moreover, the weighted averaging operation can also preserve spatial continuity of the membership distributions of each tissue class. Finally, maximum expected values principle is used to obtain the segmentation results from the weighted averaging values.
4. Results
In order to show the robustness of our method, we compared our method with GMM and Wells method on a clinical 3 T MR image, which has severe intensity inhomogeneity and noise. The acquisition parameters of the real data are as follows: slice thickness 1 mm, echo delay time (TE) 35 msec, repetition time (TR) 450 msec, and flip angle 90 deg. The size of the data is 256 × 256 × 170. From the image, we can find that most intensities of putamen are higher than those of the cortex, which also belong to the GM. GMM misclassifies putamen into WM since some intensities of putamen are closer to WM than to GM, which can be shown in Figure 3(b). The result of the Wells method [6] is shown in Figure 3(c). It can be seen from the result that the method can reduce the effect of the bias field; however, with the effect of the inhomogeneities in inner tissues, some pixels of putamen have been misclassified into WM and the method is sensitive to the noise. In contrast, MNGMM can yield satisfactory results, which are more compatible with human visual perception.
(a)
(b)
(c)
(d)
In order to quantitatively show the performances of the proposed model, we compared our method with GMM, Wells method, MCFC, and MNGMM on the synthetic data from MRI simulator with the following parameters: noise level 5% and INU level 80%. The results are shown in Figure 4. The left column shows the initial image. The end column shows the ground truth. The second column to the fifth column show the segment results of GMM, Wells, MCFC, and MNGMM, respectively. Due to the effect of the intensity inhomogeneity, many pixels of WM have been misclassified into GM in the result of GMM. Wells and MCFC can reduce the effect of the intensity inhomogeneities; however, they are affected by the noise. By using the nonlocal information and multigrid information, MNGMM can yield much more complete and continuous results, which are very similar to the ground truth.
We use Jaccard similarity value (JS) [6] to quantitatively evaluate the performance of a classification method. The value of JS ranges from 0 to 1, with a higher value representing a more accurate segmentation result. The statistical results (means and standard deviations of JS values for WM, GM, and CSF) are listed in Table 1. The results demonstrate that our method produces the most accurate results and has the best ability and robustness to the noisy images (with lower standard deviations of JS values and higher mean of JS values when the noise increases), especially in the area with abundant textures (with higher JS values for CSF tissue). We also apply the above five methods to the segmentation of 40 whole synthetic MR image data sets, in which the level of intensity inhomogeneity ranges from 20% to 100%. The segmentation accuracy is measured in terms of the average JS of WM, GM, and CSF delineation and is shown in Figure 5. Both visual and quantitative comparisons show that our method is more robust to the intensity inhomogeneity and can obtain more accurate results.
(a)
(b)
(c)
Figure 6 shows the segmentation results on a real brain MR data from the Internet brain segmentation repository (IBSR at http://www.cma.mgh.harvard.edu/ibsr/) with the name 12_3 (39th image). The intensity distribution of the basal ganglia is midway between the assumed distributions of GM and WM and the basal ganglia have low contrast. From the results, we can find that our method can obtain accurate result. In order to quantitatively evaluate the benefits, we segmented 20 standard sets of real brain MR data from IBSR by using GMM, Wells method, MCFC, and MNGMM. The average quantitative results of GM, WM, and CSF are listed in Table 2. It can be seen that our method is more accurate than others.
(a)
(b)
(c)
(d)
(e)
(f)
5. Discussion
In our work, the number of grids is set as a constant. The choices of the divided number should be based on the size of the brain region and intensity inhomogeneity level. A smaller divided number will make the brain region only divided into few patches, which makes the proposed method sensitive to intensity inhomogeneity. A larger divided number will make the method less efficient. We have studied the relationship between size of patches and segmentation accuracies. Figure 7 shows the accuracies of the segmentation to simulated images with different parameters. In fact, the larger the patch, the more the data to be clustered; the greater the similarity between the intensity distribution of the patch and that of the input image is, the more reliable the clustering results are. However, assumptions 2 and 4 require smaller patches. The left column of Figure 6 shows the accuracies of the results using different when generating the multigrid on the input data with the following parameters: noise level, 0%, and intensity inhomogeneity level, 0%, 10%, 20%, 40%, 80%, and 100%, respectively. It can be seen that with the increase of intensity inhomogeneities the accuracies decrease when is small, which also means that the patch is bigger. At this time, assumptions () and () cannot be satisfied and the accuracies decrease very quickly. When , 7, 8, 9 the results have similar accuracies. We also analyzed the influence of . The right column of Figure 7 shows the accuracies with different when . It can be seen from the results that the accuracy is decreasing with the increase of when the intensity inhomogeneities increase. This is because the bigger leads to bigger patches, which means that assumptions () and () cannot hold.
(a)
(b)
6. Conclusions
In this paper, we have presented a theoretically simple approach to automatically segment 2D or 3D human brain MRI data. To reduce the effect of intensity inhomogeneity, a multigrid Gaussian mixture model has been proposed. In order to reduce the effect of noise, we improve the Gaussian mixture by using the nonlocal information, which can preserve more detailed information. Experimental results have proved that our method outperforms other segmentation methods when segmenting images with severe intensity inhomogeneities and noise.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work was supported in part by the National Nature Science Foundation of China (61572257 and 61502206); National Research Foundation of Korea under Grant NRF-2013K2 A2S2000777; Natural Science Foundation of Jiangsu Province under Grant nos. BY2014007-04, BK20150523, and BK20150923; and Open Project of Jiangsu Key Laboratory of Meteorological Observation and Information Processing, Nanjing University of Information Science and Technology (KDXS1404).