# Automated Detection and Segmentation of Nonmass-Enhancing Breast Tumors with Dynamic Contrast-Enhanced Magnetic Resonance Imaging

Authors: Ignacio Alvarez Illan; Javier Ramirez; J. M. Gorriz; Maria Adele Marino; Daly Avendano; Thomas Helbich; Pascal Baltzer; Katja Pinker; Anke Meyer-Baese
Journal: Contrast Media & Molecular Imaging
 (2018)

---

## Abstract

Nonmass-enhancing (NME) lesions constitute a diagnostic challenge in dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) of the breast. Computer-aided diagnosis (CAD) systems provide physicians with advanced tools for analysis, assessment, and evaluation that have a significant impact on the diagnostic performance. Here, we propose a new approach to address the challenge of NME lesion detection and segmentation, taking advantage of independent component analysis (ICA) to extract data-driven dynamic lesion characterizations. A set of independent sources was obtained from the DCE-MRI dataset of breast cancer patients, and the dynamic behavior of the different tissues was described by multiple dynamic curves, together with a set of eigenimages describing the scores for each voxel. A new test image is projected onto the independent source space using the unmixing matrix, and each voxel is classified by a support vector machine (SVM) that has already been trained with manually delineated data. A solution to the high false-positive rate problem is proposed by controlling the SVM hyperplane location, outperforming previously published approaches.


---

## Body

## 1. Introduction

Accurate methods for early diagnosis of breast cancer are pivotal and contribute to an improved prognosis and survival outcomes in breast cancer patients. There is a consensus that dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is the most sensitive test for breast cancer detection and the backbone of any MRI protocol, enabling simultaneous assessment of tumor morphology and enhancement kinetics that evaluate neoangiogenesis as a tumor-specific feature. DCE-MRI has an excellent sensitivity and good specificity for lesions presenting as mass enhancement [1]. However, nonmass-enhancing (NME) lesions exhibit a heterogeneous appearance with high variations in kinetic characteristics and morphological patterns on DCE-MRI [2]. Consequently, DCE-MRI has reported lower specificity and sensitivity of 35% and 73% for NME lesions, much lower than those for mass-enhancing lesions. A set of computer-aided diagnosis (CAD) systems for breast cancer diagnosis on DCE-MRI has been developed with satisfactory performance results. However, in breast tumors presenting as NME lesions, the performance with low specificity is still suboptimal.For a CAD system to be used in breast DCE-MRI, two features are important to evaluate (i) the ability of the CAD to correctly differentiate between malignant and benign lesions and (ii) the ability of the CAD system to correctly locate malignant lesions within the 3D spatial volume. To evaluate the first feature, the diagnostic accuracy, specificity, and sensitivity are usually reported. To evaluate the second feature, commonly, the Dice similarity coefficient (DCS) is calculated between the CAD segmentation and some other ground truth segmentation. In most cases, a manual segmentation of the lesions is performed by experienced radiologists as ground truth. However, it has to be noted that even with expert’s interpretation, visual readings are prone to subjective errors [3], and specificity of DCE-MRI is limited particularly in small and nonmass-enhancing lesions, resulting in unnecessary breast biopsies [4]. In addition, CAD systems for breast cancer diagnosis have a reported high false-positive rate and, consequently, low specificity. However, this does not necessarily mean that CAD systems misclassify benign lesions as malignant. Therefore, it is not clear whether CAD systems can be optimized to improve lesion segmentation independently from lesion classification, or if irregardless, they will inherently suffer from the same limitations such as the low specificity reported in visual readings of DCE-MRI.In this work, we examined the relationship between the false-positive rate of CAD systems for breast cancer diagnosis and lesion segmentation on DCE-MRI. To achieve our aim, we obtained rich characterization of data through advanced processing techniques, combined with machine-learning paradigms intended for big data analysis, and used the resulting information to build a CAD system. We did not introduce any a priori knowledge about the disease in the workflow in order that all information may be completely data driven, which thereby also enabled us to identify new features not currently in the Breast Imaging-Reporting and Data System (BI-RADS) classification criteria that could potentially improve segmentation of visual readings. Both morphological and kinetic descriptors are considered in BI-RADS lexicons. However, in NME lesions, morphological descriptors are hard to define, and therefore, kinetic behavior can be an important source of information. Therefore, using only dynamic information of the tissue, we performed a supervised method to detect and segment nonmass-enhanced lesions on the breast.Lesion segmentation has been successfully achieved using unsupervised clustering methods [5], fuzzy c-means (FCM) [6], or improvements over FCM [7]. In unsupervised clustering, sophisticated preprocessing must be implemented to control the false-positive rate, with fine tuning of parameters and/or heuristic steps. On the contrary, it has been demonstrated that processing of dynamic signals provides relevant information for classification of tissues, such as principal component analysis- (PCA-) based decompositions closely related to the 3TP method [8].Thus, we undertook a combination of supervised segmentation and signal processing to successfully segment NME lesions with control of the false-positive rate. Independent component analysis (ICA) was used to extract a set of independent curves that described the possible dynamic behavior of different breast tissues. ICA has been shown to provide richer descriptions of underlying patterns than PCA [9, 10], and therefore, it was used for supervised classification in our work. We also incorporated machine learning, whereby we trained a classifier using the information encoded in a whole dataset of subjects, including the dynamic behavior of benign and malignant tissues. Considering features at the voxel level, the system “learned” to characterize malignant tissues with a support vector machine (SVM). A procedure was implemented to fix the SVM hyperplane location, reducing and controlling the false-positive rate. Projecting new unseen data using the unmixing matrix allowed us to obtain the features for estimating the generalization capabilities in a cross-validation scheme and compare them with visual readings of the images reported in the literature and other CAD system approaches.The methods proposed within this work demonstrate that NME lesions can be detected with kinetic information by using multiple enhancement curves, providing a promising approach for improving breast cancer diagnosis. Accurate diagnostic methods as the one we hereby present may have an impact not only on accurate diagnosis but also in reducing unnecessary breast biopsies.


