SYSTEM AND METHOD FOR AUTOMATED SEGMENTATION, CHARACTERIZATION, AND CLASSIFICATION OF POSSIBLY MALIGNANT LESIONS AND STRATIFICATION

OF MALIGNANT TUMORS

The present application claims the benefit of priority from U.S. provisional application no. 61/029,697 the contents of which are incorporated herein by reference.

BACKGROUND OF THE INVENTION

Magnetic resonance imaging (MRI) of the breast is quickly becoming an important component of breast cancer screening in the United States. The American Cancer Society revised its MRI screening guidelines in 2007 to include women with a 20 to 25 percent increased lifetime risk, including family history of breast and ovarian cancers or a history of prior treatment for Hodgkin's disease. MRI has become a tool for breast cancer screening because of its advantages over mammography and ultrasound. First, breast MRI is more sensitive than x-ray mammography and ultrasound and thus may be used to detect lesions that would otherwise remain undetected. Second, MRI has been shown to be advantageous in screening women with dense breast tissue, which is common in younger patients. Third, and likely most importantly, the dynamic contrast enhanced (DCE) component of DCEMRI provides information about vascularity of a lesion that is more difficult to obtain using x-ray mammography or ultrasound.

Triple negative (TN) breast cancer has recently gained much attention in the field of breast cancer research. TN breast cancer is a molecular subtype that lacks expression of the estrogen receptor, progesterone receptor, and the HER2 receptor. Because of its lack of receptor expression, targeted therapies are ineffective, and chemotherapy is currently the only treatment available. TN breast cancer is also particularly aggressive and accounts for 12% to 26% of all breast cancers, most often occurring in young patients and African-American patients. The research community has become interested in TN breast cancer because of its particularly aggressive clinical course and lack of targeted therapies. Studies examining the pathological phenotype have revealed a heterogeneous group of breast cancers that often present as interval cancers, presenting in the months in between annual breast cancers screenings. This clearly suggests the need for greater exploration into the detection methods and biological understanding of TN breast cancers. Studies have shown that x-ray mammography is less effective in TN breast cancer screening than DCE-MRI due to the greater sensitivity of DCE-MRI in detecting TN and basal-like breast cancer phenotypes. A similar pattern of enhancement, described as rim enhancement, was observed in both subtypes. Although,

qualitative radiologic descriptors have been developed for small pilot studies in TN breast cancer cases, the inventors are unaware of any quantitative classification of the TN subtype.

SUMMARY OF THE INVENTION The present invention is embodied in a method and apparatus for classifying possibly malignant lesions from sets of DCE-MRI images. The method includes receiving a set of MRI slice images obtained at respectively different times, where each slice image includes voxels representative of at least one lesion. The images are processed to determine the boundaries of the lesions and the voxels within the identified boundaries in corresponding regions of the images from each time period are processed to extract kinetic texture features. The kinetic texture features are used to drive a classification process wherein each lesion is identified as as malignant or benign. The malignant lesions are further classified to separate TN lesions from non-TN lesions. BRIEF DESCRIPTION OF THE DRAWINGS

Fig. 1 is a block diagram of an exemplary DCE-MRI system suitable for use with an embodiment the subject invention.

Fig. 2 is a flow-chart diagram of an example method according to an embodiment of the invention. Fig 3 is a flow-chart diagram showing details of the contour definition step of the flow-chart diagram shown in Fig. 2.

DETAILED DESCRIPTION

While the embodiments of the subject invention described below concern the detection of breast tumors, it is contemplated that they may be applied generally to detecting and classifying possibly malignant lesions in other parts of the body based on DCE-MRI data. In addition, the invention may be applied to detect and classify non-malignant regions of interest (ROI) in a body.

Due to the clinical nature of TN tumors, however, accurate and consistent identification of these specific tumors is desirable, and a computer-aided diagnosis (CAD) system that could detect the TN radiologic phenotype would assist clinicians in therapeutic decision-making, monitoring therapy response, and increasing our understanding of this aggressive breast cancer subtype.

Breast DCE-MRI is performed by first injecting gadolinium diethylenetriamine-pentaacid (Gd-DTPA) into the patient's bloodstream and concurrently acquiring MRI images of the breast. Since malignant lesions tend to grow leaky blood vessels in abundance, the contrast agent is taken up by tumors preferentially and subsequently leaks out of tumors rapidly. This provides the use

