Low grade glioma segmentation using an automatic computational technique in magnetic resonance imaging

Segmentación de glioma de bajo grado usando una técnica computacional automática en imágenes de resonancia magnética

The segmentation of anatomical structures of the human brain, present in images acquired by any imaging modality, constitutes the starting point for the diagnosis of a large number of diseases or pathologies affecting the brain. Among these are brain tumors, which originate from various cell lines and are classified according to several criteria¹. One of them is where in the body they are generated. In this sense, they can be classified into two groups: a) Primary tumors. Space-occupying lesions composed of cells (SOLC) that start in the brain and tend to remain there. b) Secondary tumors. SOLC that originate in other sites of the human body and spread and/or infiltrate the brain, as metastasis.The most frequent metastases come from cancers in skin, lungs and breast^2,3.

The World Health Organization (WHO)⁴ uses the degree of malignancy of the tumor as a criterion to classify primary tumors into four grades. According to this classification, tumors labeled with grades I and II are generally benign; while those classified in grades III and IV are considered malignant. Normally, patients with primary grade I brain tumors have a longer survival than those with grade IV tumors^1,2,3.

Gliomas are primary tumors arising from glial cells and they can also be classified by grade according to the WHO grading system⁴, grading gliomas by histopathological criteria; the grade of gliomas correlates well with prognosis⁵.

Furthermore, gliomas can be classified into high-grade gliomas (HGG) and low-grade gliomas (LGG). Thus, anaplastic astrocytomas and glioblastoma multiform tumors are included in HGG; while grade II astrocytoma (Astrocite) and oligoastrocytomas, are included in LGG⁶. The present study focuses on tumors reported in the literature as LGG, that represent approximately 27% of all primary brain tumors. The median age at time of diagnosis ranges between 43 and 48, depending on histologic subtype⁷.

According to WHO⁴, the LGG grade II astrocytoma can be located as illustrated in figure 1.

Cerebral digital neuroimages are usually accompanied by various imperfections such as noise^8,9,10 and artifacts⁸ which affect the quality of information associated with the anatomical structures that make up these images. These imperfections become real challenges, when computational strategies are implemented to generate the morphology (normal or abnormal) of the mentioned structures⁸. By way of example, figure 2, generated by magnetic resonance images (MRI), illustrates the presence of Rician noise, the stair artifact, and the low contrast between brain structures and the meningioma (structure with red cross in its center).

Reviewing the state of the art regarding tumor segmentation, the works described below were found. Kauss et al.¹¹, present a system based on adaptive template moderated classification for the automated segmentation of 3D MRI brain data sets of patients with low grade gliomas and meningiomas. In a validation study of 13 patients with brain tumors, the segmentation results of the automated method are compared to manual segmentations carried out by 4 independent trained human observers. It is shown that the automated method segments brain and tumor with accuracy comparable to the manual method. Cho et al.¹², develop an LGG and HGG classification system based on multi-modal image radiomics features. The results point at this method covering an area under the curve of 0.8870 and achieving 0.8981, 0.8889 and 0.9074 of accuracy, sensitivity and specificity, respectively.

In the present work, a computational technique (CT) is proposed for the segmentation of a LGG, present in a database formed by three-dimensional brain images of MSCT. This technique considers the stages of pre-processing, segmentation and post-processing Also, percent relative error (PrE)⁸is used to compare segmentations of the LGG obtained automatically and manually.

Description of the databases

The database (DB) used was provided by the Central Hospital of San Cristóbal-Táchira-Venezuela. It was acquired through the MRI modality and consists of three-dimensional images (3D), corresponding to the anatomical structures present in the head of a female patient. Numerical characteristics are presented on table 1.

Table 1. General characteristics of the database considered in the present work.

For comparison, a manual segmentation is available, generated by a neuro-oncologist, corresponding to the LGG present in the DB considered. This segmentation represents the ground truth that will serve as a reference to validate the results.

Description of the proposed computational technique, for the automatic segmentation of the LGG.

Figure 3 shows a schematic diagram synthesizing the methods that make up the proposed technique, to segment the tumor.

Pre-processing Stage

As shown in the block diagram, figure 3, this stage corresponds to the techniques: Thresholding, Erosion, Median and Gradient magnitude filter. They are described below.

- Thresholding:

The thresholding algorithms are, generally, simple structures that allow to classify, efficiently, the elements of an image considering one or several thresholds. Such thresholds can be selected by considering both the histogram of an image and the position, intensity of an arbitrary neighborhood of the element under study, often called the current element¹³. In the present work a simple threshold was considered, based on the choice of a value for a certain threshold.

This threshold allows discrimination between the anatomical structure of interest and the rest of the structures present in an image. Usually, the threshold is chosen considering the histogram of the image. One of the criteria applied to perform the aforementioned discrimination is the following: If the intensity or gray level of the current element is equal to or less than the selected threshold value, the gray level (GL) of the current element remains unchanged; while if such intensity is greater than GL of the current element, it is generally correlated with the lower level of gray present in the image being processed¹³.

- Morphological Erosion Filter (MEF):

Mathematical morphology is based on set theory, this allows the objects present in an image to be treated as sets of points. Generally, it is possible to define operations between two sets consisting of elements belonging to the aforementioned objects and a set called structuring element (SE). SEs can be visualized as neighborhoods of the element under study, which have variable morphology (shape) and size¹⁴.

In practice, mathematical morphology is implemented through various morphological filters whose basic operators are erosion and dilation^15,16. These operators are non-linear spatial filters that can be applied to binary, grayscale or color images.

In particular, the erosion (⊖) of a two-dimensional image (I), composed of gray levels, using a two-dimensional SE, is defined by Equation 1¹⁷.

where: min the minimum gray level contained in SE, (s, t) defines the size of the SE and (x, y) represents the position of the pixel under study.

According to equation 1, to apply the filter or morphological erosion operator the image considered with an SE or neighborhood, of arbitrary size, is covered, replacing the gray level of each of the elements of such image by the level of gray minimum, contained in the aforementioned neighborhood. For purposes of the present work, a cubic structuring element was considered and the size of said SE is left as a parameter to control the performance of the MEF.

- Median Filter (MF):

The median filter (MF) is also non-linear and normally, it is used to minimize the impulsive-type noise present in the gray levels of the voxels neighboring the voxel object of study18. This type of filter is characterized by the conservation of the edges of the objects present in the image and it has the advantage that the final value of the voxel, is a real value present in the image and not an average. In addition, the median filter is less sensitive to extreme values. One of the main drawbacks is that the computation time increases substantially, as the size of the neighborhood increases.

- Gradient Magnitude Filter (GMF):

The role of this filter is to detect the edges of the structures present in the images (I). The magnitude of the gradient is often used in image analysis, mainly to identify the contours of objects and the separation of homogeneous regions²⁰. Edge detection is the identification of significant discontinuities in the level of gray or color images²¹. This technique calculates the magnitude of the gradient using the first directional partial derivatives of an image. The classic 3D mathematical model, to obtain a filtered image by magnitude of the gradient (IGM), is presented by equation 2.

where: i, j, k represents the spatial directions in which the gradient is calculated.

In practice, the magnitude of the gradient of the image in each position of the voxel, object of study, is calculated using an approach based on finite differences. In the present work, central finite differences are used. Theoretically, the filter of magnitude of the gradient based on the intensity values is very susceptible to noise²¹; therefore, it is recommended to filter the image initially to improve the performance of the detector with respect to noise.

- Binary Morphological Dilation Filter (MDF):

In order to compensate the effect produced by the morphological erosion filter, the application of a morphological dilation filter (MDF) is considered, taking into account the binary image obtained by RG. The effect of morphological dilation is to enlarge the regions of the maximum intensity image. In particular, the dilation (⊕) of a two-dimensional binary image (Ib), using a bidimensional structuring element (B), is defined as the result of operating the Ib with the values of the SE under the logical operation OR¹⁴. For the purposes of this work, a cubic structuring element was considered and the size of said SE is left as a parameter to control the performance of the dilatation process.

Segmentation Stage Computer intelligence operators: Support Vector Machines (SVM).

Support vector machines (SVM) are paradigms that undergo training and detection processes, and are based on both the Vapnik-Chervonenkis learning theory and the minimization principle that considers structural risk. SVM can be considered as classification and functional approach tools^22,23.

