Abstract
Background
The Greylevel Cooccurrence Matrix method (COM) is one of the most promising methods used in Texture Analysis of Magnetic Resonance Images. This method provides statistical information about the spatial distribution of greylevels in the image which can be used for classification of different tissue regions. Optimizing the size and complexity of the COM has the potential to enhance the reliability of Texture Analysis results. In this paper we investigate the effect of matrix size and calculation approach on the ability of COM to discriminate between peritumoral white matter and other white matter regions.
Method
MR images were obtained from patients with histologically confirmed brain glioblastoma using MRI at 3T giving isotropic resolution of 1 mm^{3}. Three Regions of Interest (ROI) were outlined in visually normal white matter on three image slices based on relative distance from the tumor: one peritumoral white matter region and two distant white matter regions on both hemispheres. Volumes of Interest (VOI) were composed from the three slices. Two different calculation approaches for COM were used: i) Classical approach (CCOM) on each individual ROI, and ii) Three Dimensional approach (3DCOM) calculated on VOIs. For, each calculation approach five dynamic ranges (number of greylevels N) were investigated (N = 16, 32, 64, 128, and 256).
Results
Classification showed that peritumoral white matter always represents a homogenous class, separate from other white matter, regardless of the value of N or the calculation approach used. The best test measures (sensitivity and specificity) for average CCOM were obtained for N = 128. These measures were also optimal for 3DCOM with N = 128, which additionally showed a balanced tradeoff between the measures.
Conclusion
We conclude that the dynamic range used for COM calculation significantly influences the classification results for identical samples. In order to obtain more reliable classification results with COM, the dynamic range must be optimized to avoid too small or sparse matrices. Larger dynamic ranges for COM calculations do not necessarily give better texture results; they might increase the computation costs and limit the method performance.
Background
Automated statistical and structural methods applied to digital medical images have shown that the amount of quantitative information available in the image exceeds the capacity of the human visual system [1]. These methods assume that image greylevel relationships and spatial distribution are directly influenced by the properties of the underlying tissue, which themselves, are dynamic and depend on biological and chemical composition. Therefore, minor image modifications can be quantified and monitored by appropriate methods even before they are perceivable by the human eye. Such automated methods are collectively known as Texture Analysis (TA) [1].
Since the physical properties of tissues are the basis for operating imaging modalities, the reliability of the imaging output depends on the ability of the modality to provide contrast between different tissues as well as local contrast that shows early changes in the physicalchemical properties within that tissue. MRI is known to provide the best image contrast among the imaging modalities available so far; therefore, MR images are believed to be rich in digital information that can be exploited by TA and would be of important analytical and diagnostic utility. In recent years, Texture Analysis on Magnetic Resonance Images (MRITA) has been applied successfully in clinical and experimental studies and is regarded as a reliable noninvasive tool of investigation, which combines the high contrast of MRI with the good sensitivity and specificity of TA. The quantitative texture data obtained from TA are relative rather than absolute; therefore, MRITA usually has to be followed by a standard classification method.
It has been demonstrated with laboratory animals that invivo MRITA of muscles correlates with histology during degeneration and regeneration processes [2]. Direct relationships between muscle contents of fat and collagen were found using texture classification on high resolution MRI [3]. MRITA has been clinically investigated on several tissues such as breast lesions [4], and hepatic fibrosis [5]. Brain tissue also has been studied using MRITA [68]. These studies recommended MRITA as a potential tool for noninvasive investigations of cerebral tumors as well as for healthy white and grey matter. In a previous work on brain gliomas, we investigated peritumoral white matter in regions defined by the radiologists as normal nonpathological tissues, but which were in the proximity of visible tumor margins. MRITA classified these regions as a homogenous texture class, separate from the other white matter regions which clustered in one broad class [8]. We suggested that this different texture could be due to invisible proliferation by tumoral cells [8].
Since TA is based on calculations with image greylevels, it becomes crucial to understand the impact of changing image properties on the stability of texture results and on the performance of the method [9]. Investigations of such relationships eventually will lead to an optimized and more reliable method for biomedical image analysis applications, which will require less processing time and less extensive calculations.
It is commonly assumed in TA applications that increasing the image dynamic range on which texture is evaluated improves textural feature representation; and consequently, gives better classification results. However, there is no evidence in biomedical image analysis literature to confirm or reject this assumption. The objective of the current work is to investigate the dependence of a commonly used TA method, the Cooccurrence Matrix (COM), on image dynamic range and matrix calculation approach for classifying white matter regions.
Methods
Patients and MRI data
In agreement with the French ethical legislation on clinical trials, whole brain MRI datasets were acquired in the sagittal plane for ten Glioblastoma patients (age = 53 ± 18; histologically confirmed by biopsy) using a Philips 3T Achieva MR system (Philips Medical System, Best, Netherlands). The imaging sequence used was ThreeDimensional Gradient Echo (TR = 9.87 ms, TE = 4.56 ms, flip angle = 8°). Field of View (FOV) = 256 mm × 256 mm, matrix size = 256 × 256, and a slice thickness of 1 mm, gave isotropic voxel resolution of 1 mm^{3}. Transversal sections were reconstructed from the original sagittal plane. Imaging procedures and clinical diagnosis were performed in Rennes University Hospitals, Rennes, France.
Each patient showed a tumor mass developed within the brain white matter. Three Regions of Interest (ROI)s were manually outlined in the normal white matter by a radiologist on a first transversal image Slice (S1) according to relative distance of the region to the tumor: one Peritumoral White matter (PtWm) close to the visible margins of the tumor; and two Distant White matter (DWm) taken far from the tumor on both hemispheres (figure 1). Each ROI contains of about 100–200 pixels. Volumes of Interest (VOI)s were constructed by copying the ROI position to the next two adjacent transversal slices (S2 and S3) producing volumes of about 300–600 voxels each. The VOI boundaries were inspected carefully to avoid overlapping structures. Only for one patient the location of the tumor did not allow for outlining a PtWM. A total number of 89 ROIs and 29 VOIs were available for this study.
Figure 1. MR image of brain glioblastoma and the surrounding white matter. Transversal slice of MRI of brain glioblastoma showing Tumor, T; and the normal white matter regions (solid lines): PtWm, Peritumoral; and DWm, Distant White matter. An ROI and the corresponding matrix are linked (red dashed lines) to illustrate the rescaling process. Matrix A represents the original ROI which has a dynamic range from 0 to 255 greylevels. The matrix B shows the same ROI after multiplying each pixel with the ratio of the maximum greylevel value allowed in B (31 in this case) to the actual maximum greylevel value of A.
Cooccurrence Matrices
The Cooccurrence Matrix (COM) was first introduced by Haralick [10] along with 14 derived features; most of them quantitatively describe image homogeneity and greylevels correlations. Some COM features have been found to be discriminative, and therefore, can be used for texture classification [10]. In a digital image, the number of bitsperpixel (bp) coding determines the maximum number of greylevels (N) in the image (2^{bp }= N). Hence, the allowed dynamic range of greylevels is from 0 to (N1).
The Classical approach of COM calculation (CCOM) samples the probability density function P_{d,θ}(i, j), which gives the probability of finding the two greylevels i and j at a distance d (d = 1,2,3,...) in the direction of angle θ_{xy }(θ_{xy }= 0°, 45°, 90°, and 135°), on a two dimensional image defined on the x and y axes. This calculation approach ignores useful spatial information that can be obtained from relationships between slices. Therefore, recent approaches try to maximize the usefulness of COM by including data at various angles on the zaxis. One of these approaches is known as Three Dimensional Cooccurrence Matrix (3DCOM) [8]. 3DCOM is calculated on image volumes composed of several adjacent slices, and involves nine angles on the zaxis (θ_{z }= 0°, 45°, 90°, 135°, 180°, 225°, 270°, 315°, and colinear) in addition to angles θ_{xy}. More details can be found in [8]. A Direction Independent (DI) matrix results from summing COM over all angles. This indicated in the notation below by θ = DI.
In this work, both approaches are calculated: i) CCOM on ROIs of the three adjacent slices (S1, S2, S3) giving three matrices CCOMS1, CCOMS2, and CCOMS3, respectively; and ii) 3DCOM on the VOI given by the three slices. For both approaches, the resulting matrix is always symmetric about its diagonal and of NxN size with N^{2 }number of entries. Five parameters were calculated from each matrix: Angular Second Moment, Inverse Difference Moment, Entropy, Contrast and Correlation [10]. These five parameters were selected because they were found to be good descriptors of white matter texture in a previous work [8]. They provide the main information about image homogeneity and the existence of correlated patterns in the image.
The original MR images are usually digitized over 16 bitsperpixels (65536 greylevels). It is computationally extensive and time consuming to calculate COM over such a large dynamic range. Therefore, it is a standard procedure in medical image analysis to apply a quantization process in order to reduce the original range to a userdefined value of N. This is done by scaling the original pixel values with the ratio between the maximum greylevel allowed in the rescaled image and the actual maximum greylevel in the original image (figure 1). Prior to COM calculations, each ROI is rescaled for five different values of N, (N = 16, 32, 64, 128, and 256). All texture calculations and image processing methods were implemented using Matlab^{® }(ver 7.0, MathWorks Inc., Natick, MA, USA), on a PC with Intel^{® }Pentium^{® }4.0 processor and 1.24 Gb RAM.
Features Selection and Classification
Features selection aims to identify the most discriminating parameters from each matrix that separate the different classes most efficiently. Fishercoefficient (Fcoefficient) was calculated for this purpose, giving the ratio of between class variance to within class variance [11] for each parameter. The ten parameters of the highest Fcoefficient were entered to Linear Discriminant Analysis (LDA) for classification. LDA aims to find a linear transform matrix such that the ratio of withinclass scatter matrix to betweenclass scatter matrix is maximized. Such a transform is composed of eigenvectors corresponding to the largest eigenvalues of this ratio of matrices; more details about the classification method can be found in [12]. Cross validation was performed using "leaveoneout" criterion, which works by leaving one observation (i.e one ROI) out of the classifier each time the classification model is recalculated and then project this observation into the model to test its validity. This process is carried out for all observations. The percentages of False Negatives (FN) and False Positives (FP) were evaluated. The Receiver Operating Characteristic (ROC) curve was analyzed, which represents the relationship between the '1Specificity' and the 'Sensitivity' of the test. The Area Under the ROC Curve (AUC) is used to judge the separability of the two classes for the given dataset and classifier. An AUC of 1.0 represents a perfect classifier, while an AUC of 0.5 represents a random classifier.
Features Selection was performed using B11 software (version 3.2, ^{©}1999–2002 by Michal Strzelecki), which is developed under the auspices of COST action B11 European project [12]. Linear Discriminant Analysis (LDA) was followed by cross validation and was performed using the software Minitab 15 (^{© }2007 Minitab Inc). The ROC curve was analyzed and AUC were calculated using SPSS 15.0 (^{© }1989–2006 SPSS Inc).
Results and discussion
LDA classification on PtWm and DWm white matter regions always separated PtWm into a distinct homogenous class. This class was well distinguished for small as well as for large dynamic ranges for all matrices. However, the number of classification errors between the two classes depended remarkably on the dynamic range along with COM approach used. Table 1 represents the percentage of False Negatives (FN: PtWm classified as DWm) and False Positives (FP: DWm classified as PtWm) for each number of greylevels N and matrix calculation approach. The average errors and standard deviation (Mean ± SD) for CCOMS1, CCOMS2, and CCOMS3 over the three slices is also presented and will be used for comparison with 3DCOM.
Table 1. Classification results using crossvalidated LDA and for Peritumoral White matter (PtWm) classified as Distant White matter (DWm) (False Negative: FN).
Analyzing table 1 shows that the Mean FN or FP of CCOM (S1,S2,S3) decreases progressively with increasing N until reaching N = 256 for which it increases again (table 1). For 3DCOM, the lowest value of FN occurs at N = 128, which was less than those obtained from Mean CCOM for any other N. The most discriminating parameters for this analysis and their FCoefficients are presented in table 2. The percentage of FN shows a considerable increase when 3DCOM is calculated for N = 256; however, FP represents the lowest percentage obtained (table 1).
Table 2. The ten most discriminating parameters, according to the Fisher (F) coefficient, between the two white matter classes (Peritumoral white matter and distant white matter).
The bar graph of test outcomes measures (sensitivity and specificity) (figure 2) demonstrates a balanced tradeoff between the sensitivity and specificity of the 3DCOM method at N = 128 and N = 32 with higher values at the former (figure 2a). This balance is lost at other values of N (figure 2b). The Mean CCOM method on the three slices (S1,S2,S3) shows the highest sensitivity and specificity at N = 128 (figure 2b). For either CCOM or 3DCOM, figure 2 demonstrates that the specificity of the method is always higher than its sensitivity. The Area Under the ROC Curve (AUC) represents a comprehensive measure for evaluating the accuracy of the classifier (table 1). By comparing AUCs of the Mean value of the three CCOMs and those of 3DCOM, it can be shown that the highest AUC value was obtained for 3DCOM at N = 128, while the lowest was obtained for Mean CCOMs at N = 16 (figures 3a and 3b, respectively). It can also be shown that the highest value of AUC among Mean CCOMs was obtained also at N = 128 (table 1).
Figure 2. Sensitivity and specificity bar graphs. Sensitivity and specificity bar graphs for (a) 3DCOM on white matter VOIs; and, (b) The Mean value of (CCOM) on the individual slices ROIs (S1, S2, and S3). CCOM: Two Dimensional Classical Cooccurrence Matrix. 3DCOM: Three Dimensional Cooccurrence Matrix. VOI: Volume of Interest. ROI: Region of Interest.
Figure 3. ROC curves showing the highest and lowest AUC. Receiver Operating Characteristic (ROC) curves showing: a) the highest Area Under the Curve (AUC) (= 0.895) which was obtained using 3DCOM at N = 128; and, b) the lowest AUC (= 0.715) obtained using the Mean CCOMs at N = 16. CCOM: Two Dimensional Classical Cooccurrence Matrix. 3DCOM: Three Dimensional Cooccurrence Matrix.
In this study, PtWm clustering as a separate white matter region is consistent with previous findings [8]; however, we demonstrate in the current work that classification accuracy is highly dependent on the dynamic range of image quantization for both COM calculation approaches (CCOM and 3DCOM). Also, we can see that classification results among different slices might give diverse results in spite of carrying out the analyses on identical positions. This can be demonstrated for FN at N = 16 that gave 22% on CCOMS1 and 55% on CCOMS2.
It can be also shown that calculating 3DCOM on small dynamic ranges (N = 12, 32 and 64) does not enhance classification as long as the dynamic range remains relatively small. In contrast, 3DCOM on a larger dynamic range (N = 128) enhances classification remarkably, but a further increase of N worsens the method's sensitivity. Although method's specificity has increased at N = 256, the tradeoff between sensitivity and specificity remains an important factor to take into account when evaluating the method's performance. Therefore, N = 256 is probably not a good choice for 3DCOM analysis.
The relationship between the dynamic range and classification accuracy can be related to COM characteristics. This matrix, by definition, is a probability density matrix of unit sum. Decreasing the dynamic range means that the original ROI will be reduced to smaller adjoining values of greylevels as shown in figure 1; therefore, cooccurrence matrix will be smaller and the joint probabilities will be estimated for a limited number of matrix entries (eg. N = 16, COM size = 16 × 16 = 256 entry). This could be insufficient to represent texture features and may result in higher classification errors. On the other hand, increasing the dynamic range will spread the greylevels over a larger scale producing matrices with sufficient number of entries; and then, discriminating texture features would have more chance to appear; consequently, this would reduce the percentage error. Further increase of N values produces sparse matrices with probabilities broken down over a huge number of entries (65535 for N = 256); in other words, feature representation would be weakened and classification errors increased. It merits to mention that the processing time for calculating 3DCOM using N = 128 was within a fraction of a second, while it took almost 30 seconds for calculating the same matrix using N = 256. The increase in processing time is even more significant for larger VOIs.
From these results, we recommend quantizing image ROI/VOI to a number of N = 128 greylevels prior to texture analysis of white matter. This value represents a compromise for applying cooccurrence matrix calculations in white matter texture studies for the two dimensional approach as well as for the three dimensional one.
Conclusion
In this work we have demonstrated that the dynamic range on which texture features are evaluated, particularly when using the cooccurrence matrix, can directly influence the accuracy of classification of white matter regions. We found that rescaling the ROI to a dynamic range of greylevels from 0 to 127 (i.e. N = 128) gives the best classification results using two dimensional cooccurrence matrix CCOM represented by the mean value of the three slices (S1, S2, and S3). It also gives the best balance between sensitivity and specificity, using the three dimensional cooccurrence matrix 3DCOM. For both types of matrices, the AUC of the ROC curve was maximum at N = 128. We conclude that a reduced userdefined dynamic range can be faster, computationally less extensive, and more efficient in separating texture classes.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
DMG has designed the study, implemented the texture analysis methods on Matlab^{® }7.0, acquired data, analyzed results, and drafted the manuscript. MKA has participated in the programming procedures. JDC has set and supervised the protocol of MR image acquisition in Rennes University Hospitals according to rules and regulations set by Ethics Committees. FMG has participated in results analysis and critical revision of the manuscript.
Acknowledgements
This work has been achieved in coordination with COST European project Action B21 "Physiological modelling of MR image formation". It has been presenting during COST B21 Meeting in Bled, Slovenia, 2007.
A part of this work has been funded by United Arab Emirates University, Research Grant number: 0102211/07.
The authors would like to thank Biatrice Carsin, for MRI acquisition (CHRU Pontchaillou), and PierreAntoine Eliat for technical assistance (Rennes I University).
References