of DCE-MRI with an advantage over the use of other modalities to detect the tumors and contributes to the high sensitivity rates reported for breast DCE-MRI.

Both benign and malignant neoplastic tissue, however, frequently have contrast enhancement patterns that differ from normal breast tissue, and these

5 abnormalities are highlighted in the time-dependent MRIs. As such, it may be difficult for radiologists to differentiate between benign and malignant lesions simply by observing the contrast enhanced lesion on the post-contrast MRI. For this reason, clinicians have explored various methods of observing and measuring the manner in which a lesion takes up the contrast dye. It was found, for example,

I ₀ that data in the temporal MRIs could be plotted as single data points on a time series curve that is reflective of the lesion type. It has been shown that malignant lesions have a characteristic curve, showing rapid uptake of contrast (steep positive initial slope) and rapid washout (subsequent negative slope). Benign lesions, on the other hand, have slow contrast uptake (small positive initial slope) and then is plateau or do not reach peak intensity during the time series. This description of the DCE-MRI data is now considered convention in radiologic interpretation of breast DCE-MRI.

Despite making rapid strides in the interpretation of breast DCE-MRI in the past decade, the optimal accuracy in diagnosis using breast DCE-MRI has not been

2o achieved, and only highly experienced radiologists are able to accurately interpret breast MRIs. Furthermore, inter-observer variability for radiologic interpretation tends to be high. Numerous studies in the field of breast MRI have shown efficacy and improved diagnosis rates using CAD, but breast MRI CAD has yet to achieve the accuracy of that seen in CAD for x-ray mammography. There is a growing need

2 ₅ in the field to provide multimodal data interpretation methods that will be able to consistently and accurately detect, diagnose and differentiate breast tumors in general and specifically benign and malignant tumors.

Described below is a comprehensive CAD system for the discrimination of (a) benign from malignant breast lesions, and (b) triple negative from non-triple

30 negative breast lesions. The example embodiments concern a similar understanding of the typical lesion enhancement patterns to create an Expectation Maximization-driven Active Contour scheme to automatically extract the lesion contour. Quantitative features are then automatically obtained for each lesion. Although features, such as morphology and texture are considered, the

3 ₅ example embodiment uses a DCE-MRI feature called kinetic texture, which characterizes spatio-temporal changes in lesion texture. A classifier, for example, a support vector machine may be used to quantitatively classify the breast lesions in

the dataset. Graph embedding or other nonlinear dimensionality reduction techniques may then be used to reduce data dimensionality and aid in visualization of the relationships between different breast lesion classes. The example CAD system employs three components: (1) Lesion detection and segmentation combining the time-series signal intensity data with an Expectation Maximization- driven Active Contour scheme; (2) Feature extraction using over 500 different features for the identification of the most discriminatory features of the tumor types; and (3) Breast lesion classification performed in a hierarchical manner, first distinguishing malignant from benign lesions and then, within the group of malignant lesions, identifying those of the aggressive triple-negative molecular phenotype.

An example system suitable for use with the subject invention is shown in Fig. 1. This system includes a magnetic resonance imaging system 110, a processor 112, a memory 114 and a display and keyboard unit 116. The processor 112 controls the MRI system 110 to capture multiple sets of images of a breast, each image being captured at a different time. The MRI images are DCE-MRI images so that at least one of the sets of images is taken at a time before the contrast agent has been taken-up by the tumor and at least two other sets of images are taken at times after the tumor has taken up the contrast agent. In the example embodiment, the memory 114 has sufficient space to hold all of the sets of images used by the algorithm described below with reference to Figs. 2 and 3. The display/keyboard may be used to control the algorithm, to process intermediate images, as described below and to display the results of the example method. In the example method, described below, a total of 41 (24 malignant, 17 benign) breast DCE-MRI studies were used as the training and testing dataset. Immunohistochemistry was performed on the 24 malignant lesions to measure estrogen receptor, progesterone receptor, and HER2 receptor status. Of these, 13 were determined to be TN cases and 11 non-TN. Sagittal Tl weighted, spoiled gradient echo sequences with fat suppression consisting of one series pre-contrast injection of Gd-DTPA and three to five series post-contrast injection were acquired (Matrix 384x384 512x512, or 896x896, slice thickness 3cm). Temporal resolution between post-contrast acquisitions was in the range of 45-90 seconds. For each study, all of the sets of images were stored in the memory 114 by the processor 112. The data for these studies may be entered into the memory 114 through an input/output (I/O) port of the processor 112 (not shown) or through a wired or wireless network connection (not shown) to the processor 112.

In the example embodiment, the 41 studies are used in a testing mode to determine which image features are suitable for discriminating between malignant and non-malignant lesions and, within the identified malignant images to distinguish between TN and non-TN lesions. The sample studies and other sample studies may be used to determine ranges of parameter values that can be used to classify lesions as possibly malignant with known probabilities. The system shown in Fig. 1 may then be used according to the algorithm shown in Figs. 2 and 3 to assist a healthcare professional in the analysis of a set of DCE-MRI images.

With reference to Fig. 2, the operation of the example embodiment in a diagnosis mode begins at step 210 when, for a particular patient, pre-contrast MRI data is obtained. Next, at step 212, the contrast agent is administered. While this step is shown as being after step 210, it is contemplated that the contrast agent may be administered before or concurrently with step 210 as long as the set of slice images in step 210 is obtained before the contrast agent is taken-up by the lesion.

Next, at step 214, the obtained data is stored in the memory 114. At step 216, the system waits for an amount of time, for example, 45 to 90 seconds. At step 218 the algorithm determines if the last set of slices has been obtained. If not, control transfers to step 214 to obtain another set of post-contrast MRI images. If, at step 218, the last set of MRI data has been obtained, control transfers to step 220.

At step 220, using the study data - or data from a particular individual when operating in the diagnosis assistance mode - for each set of images, a radiologist selects a lesion slice most representative of each lesion. This selection may be made from set of image slices taken at any time. It is contemplated, however, that the selection may be made from a time after the contrast agent has been taken-up by the lesion. In the example embodiment, the analysis described below is performed only for that slice of the lesion volume. That is to say, for the particular study, only the identified slice in the set of slices taken at each time is considered. It is contemplated, however that multiple slices may be used, with each slice being compared to corresponding slices in succeeding and/or preceding time frames. At step 222, the selected slices are processed, according to the example method, to automatically define the contours of the lesions in the selected slice for each time period. This example process, as shown in Fig. 3 includes two steps. First, at step 310, the data is processed according to an expectation/maximization algorithm to group the voxels and second, at step 312, the grouped voxels are processed using a magnetostatic active contour model, to define the contour, or

boundary, of the lesion in the image. The lesion so obtained is then visually inspected by a radiologist and if any obvious errors in segmentation are observed, these are corrected for by an attending radiologist with expertise in magnetic resonance mammography via a medical imaging software such as MRIcro, available from www.mricro.com. In the example embodiment, the medical imaging software runs on the processor and displays the result on the display of the keyboard/display unit 116. Since the technology could be applied to identifying suspicious regions in the breast that are not necessarily lesion-like, from hereon we refer to the result of segmentation as constituting a region of interest (ROI). One embodiment of the ROI is a lesion.

In step 310, the radiologist-selected slice at each time point pre- and post- contrast is compiled into a three-dimensional matrix. A signal enhancement curve is then generated for each voxel based on the signal intensity values at each time point. The coefficients of a third order polynomial are then obtained from fitting the time series curve of each individual voxel. The image scene is defined as c = (C, f), where C is a spatial grid of voxels c e C and f is the associated signal intensity at time t e {0, 1, 2, . . . , T - 1}. The time-signal intensity vector {f°(c), f*(c), I ² Cc), ■ ■ ■ , F ^-1 Cc)) for every voxel c in the MRI image is fitted to a third order polynomial, in a least-squares sense, which is described by equation (1). f'(c) = a ₂ ^c t ³ + a ₂ ^c t ² + a°t + a ₀ ^c (1)

The Expectation Maximization (EM) algorithm groups the voxels based on a time-series coefficient matrix, [ a _c - {al, a\,a\ , a ₀ ^c ),\/c e C ]. Based on Bayes' theorem, the EM algorithm aims to compute the posterior probability P _ck of each voxel c belonging to k e {1, 2, . . . ,K}, given the priori p _ck , where p _c ι< is the priori probability that voxel c belongs to class k, and K is the number of Gaussian mixtures. In the example embodiment, K is 4, the dimension of the row vector a _c . The algorithm is run iteratively, comprising two steps: the Expectation step (E- step) and the Maximization step (M-step). The E-step calculates the posterior probability P _ck based on the current parameters of Gaussian mixture model while the M-step recalculates or updates the model parameters, ∑ _k = {μ^ σ _k , β _k } where μ _k and σ _k are the mean and covariance of each Gaussian component, respectively, and the β _k values are mixture coefficients in the Gaussian mixture model. After a pre-defined number of iterations, voxel c is assigned to one of K classes, depending on which has the highest posterior probability P _ck . The EM results obtained in step 310 of Fig. 3 are then used to provide a specific initialization to an active contour model to obtain the final ROI segmentation. In step 312 of the example embodiment, the Magnetostatic Active

Contour (MAC) model as described in an article by X. Xie and M. Mirmehdi, entitled "MAC: Magnetostatic Active Contour Model," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 30, no. 4 pp. 632-646, April 2008, is employed to extract ROI boundaries. It is contemplated, however, that other active contour models may be used to define the ROI boundaries. The MAC model implements a bi-directional force field F generated from a hypothetical magnetostatic interaction between current carrying contours and object boundaries. The field F is defined over the entire image C and drives the contour towards the object boundary. The level set equation of the contour as proposed in the referenced article takes the form described in equation (2).

(2) dφ where φ is the level set function, —*- is the time derivative of φ, α is a real constant, dn q(C) = 1/(1 + I VC\ ) and V{-) represents the 2D gradient. For a given image C, K class likelihood scenes L _k = (C, l _k ) are constructed where l _k (c) assigns each voxel c e C the probability P _c κ of belonging to a class k determined from the EM result. L _k is maximized to obtain a binarized scene L ^B _k = (C,/f) where /f (c) e {0,1} and

Il (c) = 1 if P _ck (c) is the highest probability assigned to c. The appropriate scene L ^B _k representing the ROI class is manually selected and is used to initialize the active contour. The initialization of the contour is defined as circles centered at centroids of the objects detected in L ^B _k via connected component labeling. The contour is then evolved until the difference between the contours of the current iteration to the next is below an empirically determined threshold.

It is desirable to have the boundaries of the ROIs well defined in order to obtain a true measure of the texture of the ROI. If, for example, the boundary of the ROI were not well defined and included voxels that were not a part of the ROI, these extra voxels would affect the texture determination. The example method described above, however, provides well-defined boundaries and tend to exclude voxels surrounding the ROI.

Referring to Fig. 2, after ROI segmentation in step 222, over 500 morphologic, textural, and kinetic features, h _u , u e {1, 2, . . . , 508}, are extracted in step 224. These features include gradient features and first and second order statistical features. These features may be obtained at several times, both before and after the contrast agent has been taken-up by the tumor by applying linear and spatial operations to the voxels inside the boundaries of the ROIs in the

processed MRI slices. In the example that follows, time t=0 corresponds to the time at which the contrast agent has not been taken-up and all later times, (te{l, 2, ... ,T-1}) correspond to times after the tumors have begun to take up the contrast agent. Eleven non-steerable gradient features are obtained using Sobel, Kirsch and standard derivative operations. Gabor gradient operators comprising the steerable class of gradient features are also defined for every c ε C where c = (x, y), These features, h _u (c) are described by equation (3). where ω is the frequency of a sinusoidal plane wave along the X-axis, and ψ _x and ψ _γ are the space constraints of the Gaussian envelope along the X and Y directions respectively. Filter orientation, β, is affected by the coordinate transformations: x' = z(x cos θ + y sin θ) and y' = z(-x sin θ + y cos θ), where z is the scaling factor.

Tt TU Tt 7ϊ

Gabor gradient features were calculated at 6 scales ( z e { — τ=,— , — τ=,... — } ), 8

r _n π π 3π π 5π 3π Iπ _{λ V} orientations (θ e {0, — , — , — , — , — , — , — }) and 4 window sizes (s e {3, 5,

8, 15}) .

Four first order statistical features (mean, median, standard deviation, and range) for 3 different window sizes are calculated for the gray values of pixels within the sliding window neighborhood N ₅ , s e {3, 5, 8}. Thirteen Haralick features are also included in the extracted features. To calculate the second order statistical (Haralick) feature scenes, a GxG cooccurrence matrix O _d/C/S is computed, associated with 3Vj(q), where G is the maximum grayscale intensity in C. The value at any location [e _u e ₂ ] in O _d , _c ,s, where ei, e ₂ e {1, 2, . . . ,M}, represents the frequency with which two distinct voxels q, C _j e 5V _s (c) where i, j e{ l, 2, ... , \C\ } with associated image intensities f(q) = g _i# f(C _j )= g ₂ are separated by distance d. A total of 13 Haralick features including contrast energy, contrast inverse moment, contrast average, contrast variance, contrast entropy, intensity average, intensity variance, intensity entropy, entropy, energy, correlation, and 2 information measures are extracted at every voxel ceC, based on O _d , _c ,s, for s e {3, 5, 7}, d = 1 and G e{64, 128, 256}.

The feature set includes the extracted features described above as well as kinetic features calculated from the extracted features. The kinetic features are generated from the extracted features in step 226. Each of the kinetic features models the behavior of a voxel across all of the sample times in the set of

corresponding slices. Although calculation of kinetic features is described above with respect to a ROI, it is understood that kinetic features may also be calculated for one or more non-ROI areas, by using extracted features representing the non- ROI areas. For example, kinetic features of non-ROI areas may be used to characterize and/or quantify breast parenchyma.

Kinetic signal intensity features are computed in step 226, as the coefficients \a ₃ , a ₂ , a _lf a ₀ ] of a third order polynomial obtained from fitting a curve in a least-squares sense to the signal intensity contrast enhancement curves. Hence, for each c in C ₀ , Ci, C ₂ ,..., C _τ -i, a third order curve is fitted using equation (4) : f(t) = a ₃ t ³ +a ₂ t ² +ait+a ₀ (4)

Where t e{0, 1, 2, . . . , T-l} and c e C. Note that f(c) represents the signal intensity at each spatial location, c ε C, across the pre- and post-contrast MRI scenes. To calculate the kinetic textural features, a first order statistical descriptor which could be at least one of the textural feature's mean, mode, median, variance, or standard deviation value, p _u , is plotted over time such that a kinetic texture curve is created, which is analogous to the one created for signal intensity. A third order polynomial is fitted to this curve to characterize its shape, defining four associated coefficients as shown in equation (5):

Pu(t) = r _U) 3t ³ + r _u> 2t ² + r _u ,it + i\,o- ( ₅ )

[r _u , _3ι r _u , ₂ , r _u ,i, r _U)0 ] is the feature vector, describing the kinetic texture feature, for each texture feature, u, as a function of time. The Pre-contrast Textural Features are defined as the mean texture values for each feature described above before contrast injection (t=0).

A summary of the extracted features is shown in Table 1. Table 1: Extracted Features

The features listed above are used in the testing mode of the system in which the 41 known cases are used to train and test the example system. As described below, the morphological features may not be needed when the system is operated in diagnostic assistance mode. Indeed, it is contemplated that a suitable classifier can be constructed using only the kinetic texture features.

After the kinetic features have been calculated in step 226, the next step in the process is to classify the features at step 228. In the example system, support vector machine (SVM) methods are applied to evaluate the ability of each feature class (morphology, texture, kinetic texture and kinetic signal intensity) to classify each ROI as benign or malignant. It is contemplated, however, that other classification methods may be used, such as neural networks, Hidden Markov Models and Frequent Itemset Mining. An example Support Vector Machine algorithm is described in a publication by

C. Cortes and V. Vapnik entitled "Support Vector Networks," Machine Learning, vol 20 no. 2 Springer Netherlands 273-297 (1995). An example Hidden Markov Model system is described in a publication by S. Wong, A. B. Gardner, A. M. Krieger, B. Litt, "A Stochastic Framework for Evaluating Seizure Prediction Algorithms Using Hidden Markov Models," J. Neurophysiology 97(3): 2525-2532 (2007). An example Frequent Itemset Mining system is described in a publication by A. B. Gardner, A. M. Krieger, G. Vachtsevanos, B. Litt, entitled "One Class Novelty Detection for Seizure Analysis from Intracranial EEG," J. Machine Learning Research 7 1025-1044, (2006). Although a binary classifier is described, it is understood that the classifier may include a multi- class classifier. Accordingly, the kinetic features may be classified into at least two classes.

