Gaussian mixture model for texture characterization with application to brain DTI images.

A Gaussian mixture model (GMM)-based classification technique is employed for a quantitative global assessment of brain tissue changes by using pixel intensities and contrast generated by b-values in diffusion tensor imaging (DTI). A hemisphere approach is also proposed. A GMM identifies the variability in the main brain tissues at a macroscopic scale rather than searching for tumours or affected areas. The asymmetries of the mixture distributions between the hemispheres could be used as a sensitive, faster tool for early diagnosis. The k-means algorithm optimizes the parameters of the mixture distributions and ensures that the global maxima of the likelihood functions are determined. This method has been illustrated using 18 sub-classes of DTI data grouped into six levels of diffusion weighting (b = 0; 250; 500; 750; 1000 and 1250 s/mm2) and three main brain tissues. These tissues belong to three subjects, i.e., healthy, multiple haemorrhage areas in the left temporal lobe and ischaemic stroke. The mixing probabilities or weights at the class level are estimated based on the sub-class-level mixing probability estimation. Furthermore, weighted Euclidean distance and multiple correlation analysis are applied to analyse the dissimilarity of mixing probabilities between hemispheres and subjects. The silhouette data evaluate the objective quality of the clustering. By using a GMM in the present study, we establish an important variability in the mixing probability associated with white matter and grey matter between the left and right hemispheres.

Introduction

Accurate and non-invasive methods capable of detecting and correcting localized affected areas of the brain are of substantial interest because there is significant variability in the location and extent of such areas. Moreover, there is substantial variability among individuals. The inter-subject variability of brain structures is apparent for normal subjects as well as for patients with various brain injuries [1], [2].

The statistical analysis presented in this paper is based on Gaussian models [3], [4], [5], [6], [7], [8], [9], [10]. A Gaussian mixture model (GMM) is a probabilistic model based on a Gaussian distribution for expressing the presence of sub-populations/sub-classes within an overall population/class without requiring the identification of the sub-class of interest (observational data) [11]. That is, a GMM learns to detect injured tissues using healthy patient data. Banfield and Raftery [4] considered both Gaussian and non-Gaussian models to specify the features of the clusters and to estimate which features likely to differ between clusters. These authors applied the proposed method to brain MRI images to identify similarly anatomical structures. Fraley and Raftery [6] systematically reviewed how finite mixture models provide a statistical framework for application in clustering, the effectiveness of a certain clustering method and the influence of outliers on cluster analysis. Paalanen et al. [7] investigated several estimation methods for GMMs, enabling an improvement in the representation and discrimination of patterns. Kim et al. [9] proposed a method to assess gross brain abnormalities using a GMM, with at least two Gaussian components to allocate a specific mixing probability to each subject. Then, the assigned mixture probabilities are tested between the studied groups. The authors stated that a GMM is an effective method in terms of computing resources because it does not incorporate any subject or group-specific parameters. Another popular approach is the multivariate Gaussian method, which is used to identify features and discriminate between different classes in various applications, such as hazardous chemical agents [12], moving parts of electric motors under normal conditions and those with bearing failure [13]. However, the usage of a single feature results in a single fundamental class with variables exhibiting smooth behaviour. Consequently, certain errors in the estimation of the probability density function (pdf) occur, and the discrimination between classes fails. A GMM treats analysed data as a mixture of component distributions, and the main challenge here is to correctly estimate the model parameters such as the weight of the ith component, which can be interpreted as the a priori probability, mean or covariance matrix of the normally distributed random variable. A GMM is a simple physical and data-driven model; it permits a flexible characterization of unusual distributions of pixels and provides a quantitative analysis of DTI data [14]. DTI captures vital information and plays a significant role in in vivo studies of anatomical structures in brain regions [15]. Generally, brain tissues such as grey matter (GM), white matter (WM) and cerebrospinal fluid (CSF) are considered classes; a few of their pixel features are predominantly measurable and serve as model parameters. A GMM assigns a probability to each pixel if it belongs to only one class. These parameters are learned using the expectation-maximization (EM) algorithm [16], [17]. More recently, various associations between prior knowledge on human neuroanatomy and conditional probability distributions that can characterize various brain tissues or anatomy classes were made feasible by considering a GMM associated with various neuroanatomical labels [18], [19].