### 1.1. Related Work

The use of CAD systems to improve visual readings of DCE-MRI in breast cancer ranges from purely visual methods to automatic classification. The present work combines visual comparison aspects with automatic classification techniques, thus adding a value to purely visual comparison techniques based on PCA or self-organizing map (SOM), such as in [8, 11], and complementing pure classification approaches, such as in [12, 13]. Specifically, the PCA approach in [8] extends the three-point technique (3PT) by adding an eigenvector decomposition of the time signals. However, that decomposition does not provide an independent set of sources, but only a set of uncorrelated ones. The time-intensity curve estimation in [14] also seeks for hidden kinetics, but applies them to mass lesions. Concerning the automatic classification CADs, most approaches are concentrated on the detection and classification of mass-enhancing lesions, by combining kinetic and morphological features [12, 13, 15, 16], like shape, margins, and internal enhancement distribution [17], textural kinetics [18], or more recently using deep neural networks [19, 20], among others. The detection and segmentation of lesions are usually performed as a manual or semimanual task, in which regions of interest (ROIs) are manually defined or obtained from seeds with manual inputs.For automatic lesion segmentation, keeping an acceptable false-positive rate is a common issue in DCE-MRI CAD systems of the breast [21]. In many of these cases, unsupervised methods for lesion segmentation, such as FCM algorithms in [6, 22], are used, and then the features extracted from the lesions are used for classification. Complex workflows that include vessel detection, whole-breast segmentation, and several preprocessing steps have been proposed to control false-positive detection [5, 7, 23, 24].


## 1.1. Related Work

The use of CAD systems to improve visual readings of DCE-MRI in breast cancer ranges from purely visual methods to automatic classification. The present work combines visual comparison aspects with automatic classification techniques, thus adding a value to purely visual comparison techniques based on PCA or self-organizing map (SOM), such as in [8, 11], and complementing pure classification approaches, such as in [12, 13]. Specifically, the PCA approach in [8] extends the three-point technique (3PT) by adding an eigenvector decomposition of the time signals. However, that decomposition does not provide an independent set of sources, but only a set of uncorrelated ones. The time-intensity curve estimation in [14] also seeks for hidden kinetics, but applies them to mass lesions. Concerning the automatic classification CADs, most approaches are concentrated on the detection and classification of mass-enhancing lesions, by combining kinetic and morphological features [12, 13, 15, 16], like shape, margins, and internal enhancement distribution [17], textural kinetics [18], or more recently using deep neural networks [19, 20], among others. The detection and segmentation of lesions are usually performed as a manual or semimanual task, in which regions of interest (ROIs) are manually defined or obtained from seeds with manual inputs.For automatic lesion segmentation, keeping an acceptable false-positive rate is a common issue in DCE-MRI CAD systems of the breast [21]. In many of these cases, unsupervised methods for lesion segmentation, such as FCM algorithms in [6, 22], are used, and then the features extracted from the lesions are used for classification. Complex workflows that include vessel detection, whole-breast segmentation, and several preprocessing steps have been proposed to control false-positive detection [5, 7, 23, 24].


## 2. Methods

Each voxel of the DCE-MRI image has a time signal representing the enhancement kinetics of the different contributing breast tissues. A set of DCE-MRI time signals can be analyzed in terms of the blind source separation problem, which proposes that the different dynamic behaviors can be expressed as a linear combination of a reduced set of sources, making very little assumptions on the nature of that combination. Those sources and their scores can be used as features for classification, as depicted in Figure1.Figure 1
(a) Time sequence of database images. (b) Decomposition of a sample time signalx into a linear combination of independent sources by ICA and its corresponding scores a1, a2, …, as. (c) Scatter plot of the first scores and the SVM hyperplane classifier.


### 2.1. ICA-Based Enhancement Curve Analysis