The example classifier contains two stages: (a) training and (b) testing. The features corresponding to each of the feature classes are used as inputs to the classifier individually and in combination. From the training data, a hyper-plane is created in the eigen-space that optimally separates the data into benign and malignant ROI classes. Given a set of labeled training data from two distinct classes, the example SVM classifier project the data into a high dimensional space constructed by a kernel function, ψ, operating on the training data. Testing data are then classified according to where they fall in relation to the hyper-plane when operated on by the same kernel function ψ. The objects of each class that lie closest to this hyper-plane are the "support vectors." The general form of the SVM classifier is given by equation (6):

V( ^χ ) = ∑ eτifrφ( ^χ , Xr) + b, τ=i (6) where x is the input training data, x _τ , τ e {1, 2, ^■ • ^■ ,N _S } denotes the support vectors, y e{— 1, 1} as the training labels, ψ(-, •) is a positive, definite, symmetric kernel function, b is a bias obtained from the training set to maximize the distance between the support vectors, and ξ is a model parameter chosen to maximize the objective function shown in equation (7) :

T=I τ,p-\ (7)

The kernel function, ψ(-, •), defines the nature of the decision hyper-plane. The example SVM uses a common kernel called the radial basis function (RBF). The parameters p and b are found through empirical training and testing of the classifier. After classifying the data in step 228 using the example SVM classifier, it may be desirable to reduce the dimensionality of the data in step 230 using a nonlinear dimensionality reduction technique such as graph embedding (GE) or locally linear

embedding (LLE). GE is described in an article by S. Yan et al. entitled "Graph Embedding : A General Framework for Dimensionality Reduction," Proc. 2005 Internal Conference on Computer Vision and Pattern Recognition. LLE is described in a publication by L. K. Saul et al. entitled "An Introduction to Locally Linear Embedding" which is available at http://www.cs.toronto.edu/~roweis/lle/papers/lleintro.pdf.

After dimensional reduction in step 230, the last step in the process is to display the result in step 232 using the display device of the display/keyboard 116, shown in Fig. 1. The output data may be displayed as a simple result, indicating, for example, a probability that a particular ROI is malignant and/or a triple negative lesion. Alternatively, the result may be displayed in reduced dimensions as a point relative to a reduced-dimension version of the hyper-plane that defines the boundary between malignant and non-malignant or TN and non-TN lesions. This second display may be advantageous as it may provide a readily discernable distance between the data point corresponding to the ROI and the reduced-dimension hyper-plane. The efficacy of the feature set is then evaluated based on the ability of the feature set to correctly classify each ROI, using the remaining ROIs in the dataset as the training set. Accuracy is defined as (tp + tn)/(tp+tn+fp+fn), where tp is the number of true positives, tn is the number of true negatives, fp is the number of false positives, and fn is the number of false negatives. In the example embodiment, after the system has been run to identify malignant ROIs, the data for the malignant ROIs is separated from the non-malignant ROIs and the classifier algorithm is run again on that data to separate TN lesions from non-TN malignant lesions.

Table 2 shows the ability of the example system to distinguish benign ROIs from malignant ROIs and Table 3 shows the ability of the example system to distinguish TN lesions from non-TN lesions. These tables quantify the accuracy (Ace), sensitivity (Sens.) and specificity (Spec.) of the example method using respectively different feature sets to distinguish benign and malignant ROIs, in Table 2, and TN and non-TN lesions, in Table 3. Table 2: SVM evaluation for Benign vs. Malignant ROIs

Table 3: SVM evaluation for TN vs. non-TN lesions

From these results, it is apparent that there is an advantage to using the kinetic first order textural features and kinetic second order statistical features to detect malignant ROIs and, in particular, TN lesions.

Although the invention is illustrated and described herein with reference to specific embodiments, the invention is not intended to be limited to the details shown. Rather, various modifications may be made in the details within the scope and range of equivalents of the claims and without departing from the invention.