In this study, a GMM-based classification scheme to identify the variability in the main brain tissue in DTI images (rather than searching for tumour areas) is proposed. To consider the signal intensity and contrast effects on image quality, we used multiple b-values. Han et al. [20] reported that brain imaging using a high b-value likely improves both the contrast between tissues and the capability of detecting less prominent lesions. As a probabilistic model, a GMM can characterize arbitrary mixture distributions composed of WM, GM and CSF with unknown parameters. A k-means algorithm is used to optimize the means, variances and mixture probability of the mixture distributions and to ensure that a global maximum of the likelihood function is achieved. The weighted Euclidean distance (wd) is used to validate the capability of a GMM to discriminate between the mixing probabilities across the studied classes. Moreover, a multiple correlation analysis between the left and right hemispheres based on the established mixing probabilities is performed. Finally, the silhouette plot (size and width) is used to evaluate the clustering validity. DTI images of a normal subject without a history of head injury or cerebrovascular disease (denoted H), a patient with multiple haemorrhage areas in the left temporal lobe (HA) and another with ischaemic stroke (IS) were studied.

Methodology

Finite GMM with m components

Brain DTI images contain heterogeneous sub-classes, and the analysis based on the mixture of models is adequate to model the entire distribution containing numerous sub-classes. Different tissues such as WM, GM, and CSF or lesioned tissues aggregate their intensities and contrast, which is essentially decided by the b-values in DTI, under different Gaussian curves with distinct mean and covariance parameters. The highest probability of classifying each pixel as belonging to the WM, GM and CSF or lesioned tissues is the basis of a GMM. The mean values of the weights of each brain tissue and for each b-value are projected onto a vector space with a three-dimensional feature [w1, w2, and w3]. This vector contains the intensities of the pixels for each available b-value. First, healthy subject-specific information is integrated into the algorithm using the GMM by means of a training stage. During this training stage, the weights, mean, and variance for each individual Gaussian density are determined. These factors are prior probabilities for parameters required in the initialization step of the EM algorithm. Then, the GMM parameters for the test data are established using the maximum likelihood and an iterative EM algorithm. The GMM generates a vector space, which contains the probability function of the data computed using the intensities of the pixels and their discriminative weights or mixture probability, for each available b-value and for each considered tissue class. The vector space is passed to the predictive model to capture the discriminative subject-specific information regarding the brain injuries from MRI images.

This approach is followed in an inter-hemisphere analysis, and the mixture weight functions are determined as a posteriori probability to prevent the repetition of an excessively large number of univariate analyses during the characterization of each subject. The brain segmentation into hemispheres involves the following steps: (i) skull stripping based on an irrational mask for filtration and binary morphological operations [21]; and (ii) mid-sagittal axis detection based on the location of the inter-hemispheric fissure and determination of the image centroid, as reported in [22].

In a preliminary step, an image histogram that provides raw information about the pdf of the pixel values is analysed. The number of components or Gaussian sources is established at three, according to the multimodal distribution in the histograms of the pixel distribution. The finite mixture model is based on the assumption that each finite mixture has similar probability distributions for each sub-class; however, inside the sub-class, different multivariate probability density distributions and different parameters are present [3]. For an image, let X denote a vector of pixel intensities . X is a feature vector of the observation data for a specific subject and a specific b-value. This vector describes a sub-class, which, in turn, belongs to a class . The probability function at the observation is expressed as

denotes a component of the Gaussian mixture or ‘mixture distribution’, and is the prior distribution of the pixel for each sub-class corresponding to b-values and is called mixing proportions or weights. The weights must satisfy the following conditions to be valid: and , for each sub-class. As an a priori distribution, it is obtained by observing a healthy patient. Each density function of the mixture component is characterized by a mean and a variance and is a univariate normal pdf expressed as

