A New Method to Segment the Multiple Sclerosis Lesions on Brain Magnetic Resonance Images.

Automatic segmentation of multiple sclerosis (MS) lesions in brain magnetic resonance imaging (MRI) has been widely investigated in the recent years with the goal of helping MS diagnosis and patient follow-up. In this research work, Gaussian mixture model (GMM) has been used to segment the MS lesions in MRIs, including T1-weighted (T1-w), T2-w, and T2-fluid attenuation inversion recovery. Usually, GMM is optimized by using expectation-maximization (EM) algorithm. The drawbacks of this optimization method are, it does not converge to optimal maximum or minimum and furthermore, there are some voxels, which do not fit the GMM model and have to be rejected. So, GMM is time-consuming and not too much efficient. To overcome these limitations, in this research study, at the first step, GMM was applied to segment only T1-w images by using 100 various starting points when the maximum number of iterations was considered to be 50. Then segmentation results were used to calculate the parameters of the other two images. Furthermore, FAST-trimmed likelihood estimator algorithm was applied to determine which voxels should be rejected. The output result of the segmentation was classified in three classes; White and Gray matters, cerebrospinal fluid, and some rejected voxels which prone to be MS. In the next phase, MS lesions were detected by using some heuristic rules. This new method was applied on the brain MRIs of 25 patients from two hospitals. The automatic segmentation outputs were scored by two specialists and the results show that our method has the capability to segment the MS lesions with dice similarity coefficient score of 0.82. The results showed a better performance for the proposed approach, in comparison to those of previous works with less time-consuming.

INTRODUCTION

One of the most common diseases of current century, causing disabilities in the young adults, especially women is multiple sclerosis (MS). Studies show that more than 2.5 million people suffer from MS in the world and this number is rapidly increasing.[1] In MS, the immune system attacks the central nervous system and destroys myelin sheath, which leads to lesion. Depending on which area in the brain is damaged, MS appears with symptoms such as lack of vision, abnormal gait, squint, and neurological disorder. However in the acute phase, the organ paralysis and even blindness is expected. Early detection of MS and estimation of its progression are critical for optimal treatment of the disease. Magnetic resonance imaging (MRI) is a powerful tool to diagnose MS and monitoring the disease activity and progression. However, because of the large amount of data which should be analyzed, manual segmenting of MS lesions leads to time-consuming and becomes a hard task. Hence, the automatic segmentation of MS lesions in brain MRIs may be a very good solution.

Recently, a number of segmentation methods have been presented. The suggested method by Souplet et al.[2] is one of the most significant ones. They used brain atlas to register the T1-weighted (T1-w) and T2-w images. They computed the belonging value of each voxel to the three different tissues; gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF). Resulted values were used to initiate Gaussian mixture model (GMM) parameters to segment all the three (T1-w, T2-w, T2-fluid attenuation inversion recovery [T2-FLAIR]) image sequences. Then Mahalanobis distance between voxel intensity and average intensity in each class was computed and in comparison to a constant threshold, lesions were detected. Although this approach was robust against noise and inhomogeneity, it fails when there are several lesions and also it does not show satisfying performance in brain atrophy cases. Complexity and time-consuming registration task are other drawbacks of this method.

Subbanna et al. presented a fully automated framework to identify the MS lesions in multi-channel MRIs.[3] Manual segmented images were used to extract intensity histograms of both tissue and lesions. Then multivariate Gaussian distributions were estimated from the histograms and by using Markov random fields (MRFs) the classification of brain tissue and lesions was done. This method relies on manual segmentation, which may carry human error and also it needs to be trained for new image sequences.

Khayati et al.[4] combined an adaptive mixtures method means, MRF and a Bayesian classifier to simultaneously classify the three main brain tissues and the MS lesions using only FLAIR images. In particular, they first proposed to segment the brain into four classes, means: WM, GM, CSF and ''others''. Afterward, inside the ''others'' class, lesions were detected as outliers, which have not correctly explained by the model. They used only FLAIR image while MS lesions may appear independently in different images.

METHODS

In this study, we introduce a new automatic MS detection strategy. It does not need atlas or training database and overcomes defects of previously mentioned methods.

The proposed procedure has four stages in order to segment the brain tissues and MS lesions. The first stage is pre-processing which includes intensity inhomogeneity correction and skull removing. In the second stage, which is called the brain tissue classification, the brain tissues are classified into three classes of WM, GM, CSF, and some voxels will be rejected as outliers. In the third stage, candidate lesions are detected from outliers with the Mahalanobis distance and in the fourth stage; MS voxels are separated from the candidate lesion.

Database

Our database images were acquired from two different hospitals in Isfahan to evaluate the performance of the proposed automatic segmentation of brain tissues and lesions in different scanning machines. Twenty-five patients were scanned by using the same protocol of T1-w, T2-w, and T2-FLAIR. All the images were acquired by two 1.5-T MR systems and in axial view with the slice thickness of 5 mm. The patients were between 11 and 45 years old and the same imaging protocol was used for them.

PREPROCESSING

To reach better and accurate segmentation, some preprocessing steps such as skull removing and intensity inhomogeneity correction are needed. So, in this step, these pre-processing were done as the followings.

Intensity Inhomogeneity Correction

Intensity inhomogeneity in raw MRIs leads to incorrect segmentation results. This inhomogeneity is because of the small varying biased field (BF); so that voxels with same coordinate have different intensities. Intensity inhomogeneity changes the mean and variance of image intensity in a particular area of the image which decreases the segmentation accuracy. Since BF alters slowly, it only contains low frequencies and, as a result, blurs the images (or destroys high frequencies). The real aim of inhomogeneity correction is to restore these high frequencies to the images.

Consider the following model image formation in MR:

v(x) = u(x) f(x) + n(x)      (1)

Where at location of x, v is the measured signal, u is the true signal emitted by the tissue, f is an unknown smoothly varying bias field, and n is white Gaussian noise, which has been assumed to be independent of u. The compensating problem of intensity nonuniformity is the task of f estimating. The combination of additive and multiplicative interference makes this task difficult.

We used low pass filter (LPF)[5] to compensate inhomogeneity. In this method, it is assumed that biased field can be estimated from image by itself. The following relation was used to extract original signal u.

u = exp(log(v)−lpf(log(v)))      (2)

This method can be summarized as the following:

Automatic determination of the lower level threshold to reduce the “noise”. To do this, in the first step, image histogram was calculated and smoothed by Hanning window. Then the first maximum at the histogram was determined in the lowest signal intensity. After that, the lower level of threshold was considered as 15% of this lowest signal intensity. The coefficient of 15% was found by experiment and assessment, and the results had agreement with the report of two radiologists

Determination of the average signal intensity in the nonnoise locations

Filling the “noise” locations with the average image intensity of the nonnoise locations

Smoothing the image data. Image volume was smoothed by three-dimensional (3D) Gaussian kernel. Experimentally, window size of this kernel was set to 3/8 of the image size. To keep computation time to a practical level, the convolutions were performed in k-space:

SIsmooth = SIraw (x, y, z)⊗G(x, y, z) = F−1[F(SIraw(x, y, z))F(G(x, y, z))]      (3)

Where F is the Fourier transform operation, F−1 is its inverse, ⊗ is the convolution operator, SIraw (x, y, z) is the image intensity of all pixels in the 3D volume, and G (x, y, z) is a 3D Gaussian function with the dimensions, which have been discussed previously

Normalize the signal intensity of the raw image by using the smoothed image and correct the intensity so that the average pixel intensity (in the volume) is retained after correction. Specifically, the image intensity at location i (SI[i]) becomes:

Where SIraw, SIsmooth, and SIaverage are raw, smoothed, and average image, respectively.

LPF method is effective and faster in comparison to the other similar intensity correction methods, which use an optimization phase to find the biased field. Results of all stages of this algorithm have been illustrated in Figure 1.

Figure 1

Result of intensity inhomogeneity correction, (a) input image, (b) automatic determination of the threshold lower level and segmentation of the input image to data (white) and background (black), (c) filling the noise locations with the average image intensity of the nonnoise locations, (d) smoothing image intensity, (e) output image

Skull Removing

Because of the similarity between skull intensity and other tissues of the brain, it may affect the segmentation phase and usually it is removed from brain image. Since the skull is more obvious and detectable in T1-w, we used this image to detect the exact location of skull and remove it from all three image sequences. The skull is brighter than other brain parts and usually has a component with no rupture. Therefore, we used the morphological operation to extract and remove the skull. The first operation is detecting the connected component [Figure 2c]. As it has been seen, skull is extracted as the first connected component. After removing this component from the brain image, still there are many pixels, which do not belong to GM, WM, and CFS, they are remains of the skull. Otsu's threshold algorithm[6] was applied to detect brighter pixels as remains of the skull [Figure 2d]. After these two phases, some remote and unconnected pixels remains which can be removed by morphological operation such as erosion, dilation, and filling. Results of each phase of this process have been shown in Figure 2.

Figure 2