ICA offers a solution to the blind source separation problem estimating a set of sources that maximize the statistical independence between them, measured in terms of a cost function. In the literature, several functions have been used to measure statistical independence between signals [25]. Here, we used the FastICA algorithm [26] with mutual information as a measure function. Contrary to other eigenimage decompositions based on spatial ICA, like in face recognition [27] and brain imaging [9, 28], the independent sources are obtained here in the temporal domain; in other words, we work on a voxel level.Thus, each voxel defines a temporal curvextj with t1,…,tN temporal points. A set of voxels xi,i=1,…,M forms an image and defines the N×M matrix X of observed signals. The ICA task is to find the mixing matrix A and the set of sources S:(1)X=AS.The mixing matrixA is an N×N matrix that linearly combines the independent “images.” Contrary to other related methods, such as PCA, ICA does not provide a natural way to sort the N independent components. However, it is a relevant question whether or not a reduced set p<N of components contains noisy and discardable information. The mean squared error (MSE) between the enhancement time signals and the reconstructed signals using the k source sk is calculated as follows:(2)MSEk=1Nt⋅Nr∑i,jxitj−ajk⋅ski2,and used as a parameter to measure the noise content of each sk source, with k=1,…,N.When working at the voxel level, Equation (1) can also be understood as a linear decomposition of each vector xtj into a set of temporal sources whose coefficients belong to the independent sources. Therefore, each voxel location xi has N coefficients sji to j=1,…,N, whose values are maximally independent and measure the importance of each temporal source to recover that voxel dynamics, by linearly combining them (Figure 1). In the rest of the paper, we will refer to these coefficients as the scores.It is important to stress that working on a voxel level will allow data from different patients to be included in the matrix setX. Therefore, the obtained set of sources S does not have to be restricted to represent the particular dynamic enhancement present in a single subject but can be used to model all the possible curves that independently characterize each BI-RADS category.For new unseen datax˜ at the voxel level, the scores are extracted from x˜ by projecting it onto the subspace E spanned by the signals from the matrix A. Specifically, let a1,…,ap be the basis set of temporal curves spanning the subspace E and then A denote the N-by-p matrix of which columns are a1,…,ap. Let p≤N, as some of the signals may have been removed due to their noisy nature. Since this basis need not be orthogonal, a well-known result of linear algebra stated that the projection is given by(3)PA=AATA−1AT,so that the application of that operator on a voxel signal x˜t:(4)s=PAx˜,projects it to the subspace E, obtaining its p scores s on that subspace.The independent component scoressk of the dataset are used as feature vector inputs of a SVM to learn the different enhancement patterns associated with malignant and benign tissues.


### 2.2. False-Positive Rate Control by SVM Hyperplane Translation

SVM is a machine-learning algorithm that separates a given set of binary labeled training data with a hyperplane that is maximally distant from the two classes (known as the maximal margin hyperplane). The objective is to build a functionf:IRp⟶1,0 using training data, consisting of p dimensional patterns xi and class labels yi:(5)x1,y1,x2,y2,…,xM,yM∈IRp×1,0,so that f will correctly classify new examples x˜,y. The problem of finding the maximal margin hyperplane is usually solved by quadratic programming algorithms that try to minimize a margin cost function J:(6)Jw,w0,ξ=12w2+C∑i=1lξi,subject to the inequality constraints:(7)yiwTxi+w0≥1−ξi,ξi≥0i=1,2,…,l,where the slack variables ξi incorporate to the optimization of those feature vectors that are not separable (details can be found in [29]). The solution to that problem can be expressed by a linear combination of a subset of vectors, called support vectors:(8)dx=∑i=1NSαiyiKsi,x+w0,where K.,. is the kernel function, αi is a weight constant derived from the SVM process, and si are the NS support vectors [29]. Taking the sign of the function leads to the binary classification solution.Here, we propose an SVM hyperplane translation in terms of the slack variablesξi to control the number of false positives. We add a new term gsi,ξi to the hyperplane-defining function dx so that the classification solution is now defined by(9)fx=signdx+gsi,ξi,where the function g takes the two-class average distance to the hyperplane of those support vectors with ξi>1, measured by the kernel metric K. Common kernels that are used by SVM practitioners for the nonlinear feature mapping are as follows:(i)
Polynomial function:(10)Kx,y=γx⋅y+cd.(ii)
Radial basis function (RBF):(11)Kx,y=exp−γx−y2,as well as the linear kernel, in which K.,. is simplified as a scalar product, and therefore, g in Equation (9) would average the Euclidean distance in that particular case.


### 2.3. Dataset

