Brain MRI data is typically stored in a three-dimensional format. In this study, gray matter sMRI images (96×113×94 voxels) were analyzed. If features are extracted from each voxel a dimensional disaster will occur as model performance will be greatly reduced by the extraction of a large number of irrelevant or redundant features. Therefore, here a preprocessing method that includes the slicing and calculation of weighted sums of average gray images is proposed. The detailed steps include:

(1) Slicing original image: For each subject, the size of the gray matter image was $ 96 \times 113 \times 94$ voxels. And the volume image was sliced in the $ {Z} $-axis to give 94 slices.

(2) Selecting and converting to gray image: By removing 10 slices (the first and last five slices in the image stack) which did not include feature information, and converting the remaining slices into gray images, sequentially numbered slices were obtained and denoted as $ a_{i} ~(i = 0, 1, 2, \dots, 83) $. A part of the sliced and grayed image slices from the first subject numbered NC001 are shown in Fig. 1.

(3) Weighting and averaging the gray image: Different weights are given according to the information contained in each slice. Here, the ordered slices are divided into three groups, and the weighted average gray image (Img1, Img2, and Img3) is calculated for each group. The calculations are:

(1) $ Img1=\frac{m_{0} \times 1+m_{1} \times 2+m_{2} \times 3+\ldots. +m_{27} \times 28}{1+2+3+\ldots +28},\\ $ (2) $Img2=\frac{m_{28} \times 1+m_{29} \times 1+m_{30} \times 1+\ldots. +m_{55} \times 1}{28}, \\$ (3) $ Img3=\frac{m_{29} \times 28+m_{30} \times 27+m_{31} \times 26+\ldots. +m_{83} \times 1}{28+27+26+\ldots +2+1},$

The sMRI data from each subject gives three image sets (Img1, Img2, and Img3) after preprocessing to give 600 preprocessed images from 200 subjects, which provides a data set suited for subsequent feature extraction.

Image feature extraction is a fundamental and critical step in medical image processing. The purpose is to obtain the characteristics or attributes of samples as numerical values, symbols, and feature vectors. The resulting feature extraction directly affects classification accuracy. Because texture information in the image is not sensitive to noise, light, and color, it is this feature that is chosen for the analysis.

2.2.1 Texture features based on the gray-level co-occurrence matrix

The gray-level co-occurrence matrix (GLCM) is a statistical matrix that describes the grayscale of adjacent pixels (or within a certain distance) [17, 18]. Assume the gray level of a digital image is $N$ , $p\left(q,w\right)$ represents the possibility (or frequency) of the appearance of grayscale $w$ under the condition that the starting grayscale is $q$ , where it is assumed that $w$ is along the direction $\theta$ of $q$ and the space distance is $d$ . $\theta$ denotes the angle between the position whose grayscale is $q$ and the position whose grayscale is $w$ . A GLCM contains statistical information which reflects the gray direction, interval, and amplitude variation of an image. It is calculated as:

(1) Mean:

(4) $ Mean=\bar{x}=\sum \limits_{q=0}^{N} \sum \limits_{w=0}^{N} p\left(q, w\right)\times q, $

The mean reflects the regularity of texture. The smaller the mean, the more disorganized the texture.

(2) Variance:

(5) $Variance=\sum \limits_{q=0}^{N} \sum \limits_{w=0}^{N} p\left(q, w\right)\times \left(q-\bar{x}\right)^{2},$

Variance measures deviation of the pixel value from the mean. The larger the variance, the greater the change of gray scale.

(3) Entropy:

(6) $Entropy=-\sum \limits_{q=0}^{N} \sum \limits_{w=0}^{N} p\left(q, w\right)\times \ln p\left(q, w\right),$

Entropy gives a measure of the information contained in an image. The greater the entropy, the more complex the texture.

(4) Contrast:

(7) $ Contrast=\sum \limits_{q=0}^{N} \sum \limits_{w=0}^{N} p\left(q, w\right)\times \left(q-w\right)^{2},$

Contrast reflects the total amount of local gray scale changes in an image. The greater the contrast of an image, the clearer the visual effect of an image.

(5) Correlation:

(8) $ Correlation=\sum \limits_{q=0}^{N} \sum \limits_{w=0}^{N} \frac{\left(q-Mean\right)\times \left(w-Mean\right)\times p\left(q, w\right)^{2} }{Variance},$

Correlation is a measure of the linear relationship of the gray scale. The longer the extension of the gray value in a certain direction, the greater the correlation.

(6) Homogeneity:

(9) $ Homogeneity=\sum \limits_{q=0}^{N} \sum \limits_{w=0}^{N} p\left(q, w\right)\times \frac{1}{1+\left(q-w\right)^{2} },$

Homogeneity is a measure the uniformity of the local gray level of an image. The more homogeneous the local gray scale, the greater the homogeneity value.

(7) Energy:

(10) $ Energy=\sum \limits_{q=0}^{N} \sum \limits_{w=0}^{N} p\left(q, w\right)^{2},$

Energy is the measurement of uniformity of gray distribution in the image.