Castellano G, Bonilha L, Li LM, Cenes F: Texture analysis of medical images.
Clinical Radiology 2004, 59:10611069. PubMed Abstract  Publisher Full Text

MahmoudGhoneim D, Cherel Y, Lemaire L, de Certaines JD, Maniere A: Texture analysis of Magnetic Resonance Images of rats muscles during atrophy and regeneration.
Magn Reson Imaging 2006, 24:167171. PubMed Abstract  Publisher Full Text

MahmoudGhoneim D, Bonny JM, Renou JP, de Certaines JD: Exvivo Magnetic Resonance Imaging Texture Analysis Can Discriminate Genotypic Origin in Bovine Meat.
J Sci Food Agr 2005, 85:629632. Publisher Full Text

Chen W, Giger ML, Li H, Bick U, Newstead GM: Volumetric texture analysis of breast lesions on contrastenhanced magnetic resonance images.
Magn Reson Med 2007, 58:56271. PubMed Abstract  Publisher Full Text

Kato H, Kanematsu M, Zhang X, Saio M, Kondo H, Goshima S, Fujita H: Computeraided diagnosis of hepatic fibrosis: preliminary evaluation of MRI texture analysis using the finite difference method and an artificial neural network.
AJR Am J Roentgenol 2007, 189:11722. PubMed Abstract  Publisher Full Text

