Computer-aided diagnosis for diagnostically challenging breast lesions in DCE-MRI based on image registration and integration of morphologic and dynamic characteristics

Diagnostically challenging lesions comprise both foci (small lesions) and non-mass-like enhancing lesions and pose a challenge to current computer-aided diagnosis systems. Motion-based artifacts lead in dynamic contrast-enhanced breast magnetic resonance to diagnostic misinterpretation; therefore, motion compensation represents an important prerequisite to automatic lesion detection and diagnosis. In addition, the extraction of pertinent kinetic and morphologic features as lesion descriptors is an equally important task. In the present paper, we evaluate the performance of a computer-aided diagnosis system consisting of motion correction, lesion segmentation, and feature extraction and classification. We develop a new feature extractor, the radial Krawtchouk moment, which guarantees rotation invariance. Many novel feature extraction techniques are proposed and tested in conjunction with lesion detection. Our simulation results have shown that motion compensation combined with Minkowski functionals and Bayesian classifier can improve lesion detection and classification.


Introduction
Contrast-enhanced magnetic resonance 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] which can reach diagnostic sensitivities up to 97% and specificities up to 76.5%.
Recent clinical research has shown that ductal carcinoma in situ (DCIS) with small invasive carcinoma can http://asp.eurasipjournals.com/content/2013/1/157 Figure 1 visualizes the flow diagram of the proposed computer-assisted system, including image registration, lesion segmentation, feature extraction, and evaluation.
Automatic motion correction represents an important prerequisite to a correct automated small lesion evaluation [15]. Motion artifacts are caused either by the relaxation of the pectoral muscle or involuntary patient motion and invalidate the assumption of the same spatial location within the breast of the corresponding voxels in the acquired volumes for assessing lesion enhancement. Especially for small lesions, the assumption of correct spatial alignment often leads to misinterpretation of the diagnostic significance of enhancing lesions [16].
Visual assessment of morphologic properties is a highly inter-observer variable [17], while automated computation of features leads to more reproducible indices and thus to a more standardized and objective diagnosis. In this sense, we research novel mathematical descriptors for both morphology and dynamics and will compare their performance regarding small lesion classification based on novel feature selection algorithms. In addition, we develop and test a novel descriptor, the radial invariant Krawtchouk moments, to capture the morphology of diagnostically challenging lesions.
Therefore, the goal of our paper is to evaluate quantitatively based on automated classification the utility of novel feature extraction approaches to breast lesion detection such as foci and non-mass enhancing lesions.

Patients
A total of 63 patients, all females with an age range of 42 to 73 years, with 66 solid breast tumors were examined. All patients had histopathologically confirmed diagnosis from needle aspiration/excision biopsy and surgical removal.
Histologic findings were malignant in 30 and benign in 36 lesions. Tumors were classified as diagnostically challenging lesions for both foci and non-mass enhancing lesions. Lesion size was derived from mammography images. Mean size of malignant lesions was 1.3 cm (median = 1.2 cm, range = 0.4 to 3.8 cm), and mean size of benign lesions was 1.2 cm (median = 1.0 cm, range = 0.3 to 3.0 cm).