A texture vector representing an image can be obtained by use of the foregoing seven textural features.

2.2.2 Normalization

Feature vectors must be normalized so that no particular feature dominates all the others. Feature vectors are normalized to zero mean and unit variance. The normalization is performed as:

(11) ${x\mathrm{\text{'}}}_i=\frac{x_i-\mu }{{\sigma }_x},$

where $\mu$ and ${\sigma }_x$ are respectively the mean and standard deviation of all the features $x_i$ .

In summary, for each of the three weighted and averaged gray images (Img1, Img2, Img3), the above seven statistics are calculated and normalized. Thus, each subject can be represented by a vector containing 21 features.

Machine learning can often train a model by use of the given data and can perform tasks such as classification, recognition, and segmentation. Machine learning algorithms can extract knowledge from data and make predictions using data. Here, k-nearest neighbor (KNN), SVM, logistic regression (LR), gradient boosting (GB), and an improved eXtreme gradient boosting (XGBoost) classifier are adopted as the methods by which to compare for medical image processing techniques.

The boosting method [19, 20], which combines the additive model (linear combination of basis functions) and the forward stepwise algorithm, is a widely used and effective statistical method for optimizing learning. A boosting method based on a decision tree is called a boosting tree. The GB model is established in the gradient descent direction of loss function which is derived from the above model of boosting method.

A GB classifier is used to train and classify the feature vector of the sample, and the complete set of algorithmic steps is:

(1) Determine the training and testing sets. 80% of the 200 samples were used as training sets, i.e. $T=\left(x_1,y_1\right),\left(x_2,y_2\right),\dots (x_N,y_N),$ $x_i\in X\subseteq R^n,y_i\in \left\{0,+1\right\}$ . The remaining samples were used as a testing set. Let $x_i$ denote the feature vector of each sample, $y_i$ denotes its class label. When the sample is either SCZ or NC, $y_i$ is 0 or +1, respectively.

(2) Confirm the loss function:

(12) $L\left(y,f\left(x\right)\right)={\left[y-f\left(x\right)\right]}^2,$

where $f\left(x\right)$ is the fitting function of $y$ .

(3) Initialize model variables:

(13) $ f_{o} \left(x\right)=\arg\min_{c} \sum \limits_{i=1}^{N} L\left(y_{i}, c\right),$

where $c$ is the constant that minimizes the loss function and represents a tree with only one root node.

(4) For each model do the following:

(a) For the $i$ th $(i=\mathrm{1,2},\dots ,N)$ sample, calculate the negative gradient of the loss function in the current model:

(14) $r_{\mathit{mi}}=-{\left[\frac{\partial L\left(y,f\left(x_i\right)\right)}{\partial f\left(x_i\right)}\right]}_{f\left(x\right)=f_{m-1}\left(x\right)},$

where $m$ is the model number, $m=\mathrm{1,2},\dots ,M$ . $M$ denotes the maximum value of $m$ . $J$ is the maximum value of $j$ , and j is the number of leaf node area.

(b) Fit a regression tree for $r_{\mathit{mi}}$ and obtain the leaf node area $R_{\mathit{mj}}$ ( $j=\mathrm{1,2},\dots ,J$ ) of the $m$ th tree.

(c) The number of leaf nodes $j(j=\mathrm{1,2},\dots ,J)$ is calculated from:

(15) $c_{mj} =\arg\min_{c} {\sum \limits_{x_{i} \in R_{mj} }} L\left(y_{i}, f_{m-1} \left(x_{i} \right)+c\right),$

In this step, linear search is used to estimate the value of the leaf node region and to minimize the loss function.

(d) Update regression tree:

(16) $f_m(x)=f_{m-1}(x)+\sum_{j=1}^J c_{m j} I\left(x \in R_{m j}\right)$

where if $x\in R_{\mathit{mj}}$ , $I=1$ , or $I=0$ .

(5) Obtain the final regression tree model:

(17) $\hat{f}\left(x\right)=f_{M} \left(x\right)=\sum \limits_{m=1}^{M} \sum \limits_{j=1}^{J} c_{mj} I\left(x\in R_{mj} \right),$

The XGBoost classifier is based on the GB model. The main improvement is that the loss function is constructed as a Taylor expansion and a regular term is added (L1-norm or L2-norm) in the objective function:

(18) $Obj^{t} =\sum \limits_{i=1}^{n} L(y_{i}, \widehat{y_{i} }^{\left(t-1\right)} +f_{t} \left(x_{i} \right))+\Omega \left(f_{t} \right)+a,$

where $ \sum \limits_{i=1}^{n} L(y_{i}, \widehat{y_{i} }^{\left(t-1\right)} +f_{t} \left(x_{i} \right)) $ denotes the loss function, $\mathrm{\Omega }\left(f_t\right)$ is the regularization term, and $a$ is a constant.

The complexity of the XGBoost model is controllable and overfitting can effectively be avoided by addition of the regular term. Additionally, this model also has the advantage of parallel processing and fast and high flexibility. Therefore, it is highly suited to MRI data.

As described above, the process of classifying SCZ MRI images is given in Table 1.

For binary classification, a sample is divided into positive or negative classes. Specifically, four cases will occur and are given in Table 2:

True positive (TP) denotes the classification result is positive in the case of clinical normality. False negative (FN) denotes the classification result is negative in the case of clinical normality. False positive (FP) denotes the classification result is positive in the case of clinical abnormality, and true negative (TN) denotes that the classification result is negative in the case of clinical abnormality.

The true positive rate (TPR) indicates the proportion of positive samples predicted by the classifier to all positive samples, and can be described by:

(19) $\mathit{TPR}=\frac{\mathit{TP}}{(\mathit{TP}+\mathit{FN})}\times 100\mathrm{\%.}$

Similarly, the false positive rate (FPR) is calculated from:

(20) $\mathit{FPR}=\frac{\mathit{FP}}{(\mathit{FP}+\mathit{TN})}\times 100\mathrm{\%.}$

The FPR represents the proportion of negative samples that the classifier mistakes for positive classes. The ROC is obtained by plotting FPR (X-axis) against TPR (Y-axis). Accuracy (ACC) and area under ROC (AUC) are usually used to measure the performance of the classifiers [10], where ACC is defined as:

(21) $\mathit{ACC}=\frac{(\mathit{TP}+\mathit{TN})}{(\mathit{TP}+\mathit{TN}+\mathit{FP}+\mathit{FN})}\times 100\mathrm{\%.}$

MRI data was collected by the Biomedical Image Computing and Analytics Center, Department of Radiology, and Department of Psychiatry, University of Pennsylvania (USA). 132 NC and 137 SCZ subjects were recruited. All subjects met DSM-IV criteria and were diagnosed as schizophrenic by psychiatrists. All sMRI image scans were acquired using a GE 3-T Signa scanner (GE Medical Systems, Milwaukee WI, USA) with the following protocol: slice thickness = 1 mm, TE = 3.2 ms, TR = 8.2 ms, flip angle = 12 ${\mathrm{}}^{\circ }$ , acquisition matrix = $256\times 256$ , FOV = 25.6 cm. All subjects remained quiet, without moving, eyes closed, no sleeping, and minimal cognitive activity during fMRI scanning. Subjects had no history of other neurological diseases or serious drug diseases. Written informed consent was obtained from all subjects before fMRI scanning. To obtain a balanced number of controls and subjects, and also simplify the study, 100 NC and 100 SCZ were randomly chosen.

The acquired MRI images were preprocessed using the statistical parametric mapping software package SPM (Wellcome Trust Centre for Neuroimaging, Institute of Neurology, London, UK, http://www.fil.ion.ucl.ac.uk/spm) which included the following steps: (1) skull stripping, (2) bias correction, (3) tissue segmentation into gray matter, white matter, cerebrospinal fluid, and lateral ventricles, (4) spatial registration to a Montreal Neurological Institute (MNI) template, (5) generation of the regional analysis of volumes examined in normalized space (RAVENS) [21, 22] maps of gray matter, white matter, and cerebrospinal fluid by the publicly available DRAMMS deformable registration package [23], and (6) the RAVENS maps were then smoothed using a six millimeter full width at half maximum Gaussian filter.

The statistical information of 200 age and gender matched subjects is analyzed. No statistically significant is assumed for $p<$ 0.05. Please see analytical results in Table 3.

Two different training and testing sets were analyzed in this study. 80% of the 200 samples are used as training sets, and the remaining samples are used as a test set. The extracted texture features are put into each classifier in which XGBoost, GB, SVM algorithms are related to more parameters. Parameter choice directly affected classification accuracy. Therefore, parameter tuning is particularly important. Grid search and cross validation methods were used by the classifiers for parameter tuning in which parameter K = 5 in KNN. A radial basis function (RBF) is selected as the kernel function in SVM. In XGBoost, the overall parameters were divided into three categories by the authors: general parameters that guide the overall functioning, booster parameters that guide the individual booster (tree/regression) at each step, and learning task parameters that guide the optimization performed. Parameters of GB were divided into two categories: tree-specific parameters that affect each individual tree in the model and boosting parameters that affect the boosting operation in the model. The main parameters of these classifiers are given in Table 4.

The best classification results obtained by each classifier are given in Table 5. ROC curves of the GB and XGBoost classifier are given in Fig. 3.

It can be seen from Table 4 that on the same data set the XGBoost classifier gives the highest accuracy (up to 72%) and the LR classifier returns the lowest accuracy (52%). GB and KNN classifiers are inferior to XGBoost and classification accuracy is 64% and 60%, respectively. The RBF-SVM classifier is a poor classifier for this SCZ data set. It indicates that XGBoost, which has a better generalization ability and avoids over-fitting effectively, is better than the other classifiers. XGBoost also has an advantage when dealing with irregular data.

It can be seen from Fig. 3 that, compared with the GB algorithm, a higher AUC value is obtained by using the XGBoost classifier, which shows the XGBoost classifier exhibits faster and more efficient processing. Additionally, if better experimental results are required, a more detailed analysis of the features and parameter tuning of the models must be undertaken.