The dataset used for analysis consisted of sixteen patients that presented with NME breast tumors at DCE-MRI. This patient cohort is a subset from a larger cohort undergoing multiparametric MRI using inclusion criteria described in detail in [30]. All patients underwent MRI of the breast using a 3T MRI scanner (Tim Trio; Siemens, Erlangen, Germany) with a dedicated, bilateral, 4-channel breast coil in vivo (Orlando, FL), and the imaging protocol comprised both high-spatial and -temporal resolution. Three high-spatial resolution images were taken, precontrast, peak, and postcontrast as a coronal T1-weighted (3D) FLASH sequence, with water excitation and fat suppression, with the following sequence parameters: TR/TE 877/3.82 milliseconds, FOVr 320 mm, SI 1 mm isotropic, 96 slices, flip angle 9°, matrix 320/134, 1 average, and acquisition time 2 minutes. A high-temporal resolution, contrast-enhanced, coronal T1-weighted (VIBE) sequence was obtained with the following sequence parameters: TR/TE 3.61/1.4 milliseconds, FOVr 320 mm, SI 1.7 mm isotropic, 72 slices, flip angle 6°, matrix 192/192, 1 average, and 13.2 seconds of acquisition time per volume leading to 3.45 minutes for 17 measurements. A second set of high-spatial resolution T1-weighted imaging (repeated 3D-FLASH) was acquired after these 17 low-spatial VIBE resolution images, as the peak enhancement of the lesion could be expected at the end of this time span ([30] and references therein). Finally, high-temporal resolution (repeated VIBE with 25 measurements, leading to an acquisition time of 5 minutes 35 seconds, and repeated 3D-FLASH for dynamic assessment of lesion wash-out) was performed, and then high-spatial resolution T1-weighted images were recorded. The contrast agent used was Gd-DOTA (generic name: gadoterate meglumine; Dotarem, Guerbet, France), injected intravenously as a bolus (0.1 mmol per kilogram body weight) and administered with a power injector (Spectris Solaris EP; Medrad, Pittsburgh, PA) at 4 mL/s followed by a 20 mL saline flush. The contrast agent was injected 75 seconds after starting the first coronal T1-weighted VIBE.NME breast tumors were visually assessed by three expert radiologists following the American College of Radiology BI-RADS Atlas [31] and delineated using the OsiriX software on the 3T high-spatial resolution volumes. All NME lesions were classified as BI-RADS 4: suspicious, or BI-RADS 5: highly suspicious of malignancy. Histopathology was used as the standard of reference. There were eleven invasive ductal carcinomas (IDCs), three ductal carcinomas in situ (DCISs), and two invasive lobular carcinomas (ILCs).


### 2.4. Preprocessing

All dynamic sequences were registered to the precontrast volume. This preprocessing step was required to remove any spatial misalignments on the sequence caused by involuntary movements of the patient. The algorithm employed to perform this task was the SPM12 [32] registration algorithm, which performs affine and nonaffine transformations on the data by minimizing a similarity measure cost function, selected to be the mutual information metric. Afterwards, a 3D Gaussian filter of size 2FWHM was used to smoothen the images.In spite of the existence of automatic and accurate methods for performing whole-breast segmentation [33, 34, 35, 36], we performed this task straightforwardly finding the middle chest point as in [12], and discarding the content of the image after this point, reducing the original number of 192×192×72≈2.6⋅106 voxels contained in each image to ≈1.6⋅105, and guaranteeing the exclusion of heart and other organs’ noisy signals. Concretely, the middle chest point was obtained by performing the following steps (Figure 2):(1)
Compute the cross-correlation of the convolution of the image with itself in the sagittal direction. The middle sagittal plane will lie in the symmetry plane of the body, and due to its symmetry, it will reach the maximum convoluted cross-correlation.(2)
Compute the intensity gradient of the middle chest slice in the coronal direction and find its maximummy. Remove the internal part of the image that lies in the coronal direction after the middle chest plane y=my.Figure 2
Exclusion of internal organs by detection of the middle chest plane. (a) Middle sagittal plane. (b) Middle chest plane.
(a)
(b)The described procedure ensured the removal of voxels that lie inside the thoracic cavity and the chest wall as well as background voxels.The manual delineations of the lesions were performed by three expert radiologists on 3T high-spatial resolution images using the OsiriX software, recorded as a set of axial point coordinates in mm. The Bresenham algorithm [37] was used to transform the coordinate points into 3D binary masks, and a decimation was employed to downsample the masks to the size of the low-spatial high-temporal resolution images. Thus, the downsampled masks were used to define the class labels of each voxel: 1 if the voxel was in the mask, and 0 otherwise.


## 2.1. ICA-Based Enhancement Curve Analysis

ICA offers a solution to the blind source separation problem estimating a set of sources that maximize the statistical independence between them, measured in terms of a cost function. In the literature, several functions have been used to measure statistical independence between signals [25]. Here, we used the FastICA algorithm [26] with mutual information as a measure function. Contrary to other eigenimage decompositions based on spatial ICA, like in face recognition [27] and brain imaging [9, 28], the independent sources are obtained here in the temporal domain; in other words, we work on a voxel level.Thus, each voxel defines a temporal curvextj with t1,…,tN temporal points. A set of voxels xi,i=1,…,M forms an image and defines the N×M matrix X of observed signals. The ICA task is to find the mixing matrix A and the set of sources S:(1)X=AS.The mixing matrixA is an N×N matrix that linearly combines the independent “images.” Contrary to other related methods, such as PCA, ICA does not provide a natural way to sort the N independent components. However, it is a relevant question whether or not a reduced set p<N of components contains noisy and discardable information. The mean squared error (MSE) between the enhancement time signals and the reconstructed signals using the k source sk is calculated as follows:(2)MSEk=1Nt⋅Nr∑i,jxitj−ajk⋅ski2,and used as a parameter to measure the noise content of each sk source, with k=1,…,N.When working at the voxel level, Equation (1) can also be understood as a linear decomposition of each vector xtj into a set of temporal sources whose coefficients belong to the independent sources. Therefore, each voxel location xi has N coefficients sji to j=1,…,N, whose values are maximally independent and measure the importance of each temporal source to recover that voxel dynamics, by linearly combining them (Figure 1). In the rest of the paper, we will refer to these coefficients as the scores.It is important to stress that working on a voxel level will allow data from different patients to be included in the matrix setX. Therefore, the obtained set of sources S does not have to be restricted to represent the particular dynamic enhancement present in a single subject but can be used to model all the possible curves that independently characterize each BI-RADS category.For new unseen datax˜ at the voxel level, the scores are extracted from x˜ by projecting it onto the subspace E spanned by the signals from the matrix A. Specifically, let a1,…,ap be the basis set of temporal curves spanning the subspace E and then A denote the N-by-p matrix of which columns are a1,…,ap. Let p≤N, as some of the signals may have been removed due to their noisy nature. Since this basis need not be orthogonal, a well-known result of linear algebra stated that the projection is given by(3)PA=AATA−1AT,so that the application of that operator on a voxel signal x˜t:(4)s=PAx˜,projects it to the subspace E, obtaining its p scores s on that subspace.The independent component scoressk of the dataset are used as feature vector inputs of a SVM to learn the different enhancement patterns associated with malignant and benign tissues.


## 2.2. False-Positive Rate Control by SVM Hyperplane Translation

SVM is a machine-learning algorithm that separates a given set of binary labeled training data with a hyperplane that is maximally distant from the two classes (known as the maximal margin hyperplane). The objective is to build a functionf:IRp⟶1,0 using training data, consisting of p dimensional patterns xi and class labels yi:(5)x1,y1,x2,y2,…,xM,yM∈IRp×1,0,so that f will correctly classify new examples x˜,y. The problem of finding the maximal margin hyperplane is usually solved by quadratic programming algorithms that try to minimize a margin cost function J:(6)Jw,w0,ξ=12w2+C∑i=1lξi,subject to the inequality constraints:(7)yiwTxi+w0≥1−ξi,ξi≥0i=1,2,…,l,where the slack variables ξi incorporate to the optimization of those feature vectors that are not separable (details can be found in [29]). The solution to that problem can be expressed by a linear combination of a subset of vectors, called support vectors:(8)dx=∑i=1NSαiyiKsi,x+w0,where K.,. is the kernel function, αi is a weight constant derived from the SVM process, and si are the NS support vectors [29]. Taking the sign of the function leads to the binary classification solution.Here, we propose an SVM hyperplane translation in terms of the slack variablesξi to control the number of false positives. We add a new term gsi,ξi to the hyperplane-defining function dx so that the classification solution is now defined by(9)fx=signdx+gsi,ξi,where the function g takes the two-class average distance to the hyperplane of those support vectors with ξi>1, measured by the kernel metric K. Common kernels that are used by SVM practitioners for the nonlinear feature mapping are as follows:(i)
Polynomial function:(10)Kx,y=γx⋅y+cd.(ii)
Radial basis function (RBF):(11)Kx,y=exp−γx−y2,as well as the linear kernel, in which K.,. is simplified as a scalar product, and therefore, g in Equation (9) would average the Euclidean distance in that particular case.


## 2.3. Dataset

The dataset used for analysis consisted of sixteen patients that presented with NME breast tumors at DCE-MRI. This patient cohort is a subset from a larger cohort undergoing multiparametric MRI using inclusion criteria described in detail in [30]. All patients underwent MRI of the breast using a 3T MRI scanner (Tim Trio; Siemens, Erlangen, Germany) with a dedicated, bilateral, 4-channel breast coil in vivo (Orlando, FL), and the imaging protocol comprised both high-spatial and -temporal resolution. Three high-spatial resolution images were taken, precontrast, peak, and postcontrast as a coronal T1-weighted (3D) FLASH sequence, with water excitation and fat suppression, with the following sequence parameters: TR/TE 877/3.82 milliseconds, FOVr 320 mm, SI 1 mm isotropic, 96 slices, flip angle 9°, matrix 320/134, 1 average, and acquisition time 2 minutes. A high-temporal resolution, contrast-enhanced, coronal T1-weighted (VIBE) sequence was obtained with the following sequence parameters: TR/TE 3.61/1.4 milliseconds, FOVr 320 mm, SI 1.7 mm isotropic, 72 slices, flip angle 6°, matrix 192/192, 1 average, and 13.2 seconds of acquisition time per volume leading to 3.45 minutes for 17 measurements. A second set of high-spatial resolution T1-weighted imaging (repeated 3D-FLASH) was acquired after these 17 low-spatial VIBE resolution images, as the peak enhancement of the lesion could be expected at the end of this time span ([30] and references therein). Finally, high-temporal resolution (repeated VIBE with 25 measurements, leading to an acquisition time of 5 minutes 35 seconds, and repeated 3D-FLASH for dynamic assessment of lesion wash-out) was performed, and then high-spatial resolution T1-weighted images were recorded. The contrast agent used was Gd-DOTA (generic name: gadoterate meglumine; Dotarem, Guerbet, France), injected intravenously as a bolus (0.1 mmol per kilogram body weight) and administered with a power injector (Spectris Solaris EP; Medrad, Pittsburgh, PA) at 4 mL/s followed by a 20 mL saline flush. The contrast agent was injected 75 seconds after starting the first coronal T1-weighted VIBE.NME breast tumors were visually assessed by three expert radiologists following the American College of Radiology BI-RADS Atlas [31] and delineated using the OsiriX software on the 3T high-spatial resolution volumes. All NME lesions were classified as BI-RADS 4: suspicious, or BI-RADS 5: highly suspicious of malignancy. Histopathology was used as the standard of reference. There were eleven invasive ductal carcinomas (IDCs), three ductal carcinomas in situ (DCISs), and two invasive lobular carcinomas (ILCs).


## 2.4. Preprocessing

All dynamic sequences were registered to the precontrast volume. This preprocessing step was required to remove any spatial misalignments on the sequence caused by involuntary movements of the patient. The algorithm employed to perform this task was the SPM12 [32] registration algorithm, which performs affine and nonaffine transformations on the data by minimizing a similarity measure cost function, selected to be the mutual information metric. Afterwards, a 3D Gaussian filter of size 2FWHM was used to smoothen the images.In spite of the existence of automatic and accurate methods for performing whole-breast segmentation [33, 34, 35, 36], we performed this task straightforwardly finding the middle chest point as in [12], and discarding the content of the image after this point, reducing the original number of 192×192×72≈2.6⋅106 voxels contained in each image to ≈1.6⋅105, and guaranteeing the exclusion of heart and other organs’ noisy signals. Concretely, the middle chest point was obtained by performing the following steps (Figure 2):(1)
Compute the cross-correlation of the convolution of the image with itself in the sagittal direction. The middle sagittal plane will lie in the symmetry plane of the body, and due to its symmetry, it will reach the maximum convoluted cross-correlation.(2)
Compute the intensity gradient of the middle chest slice in the coronal direction and find its maximummy. Remove the internal part of the image that lies in the coronal direction after the middle chest plane y=my.Figure 2
Exclusion of internal organs by detection of the middle chest plane. (a) Middle sagittal plane. (b) Middle chest plane.
(a)
(b)The described procedure ensured the removal of voxels that lie inside the thoracic cavity and the chest wall as well as background voxels.The manual delineations of the lesions were performed by three expert radiologists on 3T high-spatial resolution images using the OsiriX software, recorded as a set of axial point coordinates in mm. The Bresenham algorithm [37] was used to transform the coordinate points into 3D binary masks, and a decimation was employed to downsample the masks to the size of the low-spatial high-temporal resolution images. Thus, the downsampled masks were used to define the class labels of each voxel: 1 if the voxel was in the mask, and 0 otherwise.


## 3. Experiments

The dataset was divided into three subsets: training data, validation data, and testing. Training and validation data comprised half of the dataset, while the test set consisted of the other half. The data were considered at the voxel level. Therefore, after discarding nonrelevant parts of the image, a random selection ofNa≈5⋅103 benign voxel samples from the pool of all nonlesion voxels of the images was performed to balance the training set, resulting in a 2∗Na×p training and validating data matrix.The voxel data were used as input to the FastICA algorithm, obtaining a set of scores for each voxel that served as feature vectors for training and validating an SVM in a cross-validation scheme. The validation step is performed in two stages:(i)
Firstly, different parameters were optimized within a 10-fold cross-validation scheme: (i) the optimal dimensionality of the datah and (ii) the optimal kernel (linear, polynomial, or RFB). The optimal value for h was obtained by sorting the independent components by their MSE defined in Equation (2), and the feature space dimension was changed by sequentially increasing the number of components included on the scores. The optimal kernel was selected by comparing the classification performance, based on the classification error.(ii)
Secondly, once the number of components and the kernel function were fixed, the decision boundary location of the SVM was analyzed in an enlarged test dataset of size≈4⋅105, that contained all the discarded voxels in the validation step.


## 4. Results