Magnetic resonance 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 prone position. First, transversal images were acquired with short TI inversion recovery sequence (TR = 5, 600 ms, TE = 60 ms, FA = 90 • , IT = 150 ms, matrix size 256 × 256 pixels, slice thickness 4 mm). A dynamic T1weighted gradient echo sequence (3D fast low-angle shot sequence) was performed (TR = 11 ms and TR = 9 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 or 2 mm.
The dynamic study consisted of six measurements with an interval of 83 s. The first frame was acquired before injection of a paramagnetic contrast agent (gadopentetate dimeglumine, 0.1 mmol/kg body weight, Magnevist™, Schering, Berlin, Germany), which was immediately followed by five other measurements. The initial localization of suspicious breast lesions was performed by computing difference images, i.e., subtracting the image data of the first from the fourth acquisition. As a pre-processing step to clustering, each raw gray level time series S(τ ), τ ∈ {1, . . . , 6} was transformed into a signal time series of relative signal enhancement x(τ ) for each voxel, the precontrast scan at τ = 1 serving as reference; in other words, (1) x (1) . Thus, we ensure that the proposed method is less sensitive to changing between different MR scanners and/or protocols.

Motion compensation and lesion segmentation
The employed motion compensation algorithm [18] is based on the Horn and Schunck method and represents a variational method for computing the displacement field, the so-called optical flow u, in an image sequence f 1 , f 2 with movement in between the image acquisitions: (1) The method is based on two typical assumptions for variational optical flow methods: the brightness and smoothness assumption. Additionally, we enforce a divergence-free flow field. The initial image sequence f 0 is pre-processed and convolved with a Gaussian K σ of a standard deviation σ : The variational method is based on the minimization of the continuous energy functional which penalizes all deviations from model assumptions: where the spatial domain represents a spatiotemporal domain. The data term (M(D k f , u)) describes a brightness constancy assumption, the first part of the regularizer (S(∇f , ∇u)) penalizes deviations from piecewise smoothness, and the second part of the regularizer D(divu) penalizes non-divergent-free flow fields. The weight terms μ > 0 and λ represent the regularization parameters, the so called Lamé constants, where larger values correspond to more simplified flow fields. This technique is a global method where the filling-in-effect yields dense flow fields, and no subsequent interpolation is necessary as with the technique proposed in [19]. This method works within a single variational framework. This technique overcomes the aperture problem, provides sub-pixel accuracy, and can be easily enhanced and adapted. The optimal motion correction results were achieved for motion compensation in two directions for mostly small standard deviations of the Gaussian kernel and smoothing parameter. Consistent with the only study known for evaluating the effect of motion correction algorithms [16], the proposed motion compensation technique improved the tumor diagnosis by having a strong impact on the enhancement curves.
Motion compensation is followed by image segmentation where each MR image has to be segmented into two regions, the region of interest, i.e., the voxels belonging to the tumor, and the background. We are using an interactive region growing algorithm that creates a binary mask for the tumor and its surroundings. The image used for the region growing algorithm was the difference image of the second post-contrast image and the native pre-contrast image.

Features describing the enhancement kinetics
The computer-assisted interpretation of time-signal series as measured during a dynamic contrast-enhanced magnetic resonance (DCE-MR) examination for each image voxel represents one of the major steps in designing CAD systems for breast MRI. 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 DCE-MR imaging. The results indicate that the enhancement kinetics, as represented by the time-signal intensity curves visualized in Figure 2, differ significantly for benign and malignant enhancing lesions and thus represent a basis for differential diagnosis. In cancerous breast tissue, plateau or washout time courses (type II or III) prevail. Steadily progressive signal intensity time courses (type I) are exhibited by benign enhancing lesions, albeit these enhancement kinetics are shared not only by benign tumors but also by fibrocystic changes.
As a dynamical feature, we use the slope of relative signal intensity enhancement (RSIE). To capture the correct form of the RSIE curve, we are fitting a straight line l = at +b to the data of the last three time scans of every pixel. The distance that l is supposed to minimize with respect to the data is the sum of the squared errors which leads to the well-known problem of least squares. Since we only need the slope a, it is sufficient to compute for the last three post-contrast scans: where l ij is the relative enhancement l ij := f ij −f 0j f 0j , n corresponds to the number of voxels we consider, and· denotes the mean.  [2]. Type I corresponds to a straight (Ia) or curved (Ib) line; enhancement continues over the entire dynamic study. Type II is a plateau curve with a sharp bend after the initial upstroke. Type III is a washout time course SIc−SI SI where SI is the pre-contrast signal intensity and SI c is the post-contrast signal intensity. In breast cancers, plateau or washout time courses (type II or III) prevail. Steadily progressive signal intensity time courses (type I) are exhibited by benign enhancing lesions. http://asp.eurasipjournals.com/content/2013/1/157 Since this feature has already proven its descriptive power, we are using it in the evaluation of motion correction.

Features describing the lesion morphology
Morphological characteristics contain valuable information about a lesion's type. Combined with kinetic properties, one could expect a higher accuracy. Furthermore, non-mass enhancing lesions such as DCIS or ICS can be better differentiated based on morphologic properties [20]. In the previous work [21], we have considered features that describe the geometric characteristics of the shape and local moments such as that of Krawtchouk to identify the non-smooth surface.

Minkowski functionals
Another family of characteristics that capture morphologic attributes is generated by Minkowski functionals (MFs) [22]. MF can characterize geometrical and topological space concepts such as shape, convexity, and connectivity and represent a simple yet precise tool for the analysis of geometrical structures in tumors. The most important MF for three-dimensional (3D) tumors is the Euler characteristic.
Since we consider 3D objects, Hadwiger's theorem [23] states that in this case we have four MFs, namely volume V, surface area S, mean breadth B, and Euler characteristic χ. The meaning of V and S are immediately clear, and the Euler characteristic is the number of regions of connected white voxels plus the number of completely enclosed regions of black voxels minus the number of black tunnels through white regions, remembering that in our data, the tumor is shown as a white object in a black area. The mean breadth is proportional to the integral mean curvature H of the tumor A which is defined as with κ 1 , κ 2 as the principal curvatures and df as the area element of A.
Since our object is represented as an MR image, we are only computing an approximation of the MFs of the existing tumor, but due to the nature of MRI, we cannot improve that. To calculate the MFs, one can use an iterative technique that checks every voxel and updates the values if the voxel belongs to the tumor. The efficient algorithm that we have used is presented in [22]. For the classification, it can also make sense to look at normed Minkowski functionals, which arẽ The N denotes the number of voxels in the tumor when we assume a voxel size of one.

Krawtchouk moments
The (discrete) Krawtchouk polynomials represent a set of functions that are related to the binomial distribution. Their descriptive power was shown in [24]. Since they are defined over a discrete set, their application is well suited in image processing.
The Krawtchouk polynomial of nth order is the hypergeometric function and p ∈ (0, 1), x, n = 0, 1, . . . , N, N > 0. 2 F 1 is defined as Their norm ρ(n; p, N) and their weight function w(x; p, N) are associated: To avoid numeric fluctuations, it is necessary to weight the Krawtchouk polynomials: Now, the polynomialsK n (x; p, N) form a set of orthonormal functions [25]: Since the computation of the weighted Krawtchouk moments via Equations 7 and 6 requires a lot of resources, one normally uses a recursive algorithm [24]: Being able to compute the Krawtchouk polynomials, we can define the Krawtchouk moment of the (n + m + l)th order for an image of size N × M × L as Since the weighted Krawtchouk polynomials form an orthonormal system, we can even reconstruct the image as For example, the Krawtchouk moments have successfully been used for content-based search by Mademlis et al. [26].

Radial invariant Krawtchouk moments
In this section, we are describing radial Krawtchouk moments in 3D based on the radial Krawtchouk moments in 2D presented in [25]. The idea behind this is to represent the image in a radial polar system f (r, θ). In an image of a size N × N, the radius r varies from 0 to N/2 and the angle θ varies from 0 to 2π in steps of ν. Therefore, θ c = 2πc ν , c = 0, . . . , ν − 1. Given this representation, the radial Krawtchouk moments are defined as u = N 2 is the number of pixels along the radius r. If the image is rotated by an angle α, we get Hence, the φ mn are rotationally invariant.
To extend this method straightforward to 3D, we need to find a representation of all rotations/orientations in special orthogonal group in 3D (SO (3) ) that uses only angles, and then we map them on a linear space having equidistant points and no singularities [27]. Since this cannot be accomplished, we are using the Euler angles as a representation of SO (3). They allow us to represent any three-dimensional rotation R ∈ SO (3)  The rotation is uniquely defined for every choice of angles with = 0. If = 0, we get a rotation of angle + around the Z-axis which can be achieved by an infinite combination of values for and . We avoid this in our discretization by shifting the values of by half of the step size of the discretization.
Let f (r, , , ) denote the image in the Euler coordinates, which means the vector (r, 0, 0) T is rotated by R , , and its three components describe the pixel at the corresponding x, y and z coordinates (of course, the origin of the coordinate system x, y, z refers to the center of the image), then the radial Krawtchouk moments in 3D are defined as Analogous to Equation 9, we establish the rotational invariant Krawtchouk moments in 3D as because if we have an image rotated arbitrarily in steps of the given discretization, there are numbers α, β, γ such that the moment R klmn (10) becomes

Classification techniques
The following section gives a description of classification methods applied to evaluate the effect of motion compensation to breast MRI 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 x describes a K-dimensional feature vector, that there are J classes, and there are N j samples available in group j. The mean in group j is given by μ j , and the covariance matrix is given by j .

Bayes classification based on LDA and QDA
The Bayes classification [28] is based on estimating the prior probabilities π i 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 posterior probability. Thus, the classification rule becomes (12) where μ j represents the means of the classes and j is the corresponding covariance matrix. The assignment to a certain class j for a certain pattern is made based on the smallest determined value of C j .
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 classifier, while if they are identical for the different groups, it becomes a LDA classifier.

Fisher's linear discriminant analysis
The underlying idea of Fisher's linear discriminant analysis (FLDA) is to determine the directions in the multivariate space that allow the best discrimination between the sample classes [29]. 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 J − 1 components, and for visualization purposes, scores of the first component are plotted against that of the second one, thus yielding the optimal separation among the classes.
This technique identifies the first discriminating component based on finding the vector a that maximizes the discrimination index given as a T B a a T W a (13) with B denoting the inter-class sum-of-squares matrix and W the intra-class sum-of-squares matrix.

Support vector machines
Support vector machines (SVMs) [30] represent an important technique for lesion classification in medical imaging.
The key point of this technique is to determine a hyperplane H = av + b that separates the feature vectors v i , i = 1, . . . , n in their d dimensional domain in two classes v i ∈ {M, B}. First, let us assume that our data set is linearly separable and that we can find a pair ( w, b) that fulfills Furthermore, we want both inequalities to be sharp. The best hyperplane is the one that maximizes the distance (margin) of the two parallel hyperplanes defined in Equation 14. Since the distance of a hyperplane to the origin is b w , we want to maximize 2 w . This leads to the following constrained optimization problem: find ( w, b) so that L( w) = w 2 is minimally subject to the condition Of course, this only works if all the variables are linearly separable, which cannot be assumed. This problem is solved by introducing positive slack variables θ i , i = 1, . . . , n which leads us to a minimization of L( w) An additional challenge appears if there is a nonlinear function that separates the variable since the common approach would fail under these circumstances. Aizerman et al. refined this method using non-linear kernel functions instead of the scalar product which maps the variables onto another space. The optimal hyperplane computed corresponds to a non-linear function in the original feature space.

Results
In the following, we will explore the applicability of the previously described motion compensation algorithm under different motion compensation parameters and features' sets as well as different classification techniques.  The results will elucidate the applicability of motion compensation as an artifact correction method and also the descriptive power of several tumor features with and without motion correction.

Experimental results
We analyze the effect of the computer-aided system based on feature extraction and classification with and without motion compensation with different parameter settings and for several extracted features describing kinetics, morphology, or both (Tables 1 and 2). Table 3 shows the results for different classifiers and single features when applied to tumor classification. Normed Minkowski functionals (independent of tumor size) such as the Euler characteristic achieve the highest AUC value among all other features. Thus, the descriptive power of this simple morphologic parameter is almost independent of the classifier's type. The regular Krawtchouk moment scores higher than the radial invariant one, and we assume that this occurs due to the re-discretization error and the usage of only one Krawtchouk polynomial. The area under the curve -representing the kinetic features -scored lower

Discussion
Breast MRI reading is often impeded by motion artifacts of different grades. Motion correction algorithms become a necessary correction tool in order to improve the diagnostic value of these breast MR images. We applied a new motion compensation algorithm based on the Horn and Schunck method for motion correction and determined the optimal parameters for lesion classification. At the same time, we described the lesions based on morphologic and kinetic features and determined the descriptive power of single features. In addition, we developed and tested a novel descriptor, the radial invariant Krawtchouk moment. The highest AUC value was achieved for the normed Euler characteristic showing that this topological measure captures well the shape of the tumor. Different classification techniques were applied to the classification of the lesions.

Conclusions
A novel CAD system was proposed for the detection of diagnostically challenging lesions. This system comprised a novel motion compensation algorithm based on the optical flow method, segmentation, morphological and kinetic feature extraction, and classification as an evaluation method for the extracted features from the lesions. Three new morphologic features (Minkowski functional, Krawtchouk moments, and invariant Krawtchouk moments) for the discrimination of diagnostically challenging lesions have been proposed in this article and compared with the standard enhancement kinetic features. The simulation results have shown that morphologic feature descriptors such as the Minkowski functional appear to be more adequate than kinetic descriptors for descriptors for diagnostically challenging lesions.