Abstract
Automated detection and diagnosis of small lesions in breast MRI represents a challenge for the traditional computer-aided diagnosis (CAD) systems. The goal of the present research was to compare and determine the optimal feature sets describing the morphology and the enhancement kinetic features for a set of small lesions and to determine their diagnostic performance. For each of the small lesions, we extracted morphological and dynamical features describing both global and local shape, and kinetics behavior. In this paper, we compare the performance of each extracted feature set for the differential diagnosis of enhancing lesions in breast MRI. Based on several simulation results, we determined the optimal feature number and tested different classification techniques. The results suggest that the computerized analysis system based on spatiotemporal features has the potential to increase the diagnostic accuracy of MRI mammography for small lesions and can be used as a basis for computer-aided diagnosis of breast cancer with MR mammography.
1. Introduction
Breastcancer is one of the most common cancers among women. Contrast-enhanced MR imaging of the breast was reported to be a highly sensitive method for the detection of invasive breast cancer [1]. Different investigators described that certain dynamic signal intensity (SI) characteristics (rapid and intense contrast enhancement followed by a wash out phase) obtained in dynamic studies are a strong indicator for malignancy [2]. Morphologic criteria have also been identified as valuable diagnostic tools [3]. Recently, combinations of different dynamic and morphologic characteristics have been reported [4] that can reach diagnostic sensitivities up to 97 and specificities up to 76.5.
As an important aspect remains the fact that many of these techniques were applied on a database of predominantly tumors of a size larger than 2 cm. In these cases, MRI reaches a very high sensitivity in the detection of invasive breast cancer due to both morphological criteria as well as characteristic time-signal intensity curves. However, the value of dynamic MRI and of automatic identification and classification of characteristic kinetic curves is not well established in small lesions when clinical findings, mammography, and ultrasound are unclear. Recent clinical research has shown that DCIS with small invasive carcinoma can be adequately visualized in MRI [5] and that MRI provides an accurate estimation of invasive breast cancer tumor size, especially in tumors of 2 cm or smaller [6].
Visual assessment of morphological properties is highly interobserver variable [7], while automated computation of features leads to more reproducible indices and thus to a more standardized and objective diagnosis. In this sense, we present novel mathematical descriptors for both morphology and dynamics and will compare their performance regarding small lesion classification based on novel feature selection algorithms.
More than 40 of the false-negative MR diagnosis are associated with pure ductal carcinoma in situ (DCIS) and with small lesion size and thus indicating a lower sensitivity of MRI for these cases. It has been shown that double reading achieves a higher sensitivity but is time-consuming and as an alternative a computer-assisted system was suggested [8]. The success of CAD in conventional X-ray mammography [9–13] motivates further the research of similar automated diagnosis techniques in breast MRI.
In the present study, we design and evaluate a computerized analysis system for the diagnosis of small breast masses with an average diameter of <1 cm.
The automated evaluation is a multistep system which includes global and local features such as shape descriptors, dynamical features, and spatiotemporal features combining both morphology and dynamics aspects. Different classification techniques are employed to test the performance of the complete system. Summarizing, in the present paper, a multifactorial protocol, including image registration, and morphologic and dynamic criteria are evaluated in predominantly small lesions of 1.0 cm or less as shown in Figure 1.
2. Material and Methods
2.1. Patients
A total of 40 patients, all females having an age range 42–73, with indeterminate small mammographic breast lesions were examined. All patients were consecutively selected after clinical examinations, mammography in standard projections (craniocaudal and oblique mediolateral projections) and ultrasound. Only lesions BIRADS 3 and 4 were selected where at least one of the following criteria was present: nonpalpable lesion, previous surgery with intense scarring, or location difficult for biopsy (close to chest wall). All patients had histopathologically confirmed diagnosis from needle aspiration/excision biopsy and surgical removal. Breast cancer was diagnosed in 17 out of the total 31 cases. The average size of both benign and malignant tumors was less than 1.1 cm.
2.2. MR Imaging
MRI was performed with a 1.5 T system (Magnetom Vision, Siemens, Erlangen, Germany) with two different protocols equipped with a dedicated surface coil to enable simultaneous imaging of both breasts. The patients were placed in a prone position. First, transversal images were acquired with a STIR (short TI inversion recovery) sequence (TR = 5600 ms, TE = 60 ms, FA = 90°, IT = 150 ms, matrix size 256 256 pixels, slice thickness 4 mm). Then a dynamic T1 weighted gradient echo sequence (3D fast low angle shot sequence) was performed (TR = 12 ms, TE = 5 ms, FA = 25°) in transversal slice orientation with a matrix size of 256 256 pixels and an effective slice thickness of 4 mm or 2 mm.
The dynamic study consisted of 6 measurements with an interval of 83 s. The first frame was acquired before injection of paramagnetic contrast agent (gadopentetate dimeglumine, 0.1 mmol/kg body weight, Magnevist, Schering, Berlin, Germany) immediately followed by the 5 other measurements. The initial localization of suspicious breast lesions was performed by computing difference images, that is, subtracting the image data of the first from the fourth acquisition. As a preprocessing step to clustering, each raw gray level time-series was transformed into a signal time-series of relative signal enhancement for each voxel, the precontrast scan at serving as reference, in other words . Thus, we ensure that the proposed method is less sensitive to changing between different MR scanners and/or protocols.
Automatic motion correction represents an important prerequisite to a correct automated small lesion evaluation [14]. Especially for small lesions, the assumption of correct spatial alignment often leads to misinterpretation of the diagnostic significance of enhancing lesions [15]. Therefore, we performed an elastic image registration method based on the optical flow method [16]. The employed motion compensation algorithm is based on the Horn and Schunck method [17] and represents a variational method for computing the displacement field, the so-called optical flow, in an image sequence. In contrast to optical flow, we do not want to compute the displacement field in a projected image of our data, but the actual displacement in 3D space. In our work, however, we favor the original quadratic formulation, since we explicitly need the filling-in effect of a nonrobust regularizer to fill in the information in masked regions. To overcome the problem of having a nonconvex energy in the energy functional, we use the coarse-to-fine warping scheme detailed in [16], which linearizes the data term as in [17] and computes incremental solutions on different image scales.
We tested motion compensation for two and three directions and found the optimal motion compensation results in two directions [18]. Segmentation of the tumor is semiautomatic and we define an ROI including all voxels of a lesion with an initial contrast enhancement of ≥50%. The center of the lesion was interactively marked on one slice of the subtraction images and then a region growing algorithm included all adjacent contrast-enhancing voxels and also those from neighboring slices. Thus a 3-D form of the lesion was determined. An interactive ROI was necessary whenever the lesion was connected with diffuse contrast enhancement, as it is the case in mastopathic tissue.
3. Computer-Aided Diagnosis (CAD) System
The small lesion evaluation is based on a multi-step system that includes a reduction of motion artifacts based on a novel nonrigid registration method, an extraction of morphologic features, dynamic enhancement patterns as well as mixed features for diagnostic feature selection and performance of lesion evaluation. Figure 1 visualizes the proposed automated system for small lesion detection.
3.1. Feature Extraction
The complexity of the spatio-temporal tumor representation requires specific morphology and/or kinetic descriptors. We analyzed geometric and Krawtchouk moments and geometrical features as shape descriptors, provided a temporal enhancement modeling for kinetic feature extraction and the scaling index method for the simultaneous morphological and dynamics representation.
3.1.1. Contour Features
To represent the shape of the tumor contour, the tumor voxels having nontumor voxel as a neighbor were extracted to represent the contour of the tumor. In this context, neighbor voxels include diagonally adjacent voxels, but not voxels from a different transverse slice. Due to the different grid sizes in the three directions of the MR images and possible gaps between transverse slices, the tumor contour in one transverse slice does not necessarily continue smoothly into the next transverse slice. Considering tumor contours between transverse slices therefore introduces contour voxels that are completely in the tumor interior in one slice. This is illustrated in Figure 2: the dark voxels are contour voxels and the arrows indicate the computed contour chain. If voxels in the tumor having at least one non-tumor voxel as a neighbor on an adjacent transverse slice were considered part of the contour, in this example, the crossed-out voxels would belong to the contour.
Figure 3 shows an example for a tumor where the contour shifts considerably from one transverse slice to another.
The contour in each slice was stored as an 1D chain of the 3D position of each contour voxel, constituting a “walk” along the contour. The chains of several slices were spliced together end to end to form a chain of 3D vectors representing the contour of the tumor.
Next, the center of mass of the tumor was computed as where is the number of voxels belonging to the tumor, and is location of the th tumor voxel. Since the center of mass was computed from the binary image of the tumor, irregularities in the voxel gray values of the tumor were not taken into account.
Knowing the center of mass, for each contour voxel , the radius and the azimuth (i.e., the angle between the vector from the center of mass to the voxel and the sagittal plane) were computed the following way: where the subscripts and denote the position of the voxel in sagittal and coronal direction, respectively. was also extended to the range from to by taking into account the sign of .
From the chain of floating point values , the minimum value and the maximum value were computed, as well as
The entropy was computed from the normalized distribution of the values into 100 “buckets”, where is defined as follows:
From the radius, , , , , and were used as morphological features of the tumor. From the azimuth, only the entropy (computed for as in (5) and (6)) was used as a feature, since the values and are always around and , respectively, and the value is not invariant under rotation of the tumor image.
An additional measurement describing the compactness of the tumor, which was also used as a feature, is the number of contour voxels, divided by the number of all voxels belonging to the tumor.
3.1.2. Morphological Features
The spatial and morphological variations of a tumor can be easily captured by shape descriptors. We analyze two modalities as shape descriptors based on moments: the geometric and Krawtchouk moments.
Geometrical Moments
We will employ low-order three-dimensional geometrical moment invariants as described in [19] because they have a low computation time and the results are stable to noise and distortion. We will utilize the 6 low-order finite-term three-dimensional moment invariants as described in [19]. There are one second-order and fourth-order, two third-order and three fourth-order moment invariants.
Krawtchouk Moments
Global and local shape description represents an important field in 3D medical image analysis. For breast lesion classification, there is a stringent need to describe properly the huge data volumes stemming from 3D images by a small set of parameters which captures the morphology (shape) well. However, very few techniques have been proposed for both global and local shape description. We employed Krawtchouk moments [20] as shape descriptors for both malignant and benign lesions. Weighted 3-D Krawtchouk moments have several advantages compared to other known methods: they are defined in the discrete field and thus do not introduce any discretization error like Spherical Harmonics defined in a continuous field and low-order moments can capture abrupt changes in the shape of an object. The weighted 3D Krawtchouk moments [20] form a very compact descriptor of a tumor, achieved in a very short computational time. Every tumor can be represented by Krawtchouk moments since it is expressed as a function in a discrete grid .
Krawtchouk moments represent a set of orthonormal polynomials associated with the binomial distribution [21]. The th order Krawtchouk classical polynomials can be expressed as a hypergeometric function: with and the hypergeometric function is defined as and with being the Pochhammer symbol The set of the Krawtchouk polynomials has elements. This corresponds to a set of discrete basis functions with the weight function
We assume that is a 3-dimensional function defined in a discrete field . The weighted three-dimensional moments of order of are given as with . Local features can be extracted by the appropriate selection of low-order Krawtchouk moments. is given as and are defined correspondingly. Thus, every 3-dimensional function in a 3-dimensional field can be decomposed into weighted 3-dimensional Krawtchouk moments .
The tumor can be represented by Krawtchouk moments since it is expressed as a function in a discrete space .
3.1.3. Dynamical Features
Lesion differential diagnosis in dynamic protocols is based on the assumption that benign and malignant lesions exhibit different enhancement kinetics. In [2], it was shown that the shape of the time-signal intensity curve represents an important criterion in differentiating benign and malignant enhancing lesions in dynamic breast MR imaging. The results indicate that the enhancement kinetics, as represented by the time-signal intensity curves visualized in Figure 4, differ significantly for benign and malignant enhancing lesions and thus represent a basis for differential diagnosis. In breast cancer, plateau or washout-time courses (type II or III) prevail. Steadily progressive signal intensity time courses (type I) are exhibited by benign enhancing lesions. Also, these enhancement kinetics are not only present in benign tumors but also in fybrocystic changes [2].
Computing the average signal intensity of the tumor before contrast agent administration (SI) and after contrast agent administration (), the relative enhancement can be computed as
To capture the slope of the curve of relative signal intensity enhancement (RSIE) versus time in the late postcontrast time, we computed the line that approximates the curve of the RSIE for the last three time scans. The values and are the least square solutions of the overdetermined system of equations for the three last points in time (), as well as for all tumor voxels , with being the RSIE in voxel at time scan .
The solutions to these equations are given by where is the number of voxels in the tumor, and is an abbreviation for . The slope was used as a feature to describe the dynamics.
3.1.4. Simultaneous Morphology and Dynamics Representations
The scaling index method [22] is a technique that is based on both morphology and kinetics. It represents the local structure around a given point. In the context of breast MRI, such a point consists of the sagittal, coronal, and transverse position of a tumor voxel and its third time scan gray value, and the scaling index serves as an approximation of the dimension of local point distributions.
Mathematically, the scaling index represents the 2-D image as a set of points in a three-dimensional state space defined by the coordinates and the gray value . For every point with coordinates the number of points in a sphere with radius and a sphere with radius is determined and the scaling index is computed based on the following equation: where is the number of points located within an -dimensional sphere of radius centered at . As radii, we choose the bounds of the tumor shape. Thus, the obtained scaling-index is a measure for the local dimensionality of the tumor and thus quantifies its morphological and dynamical features. There is a correlation between the scaling index and the structural nature: for clumpy structures, for points embedded in straight lines, and for points in a flat distribution.
For each of the three time scans , the standard deviation and entropy were determined and used as a feature to capture the heterogeneous behavior of the enhancement in a tumor.
3.2. Classification Techniques
The following section gives a description of classification methods applied to evaluate the effect of spatiotemporal features in breast MR images.
Discriminant analysis represents an important area of multivariate statistics and finds a wide application in medical imaging problems. The most known approaches are linear discriminant analysis (LDA), quadratic discriminant analysis (QDA), and Fisher’s canonical discriminant analysis.
Let us assume that describes a -dimensional feature vector that is, there are classes and there are samples available in group . The mean in group is given by and the covariance matrix is given by .
3.2.1. Bayes Classification Based on LDA and QDA
The Bayes classification [23] is based on estimating the prior probabilities for each class which describe the prior estimates about how probable a class is.
This classification method assigns each new sample to the group with the highest a posterior probability. Thus, the classification rule becomes where represent the means of the classes and the corresponding covariance matrix. The assignment to a certain class for a certain pattern is made based on the smallest determined value of .
There are two cases to be distinguished regarding the covariance matrices: if the covariance matrices are different for each class, then we have a QDA (quadratic discriminant analysis) classifier, while if they are identical for the different classes, it becomes an LDA (linear discriminant analysis) classifier.
3.2.2. Fisher’s Linear Discriminant Analysis
The underlying idea of Fisher’s linear discriminant analysis (FLDA) is to determine the directions in the multivariate space which allow the best discrimination between the sample classes. FLDA is based on a common covariance estimate and finds the most dominant direction and afterwards searches for “orthogonal” directions with the same property. The technique can extract at most components.
This technique identifies the first discriminating component based on finding the vector that maximizes the discrimination index given as with denoting the interclass sum-of-squares matrix and the intraclass sum-of-squares matrix.
4. Results
In the following, we will explore the results of the previously described features’ sets from different classification techniques. The results will elucidate the descriptive power of several tumor features for small lesion detection and diagnosis.
4.1. Effectiveness of Krawtchouk Moments
The Krawtchouk moments describe a representation of local shape parameters and can thus describe the differences in morphology between benign and malignant tumors. Since the obtained number of Krawtchouk moments is very high (>200), we reduced their dimension based on principal component analysis (PCA). Table 1 shows the results for the Krawtchouk moments for different classifiers and number of principal components. In general, the quadratic discriminant analysis shows the best results and for PC >10 they tend to deteriorate.
4.2. Effectiveness of Combined Feature Groups
We now examine not anymore every single feature but group the features together in specific classes that contain the features described in the previous sections. Table 2 shows the results for five distinct classifiers assuming motion compensation in 2 directions. The Krawtchouk moments (reduced to a six-dimensional vector by PCA) yield the best results since they capture both local and global shape properties.
We perform receiver operating characteristic (ROC) analysis to determine the sensitivity, specificity, and area under the curve (AUC) of the CAD system. The results of the sensitivity and specificity for the current data set based on specific features selected based on their discrimination capability and also in combination are shown in Table 3. The scaling index entropy yields the highest sensitivity and the 5th geometric moment the highest specificity. This finding is not surprising since the scaling index is a spatio-temporal feature while the geometric moment is averaging over the tumor’s shape. Since benign lesions tend to have smoother surfaces than malignant, this feature can be used as a first-step discriminator between those lesions. The inclusion of geometric moments in the feature set increases the sensitivity but leaves the specificity unchanged.
The best AUC-values for single features as well as for all features combined can be found in Table 4.
The AUC-values demonstrate that the contour features are very powerful descriptors and are able to capture the spatio-temporal behavior of small lesions.
5. Conclusion
The goal of the presented study was the introduction of new techniques for the automatic evaluation of dynamic MR mammography in small lesions and is motivated to increase specificity in MRI and thus improve the quality of breast MRI postprocessing, reduce the number of missed or misinterpreted cases leading to false-negative diagnosis.
Several novel lesion descriptors such as morphological, kinetic and spatio-temporal are applied and evaluated in context with benign and malignant lesion discrimination. Different classification techniques were applied to the classification of the lesions. A surprisingly low number of eight features proved to contain relevant information and achieved for both Fisher’s LDA and LDA good classification results. Krawtchouk moments proved to capture both the local and global shape features and represent thus in term of classification the best shape descriptors. In terms of spatio-temporal features, the scaling index entropy yields the highest sensitivity demonstrating that the enhancement pattern in small lesions has to be analyzed both in terms of spatial and temporal information. The benign characteristics are best described by geometric moments. The AUC-values demonstrate that the contour features can capture very well the spatio-temporal behavior of these small lesions.
The results suggest that quantitative diagnostic features can be employed for developing automated CAD for small lesions to achieve a high detection and diagnosis performance. The performed ROC-analysis shows the potential of increasing the diagnostic accuracy of MR mammography by improving the sensitivity without reduction of specificity for the data sets examined.
Acknowledgments
The research was supported by NIH Grant 5K25CA106799-05 and by an Alexander von Humboldt Fellowship.