denotes a set of the main parameters of the GMM to be estimated by the EM algorithm. The likelihood function of the training vectors based on the probability function (2) is as follows:where N is the total number of pixels in the image. The function permits one to establish the statistical model parameters. There is incomplete data in ; therefore, these partially observed parameters in are estimated using an EM algorithm [23]. As an iterative algorithm, EM starts by using an initial model , estimates a new model Ω′, and in the next iteration, this new model Ω′ becomes an initial model to determine a new model, Ω″, etc. This process is repeated until a predefined convergence condition is achieved. In the studied problem, the mixing proportions are estimated and subsequently validated.

In the initialization step, the mean, variance and mixing coefficients of the training data were estimated for m = 3 classes. Established during the training stage, these parameters are the prior distribution. This distribution is used to initialize the value of the probability.

In the expectation step, based on the current estimated parameters , the EM algorithm computes the posterior probability using the current parameter values established in the initialization step; furthermore, the algorithm estimates the posterior probability that an observation belongs to a sub-class j over all feasible assignments of data points to Gaussian sources, as

Then, for each sub-class j, the prior probabilities are computed by averaging the posterior probabilities for each class as

In the maximization step, the values of the old estimated parameters are updated or re-estimated by computing the maximum likelihood estimates of with the expected membership values [24], [25]. The derivative of function (3) is determined and equated to zero. The determination of the global maxima of the likelihood functions implies the determination of those parameters that maximize the probability of observing the data. The new mean, variance and weight parameters are estimated, and the likelihood function is evaluated.

The iterative process continues through the expectation and maximization steps successively until convergence. Convergence implies that the changes in the parameters become small enough, i.e., it stops when , where  = 0.00001 is small enough to assure no significant changes from one iteration to the next exist. If the convergence criterion is not attained, the algorithm returns to the expectation step.

k-means algorithm for clustering

The k-means algorithm is used to assess the data clustering for the selected number of classes (m = 3) [26], [27]. Each mixture component is associated with a class or cluster based on identical estimated statistical parameters. The data produced by a GMM are clusters with centroids at the means. In the GMM, the EM algorithm converges to the local maxima of the likelihood function. However, the EM algorithm has a drawback, namely, it fails when the covariance matrix associated with a sub-class is singular or the number of observations is reduced. According to the GMM, three initial centroids are defined based on the set of parameters established in the maximization step of the EM algorithm. Thirty-five restarts of the k-means algorithm were executed to ensure that a global maximum is determined for each data set. That is, k-means optimizes the mean, variance and weight of the mixture distributions.

Weighted Euclidean distance and multiple correlation

To validate the ability of the GMM mixtures to differentiate between different subjects in a hemisphere approach, we used the wd between two j-dimensional vectors [28]:where is given by Eq. (4) for the studied sub-classes, and is the corresponding standard deviation. H denotes the normal subject, HA a patient with multiple haemorrhage areas in the left temporal lobe and IS another with ischaemic stroke. The wd balances the contributions of the variables in the computation of distance. The weight attached to the jth variable in a vector is related to the standard deviation of each distribution and [29]. A continuous image belongs to a metric space that uses metrics exhibiting the following properties: non-negativity, identity of indiscernibles, symmetry and triangle inequality. The wd is appropriate for measuring the dissimilarity of the two given mixing probabilities because other metrics fail to exhibit a few of these properties, e.g., the Kullback–Leibler distance fails to exhibit the symmetry property, and the Hellinger distance fails to exhibit the triangle inequality. The Minkowski and Mahalanobis distances are general formulations of the Euclidean distance. Hershey and Olsen [30] reported that the Kullback–Leibler divergence metric between two GMMs cannot be analytically applicable and that the algorithm is highly time consuming. Durrieu et al. [31] demonstrated that the Kullback–Leibler divergence metric can be approximated only.

Moreover, an inter-hemisphere multiple correlation analysis between the mixing probabilities was performed to characterize the association of the grey level intensities and contrast for the selected injured subjects and healthy subjects. The multiple correlation coefficients between the independent variables HA and IS and the dependent variable H are defined aswhere , m = 1, 2, 3 are the covariance between the two random variables in each of the pairs IS and H, HA and H and IS and HA, respectively [32]. Accordingly, the correlation coefficient between two sub-classes j that belong to a class m for two random variables is illustrated in the example below:where and are defined in the expectation step, and K = 10 is the number of samples (i.e., images) for each sub-class.

Clustering validation

The analysis is focused on three main brain tissues (i.e., GM, WM and CSF), and an a priori assumption of three-class clustering is considered. The goal is to examine whether these classes reflect the actual clustering structure of the data or whether these data were partitioned into a few artificial groups in the context of the GMM [33]. The quality of the clustering analysis is addressed, and the silhouette index and silhouette plots are used as the validation criteria [34]. If compact and clearly separated clusters are obtained, the targeted tissues were considered well classified. Let a multivariate data be separated into m clusters, A, , . Let us suppose that the GM tissue is described by and . We define the average dissimilarity of with all the other points k of the same cluster having the vector norm , as follows:where || || denotes a 2-norm (L2). Further, and describe the average dissimilarity of the mixing probability of H with all the points belonging to other clusters HA and IS, respectively:

The smallest average dissimilarity to another cluster is defined as . The silhouette index is

From Eq. (10), , and if , the least effective situation manifests. This method is also used for WM and CSF. Silhouette plots facilitate the interpretation of cluster analysis results because they are independent of the clustering algorithm used and rely only on the actual partition of the ‘objects’ [34].

Subjects, image acquisition, and processing

The algorithm flow is presented in Fig. 1.

Fig. 1

Algorithm scheme.

Three subjects (age range 36–60 y; one female and two males) underwent MRI scans. A subject presented multiple haemorrhage areas in the left temporal lobe (male, 48 y), and another presented with IS in the left frontal lobe (female, 60 y, median 8-mo post-stroke). The third subject was a healthy patient (male, 36 y). A series of DTI images were acquired using a pulsed gradient spin-echo sequence in 15 directions and five b-values (b1 = 250 s/mm2; b2 = 500 s/mm2; b3 = 750 s/mm2; b4 = 1000 s/mm2; b5 = 1250 s/mm2). Moreover, images without diffusion gradients (b0 = 0 s/mm2) and with otherwise identical imaging parameters were acquired. A total of 190 images were tested. A b-value encompasses information regarding the strength and timing of the gradients used to generate diffusion-weighted images. Larger b-values provide better contrast among tissues. The selection of b-value continues to be a challenge and strongly depends on the investigated anatomical features or pathology, field strength and average number of signals. In the case of the GMM, the mixing probabilities depend on the experimental conditions, i.e., the diffusion effect or b-value. Multiple b-values permit the use of a small sample size because each data set exhibits characteristics unique to it. The within-subject correlation is avoided by summarizing each mixing probability sequence with a single number. In this case, only a comparison of the statistics between the classes (see data in Table 3) is performed. Averaging repeated measurements is a reasonable choice, especially when the effect of the injury is maintained quite steadily over acquisition time.

Table 3

Correlation coefficients and multiple correlation coefficients.

	Correlation coefficient									Multiple correlation coefficient
	rHA,H1	rHA,H2	rHA,H3	rIS,H1	rIS,H2	rIS,H3	rIS,HA1	rIS,HA2	rIS,HA3	RH-(IS,HA)1	RH-(IS,HA)2	RH-(IS,HA)3
Left hemisphere	0.658	−0.421	0.214	0.654	0.515	−0.612	0.214	−0.031	0.295	0.528	0.429	0.545
Right hemisphere	0.751	−0.773	0.654	0.443	0.339	0.622	0.336	−0.214	0.345	0.714	0.699	0.564

For the data acquisition, a 1.5-T MRI scanner was operated (Philips Medical Systems, Best, Netherlands). The diffusion-weighted scans utilized a system with six-channel sensitivity encoding (SENSE) for faster scanning (FS = 1.5). The imaging parameters were as follows: 3D gradient echo with echo time ranging from 83 to 110 ms; repetition time ranging from 6500 to 7800 ms (it varies between subjects); bandwidth = 1070 Hz/pixel; flip angles (2- and 6-); voxel resolution ranging from 2.5 to 3.0 mm; and slice thickness = 4 mm. The acquisition matrix was 128 × 128. The standard Digital Imaging and Communications in Medicine (DICOM) image dataset was used.

Approval for the study was obtained from the Research Ethics Committee of the Dunarea de Jos University of Galati and Saint Andrew Hospital. Voluntary and written informed consent was obtained from each participant. The privacy policy is based on DICOM Confidential [35].

Results

The proposed GMM-based classification approach in a brain hemisphere framework is aimed at facilitating the identification of variability in the main brain tissue in DTI images and circumventing the subsequent processing for detecting tumours or lesions. For example, a DTI image (b = 500 s/mm2) of a healthy subject and the results of the GMM classification and hemisphere segmentation are shown in Fig. 2.

Fig. 2

DTI brain image of a healthy patient for b = 500 s/mm2. (a) Skull stripping of the whole brain; (b) result of GMM classification; (c) hemisphere segmentation.

The estimated weights (Eq. (5)) across the entire control group (H) and for each injured group (IS and HA) are presented in Table 1 (for the left hemisphere) and Table 2 (for the right hemisphere). The data in Table 1, Table 2 present details on the difference in the averaged weights or mixing probabilities between the left and right hemispheres for each subject and over the entire range of diffusion gradient values. There are no differences in mixing probabilities for each tissue class between the left and right hemispheres for H class. This result indicates the ‘normality’ of the healthy subject. For HA and IS, visible differences in the mixing probabilities are presented.

Table 1

GMM average mixing probability for the left hemisphere with and without diffusion gradients. The data are summarized for three mixing probabilities (w1 for GM, w2 for WM and w3 for CSF) and for three subjects H, HA and IS.

	w_1±sjH	w_1±sjHA	w_1±sjIS	w_2±sjH	w_2±sjHA	w_2±sjIS	w_3±sjH	w_3±sjHA	w_3±sjIS
b0	0.30 ± 0.012	0.26 ± 0.008	0.22 ± 0.014	0.55 ± 0.049	0.58 ± 0.048	0.55 ± 0.064	0.14 ± 0.045	0.15 ± 0.045	0.17 ± 0.061
b250	0.32 ± 0.021	0.30 ± 0.017	0.29 ± 0.020	0.52 ± 0.046	0.52 ± 0.051	0.54 ± 0.059	0.14 ± 0.051	0.15 ± 0.051	0.15 ± 0.0.047
b500	0.32 ± 0.027	0.30 ± 0.020	0.29 ± 0.021	0.52 ± 0.056	0.54 ± 0.057	0.55 ± 0.048	0.14 ± 0.045	0.15 ± 0.047	0.15 ± 0.051
b750	0.33 ± 0.023	0.29 ± 0.023	0.29 ± 0.021	0.54 ± 0.015	0.55 ± 0.052	0.56 ± 0.053	0.13 ± 0.007	0.15 ± 0.042	0.15 ± 0.041
b1000	0.32 ± 0.021	0.29 ± 0.020	0.29 ± 0.019	0.52 ± 0.048	0.55 ± 0.053	0.55 ± 0.058	0.14 ± 0.050	0.15 ± 0.052	0.15 ± 0.048
b1250	0.33 ± 0.025	0.28 ± 0.018	0.28 ± 0.016	0.55 ± 0.059	0.55 ± 0.055	0.55 ± 0.056	0.14 ± 0.044	0.15 ± 0.050	0.15 ± 0.049

Table 2

GMM average mixing probability for the right hemisphere with and without diffusion gradients. The data are summarized for three mixing probabilities (w1 for GM, w2 for WM and w3 for CSF) and for three subjects H, HA and IS.

	w_1±sjH	w_1±sjHA	w_1±sjIS	w_2±sjH	w_2±sjHA	w_2±sjIS	w_3±sjH	w_3±sjHA	w_3±sjIS