Result of skull removing, (a) is the brain image of one of the 25 patient which were studied in this research, (b) the first component of the image, (c) removing the first component of the image, (d) the processed image after applying Otsu's thresholding algorithm, (e) brain mask, (f) output image

Brain Tissue Classification

Brain can be classified into three distinct classes; WM, GM, and CSF. Assuming that intensity variation in each class has Gaussian distribution, we used GMM[7] to separate these classes from each other. In GMM, K Gaussian distributions (three classes in this work) are considered and each sample (voxel) belongs to the class, which maximizes the probability distribution function. GMM parameters (covariance matrix and mean of each class) should be determined during the classification process. The first step is registering T1-w to T2-w and T2-FLAIR so that each sequence has the same number of voxels and each voxel in one sequence exactly corresponds to the same voxel in other sequences.

Each sample data is a 3D vector, which carries intensity of all three aligned voxels in each image:

Xi = (XT1(i), XT2(i), XT2−F(i))      (5)

If each image sequence has K slices and each slice is M * N, we have K * M * N samples with 3D which should be clustered to three different classes (WM, GM, and CFS). Suppose that each class distribution is independent from those of other classes and has independent variance and mean, whole probability of each sample can be formulated as follow:

Where k is the number of classes, αj is prior probability of each class (or merging index), μj is mean vector, and Σj is covariance matrix of class j. Covariance matrix of each class is diagonal (assuming that tissue intensity is identical and independent from those of other image) as follow:

GMM parameters mean covariance matrix and the mixing parameter are merged in the parameter vector θ.

These parameters can be estimated using the MLE:

In order to obtain the MLE, we can employ the EM algorithm,[8] a technique, which is used to iteratively estimate . From a given θ, the EM algorithm obtains another θ( where L(θ < L(θ(). The algorithm is generally considered to have converged when θ and θ( are sufficiently close to each other. This method is usually chosen because it is easy to implement and there is a proof of convergence, but it has some drawbacks. The first drawback is that the EM algorithm does not ensure to reach the global maximum; different initial parameters may lead to different solutions, which make the choice of it an important issue. In order to overcome this problem, EM algorithm is applied several times with different initial point and the best result is saved as global maximum. However, this is time-consuming. We used different strategies to increase the speed and accuracy of EM algorithm.[9]

Hierarchical Initialization

To make EM robust against initial point selection, we tried to choose a reasonable start point instead of random ones. GMM was applied to segment only T1-w image with 100 different starting points with the maximum iteration of 50. Segmentation results were used to calculate the parameters of other two images. Since we limited the maximum iteration and only used T1-w image, simulation time decreased considerably and also the initial points were reasonable and not random. For the random initial parameters, the mean of each class is randomly drawn using a uniform distribution between the minimum and maximum of the image and the standard deviation of each class is set to a third of the standard deviation of intensities of the whole image.

After computing the mean and variance of each class for T1-w image, we clustered each voxel of T1-w and other two images too. Now, we should estimate initial mean and variance of each tissue (GM, WM, and CFS) in T2-w and T2-FLAR images. First, for each class in T2-w and T2-FLAR images, we calculated histogram with 256 bins and smoothed it with a Gaussian window. Initial mean of tissues were set to mode value of these histograms. Furthermore, variance of each tissue (class) for T2-w and T2-FLAR images was estimated as follow:

Where t belongs to (GM, WM, and CFS) and s is one of the T2-w and T2-FLAR images and med is median operation. Now, we have all start parameters to initiate the EM algorithm and find the best parameters so that classify each 3D voxel and be sure that EM will converge to optimum maximum more likely.

Furthermore, trimmed likelihood (TL)[10] was used instead of likelihood to make our process robust against outliers and also FAST-trimmed likelihood estimator (FAST-TLE) was applied to make it faster.

Trimmed Likelihood

The basic idea consists in maximizing the TL instead of the likelihood:

Where h is the number of voxels which will be considered as outliers, n is the number of all samples, and v(i) is a function, which sorts samples as follow:

f(xv(1)|θ)≥f(Xv(2)|θ)≥...≥f(xv(n)|θ)      (12)

In other words, the likelihood is only computed with the voxels that most likely belong to the model. Hence, we need to do another optimization in order to determine which samples are remote enough and do not fit the model. To speed up this process, we applied FAST-TLE.[9] In this method, first h samples were considered as outliers, randomly and GMM parameters were estimated. In the second step, f(xv(i)) was computed and samples with less f(xv(i)) value were considered as new outliers. These two steps were done till outlier set does not change anymore. It should be noticed that h value should be large enough to guarantee that all MS regions and tissue artifact will be detected and stored in outlier set. In this research, we set h to 10% of all pixels of the sample image. The final result of this step is an outlier set and three different classes and each sample belongs to one of these four sets. MS voxels belong to outlier set and should be extracted from this set.

DETECTION OF CANDIDATE LESIONS

In practice, the rejected voxels contain some voxels that actually fit the model reasonably well. Thus, we define the distance as the minimal Mahalanobis distance of the voxel from one of the Gaussians in the model:

The voxel is considered as a candidate lesion when the distance is greater than a threshold.

Multiple Sclerosis Voxels Separation from Other Candidate Lesion

Candidate lesions detected with the Mahalanobis distance include MS lesions, vessels, registration errors, flow artifacts, noise, etc., We define five rules in order to discriminate MS lesion voxels from the other voxels:

Intensity rule – MS lesions are known to be hyperintense compared to the WM intensity on T2-w and FLAIR sequences. We used the information given by the model to define hyperintensity. If the voxels' intensity was less than a threshold in T2-w or FLAR image, it was discarded from MS

Size rule – In order to avoid false-positives (FPs), candidate lesions smaller than 6 mm in size were rejected. These small candidate lesions are usually produced by noise or flow artifacts. In clinical practice, lesions must have a radius of 3 mm on one image slice, to be considered as MS[11]

Connectivity to WM – MS lesions' pixels in one slice usually connect to WM tissue. Hence, we removed samples from MS candidate, which did not connect to WM[9]

Neighbor slices rule – MS candidate, which appear only in one slice, was rejected. Based on expertise, MS is a volume and if there is a region of lesion in one slice it should exist at least in one or two next or previous slices. Note that this differs from size rules

Shape rule – MS lesions' regions in one slice are not distributed along a line and it almost is like a circular and or oval. Therefore, candidate regions, Which were not similar to a circular or oval should be removed from final MS set. To determine the amount of being circularity, we used region's coordinate variances along x and y-axis.

Shape rule is used to reject regions, which are similar to line or an oval with small diameter of “a” and large diameter of “b” while a/b is smaller than a threshold, means a/b< threshold. So, in this research when we say resulted shape is circular, it means that a/b> threshold (in a circular, a=b). We suppose this is the simplest way to measure how much an object is similar to line or not. Now, to find “a” and “b” we used variance in x and y-direction and set “a” to the smaller value and “b” to the larger one. Note that variance is equal to mean of distance from the center and since the mean is robust against noises and some possible outliers, we can be sure to the “a” and “b” values.

Finally, all remaining voxels in the outlier set were considered as MS. The output of each step of proposed algorithm on a typical image has been shown in Figure 3.

Figure 3

Steps of proposed algorithm on a typical image, (a) input image, (b) brain tissue classification, (c) detection of candidate lesions with the Mahalanobis distance, (d) outliers after applying intensity rules, (e) outliers after applying other heuristic rules (multiple sclerosis lesions)

It should be noted that since the smoothness of MS's boundaries is not considered by GMM, the exact location of MS's boundaries cannot be extracted and always there are some little errors. In section 2.2, it was mentioned that to find outliers, it is necessary to set h value. This h value shows how many voxels have been considered as outliers in the first step. If h value is close to one, all voxels will be considered as outliers, and if it tends to zero, there will have no outliers and all voxels are belong to brain tissue. In other words, the higher values of h makes more brain tissue and noises to be outlier and probably considered as MS in final decision and smaller h values, may consider voxels locating in MS boundary as normal brain tissue. Hence, h value selecting, is a trade-off between rejecting noise and removing MS boundaries. So, the amount of errors produced by GMM is in direct relation to h value, and this error value is clarified by dice similarity coefficient (DSC) measurement.

DSC is a spatial overlap index and a reproducibility validation metric. DSC measures consistency of results and ground-truth (extracted by experts) and it is a comprehensive measurement of the coincidence of results and ground-truth. The value of a DSC ranges from 0, indicating no spatial overlap between the two sets of binary segmentation results to 1, indicating complete overlap. It is clear that higher the DSC value, more overlap between the results and ground-truth and the more accuracy.[12]

RESULTS

As it was mentioned in the previous section, h value has a great effect on the accuracy of MS lesion detection and image segmentation. The larger amount of h means, the more tissue voxels considered as MS candidate and, on the other hand, smaller value of h result means missing some MS regions. To find the optimum value of h for our database, we changed h value from 0 to 0.5 (limit of convergence) and recorded the Dice similarity factor, which was calculated by two experts. Results showed that for patient with high lesion load (≥10 cm3), optimum h is 0.25 and for low lesion load (<10 cm3) h should be in the range of 0.1–0.15. Figure 4 shows h effects on dice factor for both groups of patients with low and high lesion load.

Figure 4

Effect of h on dice similarity coefficient factor for both groups of patients with low and high lesion load

As it was referred in previous sections after the first stage means classification, there are some voxels, which considered outliers while actually they are some part of a normal tissue. To extract these voxels from outlier set, we should adapt an optimum threshold to Mahalanobis distance. After some experiments, we set it to 4. Also, we set intensity rule threshold to 130 for images with 256 gray levels.

To evaluate the performance of the proposed algorithm, output results were compared with the detected MS regions by two experts and specialist. To do so, we used dice, sensitivity, accuracy, and specificity. Results were summarized in Table 1.

Table 1

Statistics of similarity between proposed approach and manual segmentation derived from 25 MS patients

In this table, true-positive (TP) is the number of voxels, which both the proposed method and experts were detected as MS. true-negative is number of the voxels, which were not considered as MS by any of experts and our algorithm. FP is a number of the voxels that only our method was detected as MS and false-negative is the number of those which have been recognized as MS by experts, and our method detected them as normal brain tissue.

As it has been seen in Table 1, the accuracy of the results directly depends on the amount of lesion. For a patient with less amount of lesion, if only a few numbers of voxels were missed to be detected as MS, TP decreases considerably while it never happens for images with a high amount of lesions. It should be noticed that the volume of lesions is small with respect to the whole brain volume and that is why specificity values are always close to one.

CONCLUSION

In this paper, we proposed a new strategy to initiate EM algorithm without any atlas in order to reach accurate results and make the algorithm to converge rapidly. GMM was applied to segment only T1-w image with 100 different starting points where the maximum number of iteration was considered to be 50. Also, the heuristic rules were applied to reject some false alarms and find real MS. These rules help the algorithm to reject the errors introduced by GMM. We used T1-w image sequence to extract the needed information to find initial set points, which reduce the risk of being trapped in local optimums. We used FLAIR image sequence as a standard protocol to detect and locate MS lesions. Even though this sequence is sensitive to periventricular lesions, it is less sensitive to the lesions in posterior fossa, which may increase FP errors. Hence, we used T2-w as complimentary sequence to improve segmentation performance. Our method does not need any atlas or training images. In contrary to the previous methods, which use only image intensity to find MS, we tried to take the advantage of other structural and local information such as size and shape of lesions, which improved segmentation results considerably. Intensity based segmentation can be easily affected by noise and artifact. Doing many experiments and searches, we found optimum values of all thresholds.

Unfortunately, we were not able to compare our approach to other similar methods because of two main reasons. First of all, different studies used different database and second, any MS segmentation algorithm parameters were adapted for especial MR protocol. However, in the researches of Subbanna et al.[3] and Khayati et al.,[4] which were mentioned in the literature, DSC values have been reported as 0.71 and 0.75, respectively, while in our method this parameter was averagely more than 0.82. The results showed better performance of the proposed approach, compared to those of previous works and these results were confirmed by two radiologists.

Financial Support and Sponsorship

Nil.

Conflicts of Interest

There are no conflicts of interest.

BIOGRAPHIES

Alireza Karimian received his B.Sc. degree in electronics engineering from Ferdowsi University, Mashhad, Iran. He also received his M.Sc. and Ph.D. degrees in nuclear engineering in the field of medicine from Amirkabir University of Technology, Tehran, Iran. Furthermore he has passed successfully a one year fellowship research in the field of medical physics in the La Sapienza University, Rome, Italy under grant of ICTP. In 2006, He joined the Department of Biomedical Engineering at University of Isfahan, Iran, as an Assistant Professor, and was the lecturer of some courses such as medical imaging systems, Simulation and its application in medicine, radiation shielding, Dosimetry and radiation detection, Biophysics and medical physics. Since February of 2013, He was successful to be as Associate Professor in University of Isfahan. He has published more than 150 research papers in peer-reviewed journals and conferences. His research interests are: Imaging systems - Image processing – Dosimetry - Radiotherapy- Monte Carlo simulation and its applications in medicine.

E-mail: Karimian@eng.ui.ac.ir

Simin Jafari received the B.Sc. and M.Sc. Degrees in Telecommunication Engineering from the Islamic Azad University of Najafabad, Isfahan, Iran, in 2009 and 2014, respectively. Her main research area are image processing and pattern recognition.

E-mail: siminjafari11@yahoo.com