Schad LR, Blüml S, Zuna I: MR tissue characterization of intracranial tumors by means of texture analysis.
Magn Reson Imaging 1993, 11:889896. PubMed Abstract  Publisher Full Text

HerlidouMême S, Constans JM, Carsin B, Olivier D, Eliat PA, NadalDesbarats L, Gondry C, Le Rumeur E, IdyPeretti I, de Certaines JD: MRI Texture Analysis on Texture Test Objects, Normal Brain and Intracranial Tumors.
Magn Reson Imaging 2003, 21:989993. PubMed Abstract  Publisher Full Text

MahmoudGhoneim D, Toussaint G, Constans JM, de Certaines JD: "Three Dimensional Texture Analysis in MRI: a preliminary evaluation in gliomas".
Magn Reson Imaging 2003, 21:983987. PubMed Abstract  Publisher Full Text

Collewet G, Strzeleck M, Mariette F: Influence of MRI acquisition protocols and image intensity normalization methods on texture classification.
Mag Reson Imaging 2004, 22:8191. Publisher Full Text

Haralick RM, Shanmugam K, Dinstein I: Textural Features for Image Classification.
IEEE T Syst Man Cy 1973, 3:610621. Publisher Full Text

Swets W: Using discriminant eigenfeatures for image retrieval.

Materka A: MaZda and B11 User's Manual ^{©}1999–2002. [http://www.eletel.p.lodz.pl/merchant/mazda/order1_en.epl] webcite
Prepublication history
The prepublication history for this paper can be accessed here: