Abstract
Determination of the volumes of acute cerebral infarct in the magnetic resonance imaging harbors prognostic values. However, semiautomatic method of segmentation is time-consuming and with high interrater variability. Using diffusion weighted imaging and apparent diffusion coefficient map from patients with acute infarction in 10 days, we aimed to develop a fully automatic algorithm to measure infarct volume. It includes an unsupervised classification with fuzzy C-means clustering determination of the histographic distribution, defining self-adjusted intensity thresholds. The proposed method attained high agreement with the semiautomatic method, with similarity index 89.9 ± 6.5%, in detecting cerebral infarct lesions from 22 acute stroke patients. We demonstrated the accuracy of the proposed computer-assisted prompt segmentation method, which appeared promising to replace the laborious, time-consuming, and operator-dependent semiautomatic segmentation.
1. Introduction
Cerebrovascular disease is one of the leading causes of acute mortality and chronic disability [1]. The volume of infarct is associated with severity of acute ischemic stroke and correlates with clinical prognosis and the effect of endovascular therapy [2–4]. A rapid and reliable method of determination of volume of acute infarct will help predict the prognosis and facilitate further investigation.
The diffusion weighted imaging (DWI) is more sensitive than other magnetic resonance imaging (MRI) modalities to small water diffusion changes in the acute ischemic brain, especially within 48 hours of the ictus [5–9].
Automatic algorithms for segmentation for acute infarct in MRI have been reported [10–15]. The unsupervised method developed by Li et al. was based on a multistage procedure including image preprocessing, calculation of tensor field, measurement of diffusion anisotropy, segmentation of infarct volume based on adaptive multiscale statistical classification, and partial volume voxel reclassification [11]. Bhanu Prakash et al. used a probabilistic neural network for selecting infarct slices and an adaptive Gaussian mixture model for segmentation of the infarcts [12]. Hevia-Montiel et al. developed a method for cerebral infarct lesion segmentation from DWI by applying nonparametric density estimation [13]. Gupta et al. identified the infarct slices and the hemisphere automatically in DWI based on the difference in the percentile characteristics of intensity normalized images and parameters of infarct slice identification and infarct hemisphere identification [14]. Shen et al. detected infarct lesions based on the voxel intensity segmentation and the spatial location of tissue distribution [15].
We aimed to design a DWI-based computer-assisted method to provide clinicians a prompt and accurate determination of the volumes of acute cerebral infarct. The operation of this method is based on the histographic characteristic of the output clusters of a fuzzy C-means (FCM) clustering [16]. Additional measures were taken to ensure the accuracy of infarct detection, including discriminating infarcts from artifacts due to magnetic inhomogeneity by incorporating the histographic information in the apparent diffusion coefficient (ADC) map.
2. Materials and Methods
2.1. Subjects and Image Acquisition
Landseed Hospital has been participating in the nationwide Taiwan Stroke Registry, which prospectively registered patients with stroke onset within 10 days according to a preestablished system [17]. For this study, we recruited 22 patients (11 women and 11 men, 62–83 years of age) with acute cerebral infarction and MRI examinations during January-February 2011. The protocol of this research has been reviewed and approved by the Institutional Review Board (IRB) of Landseed Hospital.
All MRIs were acquired with a Signa HDxt 1.5T Optima edition (GE Healthcare, Waukesha, WI) and consisted of a DWI scan (TR/TE/Flip angle = 6000 ms/82.8 ms/90°, FOV = 230 × 230 mm2, matrix = 128 × 128, in-plane resolution = 1.79 × 1.79 mm, 24 axial slices, 5 mm slice thickness with 1 mm gap) and an ADC map with s/mm2.
2.2. Automatic Infarct Detection Procedure
The proposed method was developed on a personal computer with Intel Core i5 CPU, 2.67 GHz processor speed, and 4 GB RAM. The infarct detection procedure was carried out mainly with a MATLAB program (The MathWorks, Inc., Natick, MA). We utilized the histographic characteristic of the DWI for infarct detection. Figure 1 illustrates one example with a large infarct volume and another with a small infarct volume. The procedure comprised the following steps.
(a)
(b)
(c)
(d)
(e)
(f)
Step 1 (coregistration and intensity normalization). The ADC map was registered to the corresponding DWI by a rigid registration (translation and rotation) and a trilinear interpolation based on the normalized mutual information method to correct for differences due to head movements [18]. The DWI and registered ADC map were normalized so that their intensities were both distributed in a standardized range . The program we used to run this step was the Statistical Parametric Mapping 8 (SPM8, Wellcome Department of Cognitive Neurology, London, UK).
Step 2 (extracting the brain mask from the whole-brain DWI). The brain mask was extracted from the whole-brain DWI based on the estimation of the inner and outer skull surfaces by using BET (Brain Extraction Tool), a software package developed at FMRIB Centre, University of Oxford, Oxford, United Kingdom [19]. Note that the fractional intensity threshold was set at 0.3, smaller than the default value 0.5, to give a larger brain outline estimate, which would completely enclose the brain. On one hand, the brain mask extraction would not eliminate any portion of the cerebral infarcts. On the other hand, no part of the brain skull portion enclosed by the brain outline would be mistaken as an infarct region in the subsequence steps, because there are obvious intensity differences between the cerebral infarcts and brain skull.
Step 3 (preclustering elimination). The histogram of DWI within the brain mask was smoothed by a third-order moving-average filter. The peak of the smoothed histogram was identified. The normalized intensity, denoted by , corresponding to this histographic peak was used as a threshold. The voxels with normalized intensities lower than or equal to would be eliminated from further processing.
Step 4 (FCM clustering). The remaining voxels after the previous step were divided into 50 clusters by an unsupervised classification with the conventional FCM clustering algorithm.
Step 5 (skimming the clusters for candidate voxels). The clusters with mean normalized intensity larger than the normalized intensity of the histographic peak plus 0.2, that is, , were selected. The voxels belonging to these selected clusters would be treated as candidate voxels of infarct in the next step.
Step 6 (eliminating labels with insufficient intensity). Each cluster of candidate voxels of infarct was further divided into one or several labels, at most the same number as the voxel number in that cluster. A label comprised connected voxels. The labels with average normalized intensity lower than or equal to the threshold () were eliminated from further processing.
Step 7 (eliminating labels with weak edge). The labels with weak edge were eliminated from further processing. The edge map was extracted from normalized DWI by using Canny edge detector [20]. The low and high threshold values were defined as . The parameter of the standard deviation of the Gaussian filter was determined to be 1.
Step 8 (eliminating candidate labels due to magnetic inhomogeneity). Let denote the intensity corresponding to the histographic peak of the ADC map. Further, let denote the average intensity of the lower-intensity half of all the voxels in a specific label on the ADC map. For this specific label, if the ratio , this label was considered an artifact caused by magnetic inhomogeneity. All artifacts thus defined were detected and eliminated in this step. Finally, all the voxels in the remaining labels in the remaining clusters were taken to be infarct.
2.3. Performance Evaluation
A voxelwise comparison between the proposed automatic segmentation and semiautomatic segmentation by the experienced neurologist (Chen) [21] gives the four parameters of each patient: true positive (TP), true negative (TN), false positive (FP), and false negative (FN). The sensitivity (Sen.), specificity (Spe.), positive prediction value (PPV), and negative prediction value (NPV) are calculated from the four parameters [22].
The similarity index (SI) is used to indicate the degree of agreement between the infarcts detected by our method and those detected semiautomatically by the neurologist. Its formula is [23]. In addition, Cohen’s kappa coefficient is also calculated. It eliminates the agreement due to random chance and is considered a conservative measure of interrater agreement. Infarct volume was calculated as the summation of the detected infarct area of axial DWIs times the slice thickness [24]. The agreement evaluation of volume measurements of the proposed algorithm was carried out by calculating the intraclass correlation coefficient (ICC) [25].
2.4. Preliminary Experiment
A preliminary experiment was conducted to find the most suitable cluster numbers for the FCM clustering in Step 4. Cluster numbers ranging from 6 to 100 were tested. For each of the tested cluster numbers, a semiprocedure of the proposed method was conducted on each of the 22 recruited patients to obtain an average SI. The semiprocedure consisted of Steps E1 to E7, whereof Steps E1 to E5 were the same as Steps 1 to 5. Steps E6 to E7 were as follows.
Step E6 (selecting infarct labels). Each cluster of candidate voxels of infarct was further divided into one or several labels, at most the same number as the voxel number in that cluster. A label comprised connected voxels. All the labels containing at least one voxel belonging to a semiautomatically demarcated infarct region were selected as infarct labels.
Step E7 (SI calculation). The SI was calculated with the voxel-by-voxel comparison between the infarct labels selected in Step 6 and the semiautomatically demarcated infarct regions.
3. Results
The exemplary images in Figure 2 illustrate the procedure of the proposed method. Figure 2(a) shows an example of the DWI slice used as the input to the proposed method in Step 1. Figure 2(b) was a slice of the whole-brain mask extracted from the whole-brain DWI in Step 2. Figure 2(c) shows a slice of the output of the preclustering elimination in Step 3. These were the DWI voxels with normalized intensities higher than within the brain mask. In Figure 2(d), different colors were used to paint the voxels of different clusters in an exemplary DWI slice after the FCM clustering in Step 4. Figure 2(e) shows the corresponding Canny edge detection map. Figure 2(f) shows an example of the final detected infarcts after further processing through Steps 5, 6, 7, and 8. Figure 2(g) is a combination of the detected infarcts and the raw DWI. Figure 2(h) shows the ADC map that was used in Step 8 to eliminate artifact-induced spurious infarcts. Figure 2(i) shows the result of the semiautomatic infarct segmentation by the neurologist on the same input DWI.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
Figure 3 shows how the histographic information of the voxel intensity was utilized in the proposed method to facilitate the identification of infarct. The blue curve in Figure 3(a) represents the smoothed histogram of the voxel intensity of a normalized raw DWI. The intensity corresponding to the peak of the histogram is referred to as . In Step 3, all voxels with normalized intensity lower than or equal to were eliminated because they very unlikely belonged to infarct areas. This preclustering elimination greatly reduced the computation load in the latter steps. Each of the 50 dots in Figure 3(b) represents the average normalized intensity of all the voxels in an individual cluster among the 50 output clusters of the FCM clustering in Step 4. In Step 5, only the clusters with mean normalized intensity higher than or equal to were skimmed (selected) as candidate infarct clusters. In the example shown in Figure 3(b), only one cluster became a candidate cluster. Each dot in Figure 3(c) represents the average normalized intensity of an individual label in the candidate cluster(s). The labels with average normalized intensities lower than were eliminated in Step 6. Figure 3(d) shows the histogram of the normalized voxel intensity of the ADC map. It was used in Step 8 to distinguish artifacts due to magnetic inhomogeneity.
(a)
(b)
(c)
(d)
Verified with the reference by the experienced neurologist (Chen), our algorithm had high sensitivity, specificity, and SI. Table 1 tabulates the performance indices of our algorithm on each of the 22 subjects. Note that we conducted the proposed method 10 times on each patient’s MRI to get a (mean ± standard deviation), shown in the supplementary table (see Supplementary Material available online at http://dx.doi.org/10.1155/2014/963032), of each performance index. Table 1 only shows the mean values of the performance indices. With total infarct lesion volume ranging from 0.155 to 482.939 mL, the sensitivity was , the specificity , and the SI . The standard deviations shown in the supplementary table reveal that the variation due to FCM clustering was low and acceptable.
Figure 4 shows the relationship in the infarct volumes obtained by the semiautomatic method and the proposed method. The infarct volumes determined by the proposed method correlated well with those determined by the semiautomatic method with an ICC of 0.991. Notice that patient number 22 was an outlier, with a much higher infarct volume than those of the other patients. This patient is not included in the plot of Figure 4. The ICC with patient number 22 included was 0.993.
The result of the preliminary experiment with various output cluster numbers of the FCM clustering algorithm is shown in Figure 5. It demonstrates that the number of the output clusters has substantial influence on the average SI value. Figure 6 illustrates the result of the FCM clustering by using different colors to represent the 50 different clusters in a whole-brain DWI.
To illustrate the effect of the proposed algorithm, the semiautomatic and the proposed automatic segmentation results on patient #9 are displayed in Figure 7.
(a)
(b)
(c)
(d)
It took our personal computer less than 90 seconds to execute from Step 1 through Step 8. However, because the SPM program in Step 1 and the BET program in Step 2 required user intervention, the whole procedure actually took nearly 5 minutes. In the future, if the functions currently executed by the SPM and BET programs are integrated into the main MATLAB program, the whole procedure will finish within 90 seconds.
4. Discussion
We reported a high SI of our algorithm. The method proposed by Bhanu Prakash et al. attained 60% in SI [12]. The SI of their proposed method two years later [10] was improved to 67%. Another method [14] proposed by the same team attained 67.6% in SI. The method [15] by Shen et al. could attain 87.9% in the average SI, obtained with simulated lesions. The method [11] proposed by Li et al. attained SI above 92%. In comparison, the performance of our method in terms of SI is higher than those of Prakash and Gupta’s team, similar to the method of Shen et al., and only next to that of the method of Li et al. It is worth noting that the SI values reported by different research teams were based on different gold standards and calculated from different experiment setups. Hence, there is no fair comparison among the performances of these methods.
We attribute the high SI of our method to the following key points.
First, we used a self-adaptive threshold, that is, , in the preclustering elimination step (Step 3). The purpose of the preclustering elimination was to increase the efficiency of fuzzy clustering and reduce the computation time. The histographic peak corresponded to the most abundant intensity of the DWI image. The voxels with intensity lower than were very unlikely to belong to an infarct, so it is safe to eliminate them from further processing. The threshold was self-adaptable in that the of every patient served the purpose well.
Second, we used an optimal value for the output cluster number of the FCM clustering algorithm in Step 4. It was chosen to be 50 since this was a value that would lead to high SI values, as was demonstrated in Figure 5.
Third, in Step 5, we selected the mean demarcation threshold, , used by the neurologist as the threshold for skimming the clusters for candidate voxels. As illustrated in Figures 8(c), 8(d), and 8(e), different extents of skimming led to different SI values. We found that the difference between and the minimum normalized intensity of the infarct regions demarcated semiautomatically by the neurologist averaged around 0.2. In the proposed algorithm, the normalized intensity value that was 0.2 higher than was selected for the value from which the skimming started. Note that the demarcation threshold had been determined based on the statistics of the experienced neurologist’s semiautomatic demarcation results. The level of this threshold can be changed to accommodate for different neurologists, scanners, and acquisition parameters. A suitable new threshold can be obtained by statistical analysis on a training set, making the proposed algorithm adoptable to all situations.
(a)
(b)
(c)
(d)
(e)
Fourth, we eliminated false-positive labels in the candidate clusters. The results before adopting Step 6 had shown that there would be some false-positive infarct regions in the result if all the voxels of the skimmed clusters were taken as infarct. In Step 6, by dividing each skimmed cluster into labels and eliminating labels with low intensity, the false positives could mostly become true negatives. As illustrated in Figure 9(a), the labels with average normalized intensity higher than or equal to were retained in the candidates and the others were eliminated. Figure 9(b) exemplifies a true-positive infarct label, which looked white and might have a faint suburb. Such a label was classified as an infarct region, painted green in Figure 9(c). Figure 9(d) shows some labels with average normalized intensity lower than ; these labels looked faint all over. Such labels were classified as noninfarct regions, painted red in Figure 9(e).
(a)
(b)
(c)
(d)
(e)
Fifth, edge detection was used to further eliminate labels with weaker edges. The results of Step 6 still contained a few false-positive labels. That is to say, some noninfarct labels had sufficient normalized intensity to pass the decisions in Steps 5 and 6. However, these false-positive labels had weaker edges than real infarct labels had. Step 7 used this edge for eliminating the false-positive labels. Figure 10 demonstrates the function of this step with real examples.
(a)
(b)
(c)
(d)
(e)
(f)
Sixth, we detected and eliminated the false-positive labels due to magnetic-inhomogeneity artifacts. The result of Step 7 still contained false-positive labels, which were actually artifacts caused by the magnetic susceptibility differences between adjacent air and cerebrospinal fluid structures and the surrounding soft tissues with echo-planar imaging techniques [26]. In the DWI, the magnetic inhomogeneity created artifacts with intensities commensurate with those of infarcts. It had been a difficult task to eliminate the artifacts [12]. However, in the ADC map, the artifact intensity was higher than the infarct intensity. This property was used in Step 8 to detect and eliminate the artifacts caused by magnetic inhomogeneity. Figure 11 illustrates this phenomenon.
(a)
(b)
(c)
The proposed method had the second lowest SI (79.722%), despite of 100.000% sensitivity, for patient number 13 among the 22 patients. That was actually due to inconsistent selection of threshold by the neurologist. Notice from Table 1 that the demarcation boundary selected by the neurologist for patient number 13 was , which was much higher than the average value. In other words, the neurologist included fewer voxels into infarct than usual. The demarcation boundary selected by the neurologist for patient number 14 was , much lower than the average value. That is to say, more voxels than usual were included into the infarct by the neurologist. That was why the proposed method had the lowest SI (74.359%) and the lowest sensitivity (59.444%) for patient number 14 among all the patients. These cases demonstrated the inconsistency in the semiautomatic way of segmentation. In contrast, the proposed method attained consistent segmentation, because it used the same value as the threshold in all cases.
Although the performance of the proposed method has been satisfactory, an even higher SI is still desirable. At present, the magnetic inhomogeneity does pose a limitation in further increasing the accuracy of infarct detection. To attain higher SI, it will be necessary to find a smarter way than the method we use in Step 8 to identify artifacts due to magnetic inhomogeneity.
This proposed infarct detection method will also be useful for the development of the automatic detection of white matter lesions. It can increase the accuracy of white matter lesion detection by excluding infarct lesions, which could be mistook for white matter lesions easily.
5. Conclusion
This proposed algorithm for acute infarct segmentation provides a prompt calculation of acute infarct volume from the DWI and ADC map. The infarct detection was achieved with fuzzy clustering that divided the DWI into ordered clusters based on voxel intensities. Under careful scrutiny on the intensity spectrum, candidate clusters were skimmed and deceptive labels further eliminated. Additionally, the ADC map was used to help identify spurious infarcts that were actually artifacts caused by magnetic inhomogeneity. This method attained high similarity indices. With the semiautomatic segmentation by the experienced neurologist for comparison, this automatic method attained high similarity indices.
Our future research will emphasize finding a more effective way to deal with the magnetic inhomogeneity to attain a higher accuracy in infarct segmentation. Applying the proposed infarct detection method to facilitate accurate computer-assisted segmentation of white matter lesions will also be our future effort.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The study is supported by a grant (NCU/LSH-101-A-005) from 2012 Annual National Central University and Landseed Hospital Joint R & D Project, Taiwan and in part by a grant (NSC 101-2221-E-008-015-MY3) from the National Science Council, Taiwan.
Supplementary Materials
Supplementary Table: The statistics of computer-assisted segmentation of cerebral infarct of individual patient after repeating the proposed method for 10 times. The values of mean±standard deviation of each statistical method for individual patient were shown to demonstrate the consistence of results of our proposed algorithm.