b0	0.32 ± 0.016	0.33 ± 0.018	0.28 ± 0.013	0.54 ± 0.041	0.53 ± 0.050	0.57 ± 0.039	0.12 ± 0.032	0.13 ± 0.041	0.14 ± 0.048
b250	0.34 ± 0.022	0.33 ± 0.041	0.30 ± 0.021	0.51 ± 0.053	0.53 ± 0.051	0.53 ± 0.059	0.14 ± 0.045	0.13 ± 0.044	0.15 ± 0.047
b500	0.34 ± 0.022	0.34 ± 0.029	0.30 ± 0.028	0.50 ± 0.050	0.51 ± 0.059	0.52 ± 0.055	0.14 ± 0.045	0.14 ± 0.044	0.15 ± 0.050
b750	0.35 ± 0.016	0.34 ± 0.029	0.28 ± 0.023	0.52 ± 0.007	0.48 ± 0.042	0.52 ± 0.041	0.13 ± 0.039	0.13 ± 0.039	0.16 ± 0.043
b1000	0.33 ± 0.018	0.33 ± 0.029	0.30 ± 0.029	0.51 ± 0.054	0.51 ± 0.060	0.53 ± 0.055	0.14 ± 0.044	0.14 ± 0.045	0.15 ± 0.047
b1250	0.34 ± 0.019	0.33 ± 0.026	0.30 ± 0.023	0.54 ± 0.053	0.51 ± 0.057	0.51 ± 0.061	0.14 ± 0.044	0.14 ± 0.045	0.15 ± 0.046

In Fig. 3, the estimated average wds (Eq. (6)) for all diffusion gradients and for each brain hemisphere and subject are presented.

Fig. 3

Average weighted Euclidean distances for pairs of probability density function distributions of mixtures probability of GMM. Estimation is performed for all diffusion gradients and for each brain hemisphere. L denotes the left hemisphere, and R denotes the right hemisphere.

Fig. 3 indicates that the proposed approach exhibits the ability to highlight the differences between brain tissues in the right and left hemispheres for each level of diffusion weighting and subject category. This distance balances the contributions of the variables by considering the standard deviation of each distribution.

The correlation matches images characterized by various intensities and contrasts, albeit with largely similar local intensity variations. The results of the correlation analysis (Eqs. (7), (8)) are presented in Table 3. First, the correlation between each pair of classes has been investigated. The results indicate that classes HA and IS are not correlated because the correlation coefficient is near zero. This observation leads to the following hypothesis: H is the dependent variable, and HA and IS are not correlated and are the independent variables. Therefore, the multiple correlation coefficient is computed according to Eq. (7).

As the data in Table 3 indicate, for the CSF class (index 3), HA and IS do not correlate with H for neither the left or right hemisphere. The results for the WM class (index 2) illustrate that, for the left hemisphere, HA and IS are marginally correlated with H. The correlation increases by approximately 50% for the right hemisphere. For the GM class (index 1), HA and IS correlate well with H for the right hemisphere and do not correlate with H for the left hemisphere.

The resulting silhouette plots (Eq. (10)) for the whole brain and the left and right hemispheres are displayed in Fig. 4.

Fig. 4

Silhouettes of a data set for three clusters (line 1 on the silhouette plot corresponds to healthy subjects, line 2 for HA and line 3 for IS). Row 1: whole brain; Row 2: right hemisphere; Row 3: left hemisphere.

The average silhouette width is approximately 0.9, i.e., 90% of the selected clusters are considered the optimal number of clusters (Table 4). The a priori selection of the three main brain tissues or ‘natural determination’ is validated and performs best with respect to the hemisphere approach. The width of cluster 2 (HA subject) is not significantly high for CSF and GM in the left hemisphere. This narrow silhouette is interpreted as a spread of the point inside the cluster and as a slightly inadequate separation of the cluster.

Table 4

Average silhouette width for evaluating clustering validity.

Class	Whole brain	Left hemisphere	Right hemisphere
H	0.9176	0.935	0.9284
HA	0.9829	0.9774	0.9326
IS	0.9989	0.8578	0.9296

Discussion

