The study introduces a novel method for automatic segmentation of vertebral column tissue from MRI images. The paper describes a method that combines multiple stages of Machine Learning techniques to recognize and separate different tissues of the human spine. For the needs of this paper, 50 MRI examinations presenting lumbosacral spine of patients with low back pain were selected. After the initial filtration, automatic vertebrae recognition using Cascade Classifier takes place. Afterwards the main segmentation process using the patch based Active Appearance Model is performed. Obtained results are interpolated using centripetal Catmull–Rom splines. The method was tested on previously unseen vertebrae images segmented manually by 5 physicians. A test validating algorithm convergence per iteration was performed and the Intraclass Correlation Coefficient was calculated. Additionally, the 10-fold cross-validation analysis has been done. Presented method proved to be comparable to the physicians (). Moreover results confirmed a proper algorithm convergence. Automatically segmented area correlated well with manual segmentation for single measurements () and for average measurements () with . The 10-fold cross-validation analysis () confirmed a good model generalization resulting in practical performance.

1. Introduction

Pathology of the intervertebral disk is one of the common causes of pain in the lumbar spine. In 40% of cases, pain of the lumbosacral spine is diagnosed as a discogenic [1]. What is more, 80% of the general population will have or already have had pain of the lumbosacral spine [24], in 5–10% of them a chronic pain develops [1, 5].

In contemporary diagnostics Magnetic Resonance Imaging (MRI) is the modality of choice for intervertebral disc visualization. Magnetic Resonance Imaging, for almost all spinal disorders, provides robust images of the spine [6] with high quality soft-tissue visualization, much more detailed than results obtained with other modalities [7]. The additional advantage of the MRI is the lack of radiation.

Automatic tissue segmentation from Magnetic Resonance Imaging data is a challenging task, because the quality of the data affects the process; what is more, the differences between medical facilities, used protocols, and imaging machines force the necessity of universality.

Till now multiple approaches have been presented. Dong and Zheng in [8] divided the common solutions into methods that rely on graphical model [9], probabilistic model [10], watershed algorithm [11], atlas registration [12], graph cuts [13], Statistical Shape Model [14], anisotropic oriented flux [15], and random forest regression and classification [16].

The methods mentioned above are based on discrete classification returning a limited and inaccurate information about the tissue. The paper describes a method that combines multiple stages of Machine Learning (ML) [17] techniques to recognize and separate different tissues of the spine.

The objective of this study is to introduce a novel method for automatic segmentation of vertebral column tissue from MRI images.

2. Materials and Methods

2.1. The Data

For the needs of this paper 50 MRI examinations presenting lumbosacral (LS) spine of patients with low back pain were selected. The examinations were made with Siemens MAGNETOM Spectra 3T MR device. For vertebral body recognition T1 TSE (Turbo Spin Echo) Sagittal sequences, with Echo Time 9.3 ms and Repetition Time ranging from 550 ms to 700 ms, were chosen. The image sets consisted of 17 to 31 images with 4 mm slice thickness, 4.8 mm slice distance, and  px resolution.

2.2. General Procedure for Segmentation

Presented solution is based on well-known Machine Learning (ML) [18] techniques combining Cascade of Boosted Classifiers [1921] with patch based Active Appearance Model (AAM) [22, 23] algorithm and Principal Component Analysis (PCA) [24] (Figure 1).

At the beginning DICOM images are read. After that initial filtration is made to increase the quality of the data. Afterwards automatic vertebrae recognition using Cascade Classifier [21, 25] takes place. After the initial recognition the main tissue segmentation process is made using the patch based Active Appearance Model [22, 23]. Combined information about location, shape, and appearance provides a high quality model used for search and extraction of desired tissue. The results are afterwards interpolated using centripetal Catmull–Rom splines [2628].

2.3. Initial Filtration

Due to low quality of the data (low resolution, intensity inhomogeneity, and high noise), initial filtration is needed. At the beginning the images are being resized to increase resolution. For the needs of presented method a high-resolution cubic spline has been chosen [29] (Figure 2).

After that the developed intensity inhomogeneity (IIH) correction method is performed. The method is based on recalculating local intensities in such a way to fit the global exponential function defined from the boundary fat-skin tissue intensity contrast. After calculations a nonlinear selective Gaussian Blur [30] using the same global exponential function for parameterization is performed to remove the noise amplified through the correction process. As a result of this method, an intensity inhomogeneity correction is achieved (Figure 3).

2.4. Preliminary Vertebrae Recognition

At the beginning, to achieve accurate segmentation results and reduce number of MRI examinations needed for training, vertebrae recognition is made. The goal of this action is to extract each vertebra from the whole image containing spine MRI examination. To achieve this the Machine Learning [18] training of Cascade of Boosted Classifiers [1921] based on extended set Haar-like features [31] was made.

The vertebrae recognition consists of two major stages: training the classifier and vertebrae detection. Both were done using OpenCV library [25]. For the training two types of information are needed: positive examples presenting desired object that one is looking for and negative examples presenting background. To prepare the data, special software allowing fast cutting, artificial data generation, and automatic background reconstruction was developed. For the training process 50 MRI examinations were used. From those examinations over 1000 vertebrae images were extracted manually and used for automatic creation of 10,000 artificial positive examples with Thin Plate Splines (TPS) transformations [32, 33]. Afterwards negative examples were reconstructed from the same examinations by covering previously cut out vertebrae using Image Inpainting method [34].

Both positive and negative examples are afterwards used for classifiers training based on the AdaBoost algorithm [35]. Multiple weak classifiers are then combined in a cascade resembling Decision Tree [36] creating a strong classifier [37]. To achieve best performance, after the recognition, additional size constraints were introduced, removing the false positive hits. Obtained model allows proper vertebrae recognition (Figure 4).

2.5. Tissue Segmentation

After the vertebrae recognition the main tissue segmentation is made. The solution is based on Active Appearance Model (AAM) [22, 23, 38] algorithm and combines a Statistical Shape Model based on Principal Component Analysis [39], with a gray-level Appearance Model. The method focuses on recognizing the predefined characteristic features from previously extracted vertebrae images by combining the information about each pretrained characteristic feature appearance with the information about features’ mean position, their arrangement, and possible deviation. Similarly to preliminary vertebrae recognition, the tissue segmentation procedure consists of two stages: training and detection; however, contrary to previously trained classifiers, the built model is used for recognition of small patches instead of a whole vertebra. Each image used for the training originates from the prepared vertebrae database and was previously manually labeled by the group of five experts (physicians trained in MRI images assessment) with 16 characteristic points corresponding to vertebra features. Introduced information is used for building Point Distribution Model and creating training examples for the Appearance Model. The Point Distribution Model is used in a PCA [39] analysis to obtain the Shape Model containing information about the mean shape, eigenvectors, and eigenvalues. The positive and negative training examples are used for training to obtain the Appearance Model. Trained AAM model is afterwards used for spine tissue detection and classification. The detection procedure starts with an initial guess based on a perturbed ground truth shape. For this study a patch based AAM approach [22, 23, 38] has been chosen, representing the appearance of features as a rectangular patches distincted around each landmark. Finally optimization of the cost function is solved by Lucas–Kanade Optimization [40, 41] method with Wiberg Inverse Compositional algorithm [4244] (Figure 5).

2.6. Shape Interpolation

Automatically extracted 16 feature points for each vertebra image visible in the MRI examination are afterwards used for spine tissue segmentation. The information between the points is interpolated with centripetal Catmull–Rom Splines [2628] (Figure 6), ensuring C1 continuity, proper tightness with no self-intersubsubsections, and knot parameterization, leaving an area for further curve optimization.

3. Results

The method was tested on a set of 50 previously unseen vertebrae images. The spine tissue was manually segmented by 5 physicians and compared with Machine Learning results. For the numerical evaluation three measures were used [4547]: True Positive Fraction (TPF) (1), False Negative Fraction (FNF) (2), and False Fraction (FF) (3):where true positive area , is a manually segmented (by an expert) tissue area and is an automatically segmented (by a computer) area.where false negative area .where false positive area .

To achieve better performance five different optimization algorithms, available in Menpo Framework [22, 38], were tested.

Only certain Inverse Compositional algorithms [38, 4143, 4852] were chosen to be tested: Wiberg Inverse Compositional (WIC) [38, 5355] algorithm, Simultaneous Inverse Compositional (SIC) [49, 52] algorithm, Project-Out Inverse Compositional (POIC) [41] algorithm, Alternating Inverse Compositional (AIC) [42, 48] algorithm, and Modified Alternating Inverse Compositional (MAIC) [42, 48] algorithm. Three algorithms (WIC, AIC, and MAIC) achieved almost identical results (Table 1) (Figures 7, 8, and 9). Because of the best stability and lowest standard deviation (Table 1) Wiberg Inverse Compositional algorithm was chosen for further calculations.

What is more, to achieve reliable results, a mean value obtained from 100 procedure passes with 25 algorithm iterations each was computed and compared to results obtained manually by five experts (Table 2). The False Fraction is a general segmentation evaluation measure and is defined by the difference between manually segmented area and automatically segmented area, divided by the total area resulting from the manual segmentation. In this case the AAM algorithm () proved to be almost identical to expert 5, has almost the lowest standard deviation (), and is almost unnoticeably worse than other experts. False Negative Fraction provides information about percentage of nonselected pixels classified by the investigators as a spine tissue and is an amount of manually segmented area not indicated by the automatic segmentation, divided by the total area resulting from the manual segmentation. The AAM has the highest False Negative Fraction (, ) of all investigators. The True Positive Fraction provides information about percentage of properly segmented pixels and is an amount of automatically segmented area consistent with manual segmentation, divided by the total area resulting from the manual segmentation. The AAM method has the True Positive Fraction value of and standard deviation of .

What is more, a test validating algorithm convergence by comparing automatic segmentation per iteration results with manual segmentation results was performed. A single algorithm pass with 25 iterations for a set of 50 previously unseen vertebrae images was executed and a change of TPF, FNF, and FF values for each iteration was calculated (Table 3). The True Positive Fraction increases every iteration (Figure 10), while the False Negative Fraction decreases (Figure 11) simultaneously leading the False Fraction to increase per iteration (Figure 12), confirming proper functioning of presented algorithm—its convergence to manually segmented data.

The Intraclass Correlation Coefficient (ICC) was calculated to evaluate the consistency of the vertebral bodies area determined by the experts and the computer (Tables 4 and 5). The high ICC results for single measurements () and for average measurements () with the confirmed that the automatically and manually obtained segmentation results are comparable.

Additionally, the 10-fold cross-validation [5658] analysis has been done. The database of 1000 training images was divided into equal parts and tested iteratively 10 times by training the model from 90% of images and performing a test on the remaining 10% with ground truth annotations. The ground truth landmarks have been used for shape interpolation using Catmull–Rom splines and compared with automatic segmentation results using TPF, FNF, and FF measures (Table 6). The False Fraction mean value of with standard deviation confirmed a good model generalization to an independent dataset and resulting practical performance.

4. Discussion

Low resolution of presented data, high noise, and nonhomogeneous information about the tissues enforced the increasing of the quality of input data by initial filtration. From multiple interpolation methods [29, 59] widely used for image resampling, a high-resolution cubic spline has been chosen to increase the resolution of the input data, because of its good high-frequency response and high-frequency enhancement. What is more, a novel method of intensity inhomogeneity correction useful for sagittal MRI spine images has been presented. In last years, multiple methods used for intensity inhomogeneity correction emerged [60, 61]; however, they were mostly used for and tested with brain MRI scans. Because of a different application, segmentation of bone tissue instead of brain tissue, an additional method of initial filtration was developed.

For MRI images, a robust method for spine segmentation was prepared. The procedure combined well-known and widely tested Machine Learning methods [18]: Cascade of Boosted Classifiers [1921] based on extended set Haar-like features [31] for preliminary vertebrae detection, with patch based Active Appearance Model [22, 23, 38] and Principal Component Analysis [39] for precise tissue segmentation. Usage of feature localization method and interpolation of the resulting information with centripetal Catmull–Rom splines [26, 27] omitted the problem of low quality. Due to the nature of Catmull–Rom splines [28] further optimization can be done to achieve better interpolation results. The paper [62] presents multiple recent methods for intervertebral disc segmentation, which can be treated as a similar task, including Machine Learning and deep learning based approaches. The segmentation results presented in the paper [62] were measured with dice overlap coefficients and varied from 81.6% to 92% for different methods. Comparing those results with obtained segmentation and generalization results of and , one can conclude that presented AAM approach provides a good segmentation performance and moreover can be applied for intervertebral discs localization and segmentation.

In the future, automatically defined landmark localizations could be used for automatic creation of a discrete (Figure 13) and continuous (Figure 14) 3D spine models, which can be easily used in Finite Element Analysis [6365], contrary to standard voxel representation.

Obtained three-dimensional model (Figures 13 and 14) contains information about the size and shape of the intervertebral disk and the adjacent vertebral bodies. Based on it, one can determine morphology of the intervertebral disk including direction, dimensions, and volume of the herniation of intervertebral disc, giving the clinicians a tool for better understanding of the pathology.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.


This work was supported by the National Centre for Research and Development under Grant no. PBS3/B9/34/2015.