The scores defined in Equation (1) are depicted in two different spaces: the 3DDCE-MRI space coregistered with the original data (Figure 3) and the E subspace spanned by the first two temporal sources a1 and a2 (Figure 4), sorted according to the MSE-defined criteria. The representation in the 3DDCE-MRI space shows that similar score values are grouped together around tissues that have a similar enhancement. On the bottom, voxels belonging to the lesions present a high score value, revealing that the associated independent component encodes the malignant dynamic information. On the top, the distribution of score values does not concentrate on specific regions but spreads over the breast tissues revealing a relation with normal tissue enhancement dynamics. That information complements the representation on the E subspace, where a clear separation between tumor tissues represented in blue and normal tissues in red can be inferred, although some regions of overlapping are present. Also, the independent components s1 and s2 are shown in Figure 5, together with other extracted sources. It is interesting to note that being automatically data-driven extracted, these independent components take the form of enhancement curves: curve IC1 is a normal enhancement, while curve IC2 has a “typical” malignant behavior, according to model-based descriptions [1]. The remaining set of independent components cannot be assigned to any particular dynamic nor tends to form clusters of similar enhancement when depicted in 3D, therefore not possessing an obvious interpretation. However, the common classification into wash-out, plateau, and permanent enhancement of dynamic curves is reduced by ICA to only two clearly identifiable curves. Therefore, the ICA-based signal processing analysis reveals that dynamic enhancement curves reaching a plateau do not behave independently in the ICA sense from wash-out curves, while permanent enhancement curves do.Figure 3
One middle axial slice for the 16 patients from the projected independent component space of the dataset. Intensity represents voxelwise scores of the first and second independent components in the 3D MRI space. The scores of the first IC, on the top, correspond with the IC1 in blue in Figure5. The scores of the second IC, on the bottom, correspond with the IC2 in red in Figure 5. High values on the IC2 (in yellow) can be related to malignancy.Figure 4
Scatter plot of the scores corresponding to the first two independent components of the training data, together with the linear decision SVM function (in black) and the support vectors (SVs).Figure 5
First four independent components sorted by MSE. In red, the IC2 shows typical “malignant” dynamics, while in blue, IC1 shows a persistent enhancement curve, characteristic of benign tissues.The results of the cross-validation are shown in Figures4 and 6 and also in the left part of Table 1. In Figure 4, the 2∗Na×p training data are shown after the SVM is trained, and the obtained support vectors are marked with circles. From Figure 6, the optimal number of components used to reconstruct the signal is above 5, revealing that a simple decomposition of signals into benign and malignant behaviors can be enriched with other significant components reaching ROC values over 0.90.Figure 6
ROC and area under the curve (AUC) values on the cross-validation scheme for the RBF kernel by varying the number of components on ICA. ICA components are sorted according to MSE.Table 1
Performance parameters on training and validation data.
Training
Validation DSC
Hinge loss
Accuracy
Specificity
Sensitivity
d
=
0
Max.
DSC(μd) [DSC(μd±σd)]
PCA + linear SVM
0.9764
0.7263
0.6581
0.7944
0.31±0.01
0.3382±0.0005
0.3310 [0.3039 − 0.2169]
PCA + RBFSVM
0.9529
0.7263
0.6581
0.7944
0.31±0.01
0.3382±0.0005
0.3310 [0.3039 − 0.2169]
ICA + linear SVM
0.1254
0.9501
0.9410
0.9593
0.31±0.01
0.53
±
0.01
0.5295 [0.475 − 0.484]
ICA + RBFSVM
0.1083
0.9515
0.9573
0.9457
0.29±0.01
0.44±0.04
0.1085 [0.3711 − 0.0559]
Raw + linear SVM
2.4429
0.8026
0.8446
0.7605
0.15±0.07
0.30±0.05
0.2325 [0.1373 − 0.3058]Figure7 shows the NME lesion delineated by the expert radiologist (in red), together with a distance-to-hyperplane map (distance d=0 is represented by a black contour). The value of each voxel in the map is defined in Equation (8). It can be seen that hyperplane location (value d=0) produces big regions of false positives. Those regions are mostly concentrated around the delimited lesion, but extended regions can also be found in nonconnected regions where benign dynamics are expected. The false-positive rate can be controlled by modifying the defining value of the hyperplane location, set to 0 by definition in SVM. Translating the hyperplane towards the positive values produces a more conservative definition of feature vectors belonging to the +1 class. Therefore, only score values high above the hyperplane would be considered as malignant, while intermediate values not clearly projecting malignant-related score values will not be classified as lesion, decreasing the false-positive rate and increasing specificity. However, there must be a compromise between specificity and sensitivity, since increasing the defining value of the decision function also has an impact on the false-negative rate.This trade‐off requires to be very finely tuned, as the number of benign samples is several orders of magnitude bigger than the number of malignant samples, producing an imbalanced classification problem. In Figure 8, the influence of the imbalanced classes can be perceived if compared with the scatter plot of the scores considering only the reduced training data of Figure 4. Although other solutions exist to the problem of imbalanced dataset in SVM classification, we propose here a very conservative approach, in which the hyperplane-defining value is translated into the +1 class region, guaranteeing that only very distant scores from the hyperplane are considered as malignant.The hyperplane-defining value g is given in Equation (9). Other values could be used to make this transformation but are prone to be affected by outlier support vectors that uncontrollably increase the false-negative rate. By averaging the support vector’s distance to the hyperplane with the condition ξi>1, we are smoothing the effect of possible outlier support vectors, while translating the hyperplane to actual relevant values. Alternatively, we calculate the decision-defining value experimentally, in the second validation on the training data, and test both on the test set: the theoretically derived value and the experimentally adjusted one. In the special case in which all ξi are less than 1, we average the support vector’s distance to the hyperplane with the condition 1>ξi>0.Figure 7
Five representative axial slices of an NME. The values on voxels represent the distance to the hyperplane after classification on a trained SVM. The black contour represents the location of the hyperplane atd=0, and the red contour is the manual delineation of the lesion.Figure 8
Scatter plot of the scores corresponding to the first two independent components of the validation data.To evaluate the lesion detection performance, the DSC is calculated as follows:(12)DSC=2∗A∩MA∪M,and the amount of overlap between segmentation algorithms (A) and manually generated (M) segmentations is measured with respect to the size of the segmented region.Table1 shows the validation values obtained by default SVM at d=0, at empirical maximum, and at the proposed value, for 2-component PCA, ICA, and raw data using 2 kernels. Raw data are displayed for reference and correspond to the use of dynamic curves as feature vectors for SVM, without multicurve extraction. The PCA method [8] shows higher DSC at d=0 than the proposed ICA approach. Hyperplane translation has a lower effect in the PCA case since all support vectors lie in the condition ξi<1. In the ICA with a linear kernel case, the false positives are reduced significantly reaching the maximum DSC values, in agreement with the interval of maximum empirical values.Figure9 reports a free-response receiver-operating characteristic (FROC) curve analysis [38] at the voxel level. Although in mass lesions FROC analysis is usually reported at the lesion level, in NME lesions FROC analysis at the lesion level can be misleading, as can be seen from Figure 7: increasing the confidence threshold increases the number of false-positive lesions due to lesion fragmentation, although false positives at the voxel level decrease. Two reference methods are shown for comparison: the signal enhancement ratio (SER) method, based on the following SER=(SI(t=1st postcontrast time point)–SI(t=precontrast time point))/(SI(t=final postcontrast time point)–SI(t=precontrast time point)), with a varying threshold; and the derivative SER, a modified version of the method that uses the Laplacian of the image to obtain the SER, as defined in the work of Levman et al. [21]. The FROC curve for the ICA-SVM method proposed in this paper is obtained on the test set by adding a varying threshold k to the SVM output in Equation (8) and computing the sign signdx+k.Figure 9
FROC curves for the proposed algorithm (SVM-ICA) in comparison with the references (SER and SER derivative [21]).