A variant of the SVM, called the least squares vector support machine (LSSVM), can be obtained using robust statistics, Fisher discriminant analysis and replacing the system of inequations that govern the SVM, by an equivalent system of linear equations, which can be solved more efficiently²⁴. Additionally, unlike other learning-based classification systems such as artificial neural networks (NN), LSSVMs use the criterion of minimization of structural risk, which raises the generalization capacity of the mentioned machines to optimum levels, making it possible for LSSVM perform adequately in the validation process, surpassing NN in this aspect, which uses empirical risk^24,25.

In this work, the location of the seed voxel, to initialize the segmentation technique called region growth (RG)⁸, is calculated using LSSVM. There are several functions that can be considered to construct the decision surface that allow the vector support machines to identify the seed. For purposes of the present work, a Gaussian radial base function (RBF) is considered and, therefore, a formulation is obtained that depends on the hyperparameters, identified as: a) Error penalty parameter (γ). b) Parameter to control the selectivity (σ²) of the LSSVM.

In this sense, the LSSVMs call for a process of tuning of such hyperparameters. Theoretically, both parameters can assume values belonging to the range of real numbers comprised in 0 and infinity²⁴. This tuning process is necessary because it is very difficult to know, a priori, the combination of values that will generate optimal results when the LSSVM carry out the training and validation processes.

Additionally, to automatically identify the coordinates of the seed voxel, the following procedure was implemented:

i) A size reduction technique, based on bicubic interpolation, optimal reduction factor, is applied to match the one obtained in⁶. This allows to generate sub-sampled images of 64x64 pixels from filtered images of 512x512, that is, the mentioned factor is 8.

ii) A neurosurgeon selects, on the sub-sampled image, a reference point (P1) given by the centroid of the layer containing the maximum blood pool occupied by the LGG. For this point, the manual coordinates that unambiguously establish their spatial location in each considered image are identified.

iii) An LSSVM is implemented to recognize and detect point P1. For this, the processes of:

a)Training. Training circle circular neighborhoods of 10 pixels, manually traced by a neurosurgeon, containing both point P1 (markers) and regions not containing P1 (no markers) are selected as a training set. For the markers, the center of their respective neighborhoods coincides with the manual coordinates of P1, previously established.

Such neighborhoods are constructed on the axial view of a sub-sampled image of 64x64 pixels. The main reason why a single image is chosen, for each reference point, is because it is desired to generate a LSSVM with a high degree of selectivity, which detects only those pixels that have a high degree of correlation with the training pattern.

Then each neighborhood is vectorized and, considering its gray levels, the attributes mean, variance, standard deviation and median, are calculated. Thus, both markers and non-markers are described by vectors (Va) of statistical attributes, given by: Va = [mean, variance, standard deviation and median].

Additionally, the LSSVM is trained considering the vectors Va as a training pattern and intoning the values of the parameters that control its performance, γ and σ². This approach, based on attributes, allows the LSSVM to work with greater efficiency, than when using the larger vector-based approach, which only considers the gray level of the elements of an image.

The training set is constructed with a ratio of 1:10, which means that 10 non-markers are included for each marker. The tag +1 is assigned to the class made up by the markers; while the -1 tag is assigned to the class of non-markers, that is, the training work is done based on a binary LSSVM.

During training, a classifier with a decision boundary is generated to detect LSSVM entry patterns as markers or non-markers. Subsequently, due to the presence of false positives and negatives, a process is applied that allows incorporating into the training set the patterns that the LSSVM initially classifies inappropriately.

In this sense, it was considered a toolbox called LS-SVMLAB and the Matlab15 application to implement an LSSVM classifier based on a radial base Gaussian kernel with parameters σ² and γ.

b) Detection. The trained LSSVMs are used to detect P1, in images not used during training. To do this, a trained LSSVM looks for this reference point, in the axial view, from the first to the last image that makes up each of the 7 databases considered.

The validation process carried out with LSSVM allows automatic identification of the coordinates for P1 which are multiplied by a factor of 8 units, in order to be able to locate them, in the images of original size. In this way, the aforementioned coordinates are used to establish the exact location of the seed voxel required by the RG for its initialization.

Finally, as a synthesis, figure 4 illustrates the process followed to locate the seed voxel in the databases considered.

· Region growing (RG).

The region growing is an unsupervised clustering technique, which performs an iterative process that attempts to characterize each class according to the similarity between the voxels that integrate them and thus perform the segmentation⁸. The RG method allows you to group the pixels or voxels belonging to the objects that make up an image according to a predefined criterion. The RG requires a "seed" point which can be selected, manually or automatically, to extract all the pixels connected to seed⁸.

To apply the RG, to the pre-processed images, the following considerations were made: a) The initial neighborhood, which is constructed from the seed, is assigned a cubic shape whose side depends on an arbitrary scalar .. The . parameter requires a tuning process. b) As a pre-defined criterion, modeling is chosen through Equation 3.

|I(x) − µ| <mσ(3)

Where: I(x) is the intensity of the seed voxel, μ and σ the arithmetic mean and the standard deviation of the gray levels of the initial neighborhood and . a parameter that requires tuning.

· Tuning process: Obtaining optimal parameters

The adequate performance of the proposed technique requires obtaining optimal parameters for each of the algorithms that comprise it. To do this, using DB1 as a reference, modify the parameters associated with the technique to be tuned by systematically going through the values belonging to certain ranges, as described below.

o Erosion, Median and Dilation filters have the size of the observation window as a parameter. In order to reduce the number of possible combinations, an isotropic approach was considered to establish the range of values, which control the size of the aforementioned window, in turn given by the odd combinations set by the following ordered lists: (1,1,1), (3,3,3), (5,5,5), (7,7,7) and (9,9,9).

o The parameters of the LSSVM, σ² and γ, are toned assuming that the cost function is convex and developing tests based on the following steps:

Ø For the tuning of parameter γ the value of σ² is arbitrarily set and values are systematically assigned to the parameter γ. The value of σ² is initially set at 2.5. Then, γ is varied considering the range [0,100] by steps of 0.25

Ø An analogous process is applied to intone the parameter σ², that is, γ is assigned the optimal value obtained in the previous step and, a step size of 0.25 is considered to assign to σ the range of values contained in the interval [0.50].

· The optimal parameters of the LSSVM are those values of γ and σ² that correspond to the relative minimum percentage error, calculated considering the manual coordinates of the reference seed, established by the neurosurgeon and the automatic ones generated by the LSSVM.

o During the RG parameters tuning process, each one of the automatic segmentations of the LGG corresponding to the DB1 described, is compared with the manual segmentations of the LGG generated by a neurosurgeon, considering the PrE. The optimal values for the parameters of the RG (. and .), are matched to that experiment that generates the lowest value for the PrE.

Quantitative results

The optimal sizes for the Erosion, Median and Dilation filters were (3,3,3), (3,3,3) and (5,5,5), respectively. Regarding the trained LSSVM, values of 1.00 and 1.25 were obtained as optimal parameters for γ and σ², respectively. These parameters were arrived at by means of the error present when comparing the manual and automatic coordinates, considering the percentage relative error (PrE). For this minimum value of PrE, the optimal values of RG’s parameters, r and m, were 1 and 5, respectively.

The best automatic segmentation of the LGG yielded a volume of 43.5375 cm³; while the volume associated with manual segmentation of the tumor, was 44.1677 cm³. This implies that the minimum PrE was 1.43%.

Qualitative results

Figure 5, shows a 2-D view of both the original LGG and the processed versions after applying the proposed technique.

A computational technique has been presented and its tuning process shown to allow an accurate segmentation of LGG grade II astrocytoma, present in magnetic resonance images. This statement is based on the fact that the PrE obtained was very low.

The use of intelligent operators, represented by the least squares vector support machines, allowed the automatic identification of the coordinates corresponding to the seed voxel which plays a crucial role in the adequate initialization of the unsupervised grouping algorithm based on region growing.

The LGG volume is vital when deciding whether a patient is surgically treated or not. Both the size of a LGG and its location can seriously compromise the health of a patient.

Computational techniques such as the one developed in the present investigation are precise and useful to generate the exact location and precise quantification of the volume occupied by the tumor. Additionally, the LGG segmentation is useful to plan surgery and to remove as much of the tumor as possible; while avoiding parts of the brain that control vital functions.

In the future, it is planned to carry out an inter-subject validation of the developed computational technique considering a significant number of three-dimensional images, linked to patients with this type of disease.

ACKNOWLEDGEMENTS

The authors are grateful for the financial support given by the Universidad Simón Bolívar-Colombia through the 2016-16 code project.