A different classification scheme based on the GMM for identifying the variability in the main brain tissue through a hemisphere approach (rather than by searching for tumours or lesion areas) was presented. A whole-brain imaging analysis is labour intensive and tends to be biased towards structural anatomical boundaries. Brain asymmetry analysis is a tool for analysing the neuroanatomical basis of disorders with an assumed developmental aetiology, such as dyslexia, autism and schizophrenia, in men and women [36], [37], [38], [39]. Most studies have focused on exploring the asymmetry of the WM structure. Furthermore, these studies are mostly based on a region of interest in an image data set that is specified by users. Our results follow these observations and enlarge the applicability of the research of hemispheric specialization to the apparent difference in statistical features to reveal abnormal asymmetries of the statistical distribution of the main brain tissues. By using multiple b-values, we constructed a tool to evaluate Gaussian diffusion based on the decreased degree of diffusion-related signal attenuation with the increased b-value.

Mixture distribution models such as a GMM expresses the presence of the sub-class in a class without requiring that the sub-class of interest (observational data) be identified [11]. That is, a GMM expresses the probability distribution of the observational data in a class. We are focused on the three main brain tissues; thus, a GMM extracts a class's characteristic from a sub-class. As a mixture distribution model, a GMM does not seek the sub-class’s information identification; since a GMM can simultaneously provide the observational data about the class, it also provides a statistical inference about the characteristic of the sub-class. Generally, a GMM requires the number of components to be specified in advance for analysing the data, i.e., inputting the number of components m (Eq. (1)) present in the mixture is necessary [40]. For ten classes of univariate distributions (including Gaussian distributions), Khalafal–Hussaini and Ahmad [41] established that all the finite mixtures generated by the family of parameters are identifiable. Chen et al. [42] reported that a finite mixture model with k components (k = 2 and k ≥ 3) appears to provide consistent results in both cases (k = 2 and k = 3). Moreover, these authors claimed the absence of evidence that indicates k ≥ 4.

The images contain multiple regions with different intensity distribution characteristics. Pixels with similar characteristics will cluster together. However, pixel classification as either CSF, GM, or WM can have a < 100% probability of belonging to a certain brain tissue. In this case, a low mixing probability can be interpreted as a possibility that a pixel has lower percentages of content of the various tissues, as data in Table 1, Table 2 showed for GM and CSF.

Specifying that the pdfs were estimated and that the mixing probabilities were computed based on the individual pixel distribution is necessary. Therefore, the highlighted differences in the probability density function originate from the particular feature of each mixing probability. The GMM analysis through the hemisphere approach evidently indicates the ‘normality’ of the healthy subject. There is no difference in the mixing probabilities between the left and right hemispheres for any of the classes. In contrast, a GMM with mixture probabilities tested between the left and right hemispheres for the injured subjects (both HA and IS) indicated differences and permitted the estimation of the effect of disease on the pixel distribution. A higher difference is captured for the w2 mixing probability characteristic of WM, according to the restricted diffusion mechanisms. The difference in the diffusion coefficient between normal WM water diffusion and diseased tissue indicates the loss of WM microstructural elements. w1 exhibits smaller albeit visible differences. The proposed approach transcends limitations identified by Schmithorst et al. [43], which have determined developmental differences between males and females in the brain structure or various pathologies such as alcoholism [44] and schizophrenia [45].

To validate these observations, we computed the wd based on the mixing probability between pairs of subjects. One of each pair’s members is a healthy subject (H), whereas the other is a patient with one of the studied diseases (HA and IS). Fig. 3 shows the results of significant variation in the wd values. An evident difference between the brain hemispheres is apparent when a paired comparison is performed between the parameters of the healthy subject and the patients with HA and IS diseases. The differentiation of mixing probabilities for the left and right hemispheres by wd provides a simple tool for assessing the variability in the main brain tissue in DTI images. To cross-verify this hypothesis, we performed a multiple correlation analysis. This analysis highlighted that HA and IS do not correlate and are the independent variables and that H is the dependent variable. The multiple correlation analysis reinforced the conclusion that for an uninjured right hemisphere, the mixing weights for HA and IS correlated well with those for H, and for the left hemisphere, the correlation weakened and indicated brain injuries. A weak correlation between the injured subjects and the control for the left hemisphere validates the variability of mixing probabilities inside the class. Furthermore, because we have a priori grouping of our data into three clusters, the silhouette plots graphically validate that the analysed ‘objects’ are grouped into three natural clusters. The explanation for the less wide silhouette of cluster 2 (HA subject) for the CSF and GM in the left hemisphere (Fig. 4) lies in the spread of the points inside the cluster. The data reported in Table 4 validate the initially assumed hypothesis for the three classes used in the GMM.