## 5. Discussion

The contributions of this work are twofold: first, visual interpretations of the DCE-MRI image can be enriched by using the proposed ICA-based processing of time signals, which produces a data-driven decomposition of dynamic enhancement signals into multicurve description signals, that are statistically independent and disease specific. The idea of producing multiple curves to characterize lesions has also been explored by Liu et al. [14], but from the total variation perspective, it is not data-driven but based on assumptions on the data. Other visual methods based on CAD techniques, such as PCA in Eyal et al. [8]or PCA-SOM-LDD in Varini et al. [11], have been proposed in the literature to enrich the well-known 3TP method. Thus, visual support is an important characteristic to evaluate in aiding diagnosis of breast cancer by computer systems. It is also important to stress that the ICA extraction must be done only in the training phase of the algorithm. The CAD system will then benefit from an online response, once the CAD is conveniently trained. The presented approach outperforms PCA-based methods as shown in Table 1 in terms of automatic segmentation performance and provides a meaningful visual support for experienced and unexperienced readers.The low incidence of NME lesions reduces the available testing data, therefore limiting the validation of the presented method. Moreover, the heterogeneous nature of NME lesions also limits the accuracy in lesion annotation performed by experts when compared to CAD segmentations. Therefore, the reported DSC values when comparing ground truth and CAD results must be understood as a lower bound estimation of the segmentation capabilities of the presented CAD, since a semiautomatic annotation can potentially boost the DSC values.The second contribution is the supervised nature of the detection and segmentation method, which allows control of the false-positive rate. Most CAD systems for lesion classification start from a manual or semimanual ROI delineation [8, 12, 16], that limits control of the false positives. The baseline approach to lesion segmentation is the FCM unsupervised method, which in Liang et al. [39]is reported to have a 6%±9% of overlap with manually defined ROIs, and is commonly used in many CAD systems for breast cancer diagnosis in DCE-MRI. In Jayender et al. [7], an enhancing preprocessing step is added to the usual FCM algorithm using linear dynamic system modeling. The overlap of the algorithm output with the radiologists’ segmentation and CAD stream output, computed in terms of DSC, was 0.77 and 0.72, respectively. In the unsupervised approach of Cui et al. [23], a combination of Gaussian mixture modeling and marker-controlled watershed transform was used to segment the lesions. The overall overlap ratio between the two radiologists’ manual segmentations and the proposed algorithm was 64.3%±10.4%. The supervised method of Liang et al. [39]shows overlap rates with the ground truth of 51%±26% and 48%±25%. This method required a robust intensity normalization method to make intrapatient comparisons, while the ICA method presented here characterizes the form of the curve, thus not requiring intensity normalization. Moreover, we report higher or comparable DSC values than those in the literature, even in the more challenging case of NME breast lesions. We also report better control of false-positive rate than the method proposed by Levman et al. [21], with sensitivity greater than 75% at 105 false-positive voxels. Derivative SER reaches sensitivity 40% at that level, outperforming SER as already proved.


## 6. Conclusions

This paper presents promising results for challenging NME breast lesion detection in DCE-MRI. We propose an approach that develops a linear expansion of features for every voxel in the image based on ICA, allowing for a multicurve characterization of the enhancement behavior, in contrast with usual single-curve voxel characterization. The data-driven obtained features are used to train and test an SVM with satisfactory performance. In addition, previously, the imbalanced nature of the interest class features limited automatic detection by supervised methods such as SVM. In this work, we propose parameter optimization on the SVM hyperplane location, such that the false-positive rate is controlled, thus providing a solution to the low specificity problem in CAD of breast cancer. With that optimization, the DSC value is increased approximately a 50% from the defaultd=0 margin value, reaching a peak value of 0.5295.



---

*Source: 5308517-2018-10-24.xml*