To our knowledge, no studies have been reported on a GMM based on pixel intensities and contrast and following a hemisphere approach to assist with brain injury diagnosis. A GMM approach with fMRI data has been proposed to explore hemispheric lateralization for language production or the human visual system in genome-wide association [46], [47], [48], [49]. The hemispheric functional lateralization index probability density function was modelled separately for both the hemispheres using a mixture of n Gaussian components with fMRI data [47]. Furthermore, recent developments used an alternative thresholding approach based on model fit as part of mixture distribution to demonstrate that mixture modelling provides satisfactory results for the human visual system [48]. A study carried out by Kherif and Muller [49] on subjects with aphasia caused by IS demonstrated that GMMs are capable of dissociating between the sub-groups of the subject based on the main sources of variability in fMRI (i.e., handedness, sex, and age). Moreover, the authors reported that the GMM in combination with fMRI and automated lesion detection techniques is a reliable method for analysing how a normal language function is sustained notwithstanding brain injuries in the critical area. A recent study performed by Pepe et al. [50] investigated the local statistical shape analysis of gross cerebral hemispheric surface asymmetries through the brain’s morphological features (i.e., surface vertices) to establish the correspondence between the hemispheric surfaces. The proposed statistical method was tested on a small sample of healthy patients and first-episode neuroleptic-naïve patients with schizophrenia.

A few limitations of the proposed approach are the following: (i) an important topic for further studies is the monotone change of relative pixel numbers with the age of the patients. A decrease of the relative number of pixels from the brain tissues as the age of the patient increases exists and age-related changes were found in the mean and variance. GMM can be affected by this finding and further examinations with more extensive age classes are required and (ii) in the current study, the optimal number of Gaussian components (m) for GMM was determined based on histogram distribution. Other methods like the Akaike Information Criterion or Bayesian Information Criterion can be used to determine this number.

Employing a GMM provides flexibility in terms of pixel spatial distributions that can be associated with a specific pathology. This method may be used to automatically detect brain microstructural differences, which exhibit statistical characteristics different from those for the same hemisphere in a normal subject, when mixing weights are considered. A major advantage is that the statistical approach over hemispheres accurately identifies the structural variations in the brain tissues by using a small number of data samples to estimate the GMM parameters. Another advantage of using the GMM for this application is that it is based on unsupervised learning; in addition, it can be rapid and to a certain extent, capable of circumventing the subsequent processing for detecting tumours or lesions. Moreover, the proposed approach is unbiased, not operator dependent and circumvents the region of interest. The main drawbacks of a DTI acquisition system (such as noise, vibration and movement artefacts) can be overcome by this multimodal approach.

Conclusions

This study, based on the asymmetries of mixture distributions between the left and right hemispheres in the human brain, can improve and more effectively assist in early diagnosis. This study is a collection of cross-sectional data samples of different subjects. The main advantage of this approach is that it is very simple, fast and can summarize the existing differences between subjects. An important source of variability in the probability density function distribution for the w2 (associated with the WM) and w1 (associated with the GM) mixing probabilities between the left and right hemispheres was established. The differences between the subjects in terms of mixing probabilities were also reflected by the variation in the wds. The GMM approach, mixing probabilities and wd measure represent practical and convenient tools for large-scale meta-analysis of DTI data without searching for delimitation of tumour/affected areas. Specifically, two advantages were identified. The statistical approach over the hemispheres accurately identifies structural variations in brain tissues by using a small number of data samples to estimate the GMM parameters; moreover, it is unbiased, operator independent and circumvents the region of interest.

Conflict of interest

The authors have declared no conflict of interest.

Compliance with Ethics requirements

All procedures followed were in accordance with the ethical standards of the responsible committee on human experimentation (institutional and national) and with the Helsinki Declaration of 1975, as revised in 2008 (5). Informed consent was obtained from all patients for being included in the study.