Evaluation of a Nonrigid Motion Compensation Technique Based on Spatiotemporal Features for Small Lesion Detection in Breast MRI

Abstract

Motion-induced artifacts represent a major problem in detection and diagnosis of breast cancer in dynamic contrast-enhanced magnetic resonance imaging. The goal of this paper is to evaluate the performance of a new nonrigid motion correction algorithm based on the optical flow method. For each of the small lesions, we extracted morphological and dynamical features describing both global and local shape, and kinetics behavior. In this paper, we compare the performance of each extracted feature set under consideration of several 2D or 3D motion compensation parameters for the differential diagnosis of enhancing lesions in breast MRI. Based on several simulation results, we determined the optimal motion compensation parameters. Our results have shown that motion compensation can improve the classification results. The results suggest that the computerized analysis system based on the non-rigid motion compensation technique and spatiotemporal features has the potential to increase the diagnostic accuracy of MRI mammography for small lesions and can be used as a basis for computer-aided diagnosis of breast cancer with MR mammography.

1. Introduction

Breast cancer is one of the leading death cases among women in the US. MR technology has advanced tremendously and became a highly sensitive method for the detection of invasive breast cancer [1]. While only dynamic signal intensity characteristics are integrated in today's CAD systems leaving morphological features to the interpretation of the radiologist, an automated diagnosis based on a combination of both should be strived for. Clinical studies have shown that combinations of different dynamic and morphologic characteristics [2] achieve diagnostic sensitivities up to 97% and specificities up to 76.5%.

However, most of these techniques have not been applied to small enhancing foci having a size of less than 1 cm. These diagnostically challenging cases found unclear ultrasound or mammography indications can be adequately visualized in magnetic resonance imaging (MRI) [3] with MRI providing an accurate estimation of invasive breast cancer tumor size [4].

Automatic motion correction represents an important prerequisite to a correct automated small lesion evaluation [5]. Motion artifacts are caused either by the relaxation of the pectoral muscle or involuntary patient motion and invalidate the assumption of 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 [6]. Therefore, spatial registration has to be performed before enhancement curve analysis. Due to the elasticity and heterogeneity of breast tissue, only nonrigid image registration methods are suitable. Although there has been a significant amount of research on nonrigid motion compensation techniques in brain imaging, few methods have been so far proposed for breast MRI. Most proposed techniques employ physically motivated deformation models [6, 7], transformations based on the deformation of B-splines, [8, 9], elastic transformations [10], and more recently adaptive grid generation algorithms [11]. We present in this paper a novel elastic image registration method based on the variational optical flow computation and will study its impact on the shape of the enhancement curves for small lesions.

The automated evaluation is a multistep system which includes a non-rigid motion compensation technique based on the optical flow, global and local feature extraction methods as shape descriptors, dynamical features, and spatiotemporal features combining both morphology and dynamics aspects to determine the optimal motion compensation parameters.

The paper is organized as follows. Section 2 describes the materials and methods, Section 3 the motion compensation technique, Section 4 the feature extraction techniques, and Section 5 the results. In Section 3, we present the non-rigid motion compensation technique, the feature extraction based on geometrical features, morphological features, dynamical features, and scaling index method for simultaneous time and space descriptions, and as classification technique the Fisher's discriminant analysis.

2. Material and Methods

2.1. Patients

A total of 40 patients, all female and age range 42–73, with indeterminate small mammographic breast lesions were examined. All patients were consecutively selected after clinical examinations, mammography in standard projections (craniocaudal and oblique mediolateral projections) and ultrasound. Only lesions BIRADS 3 and 4 were selected where at least one of the following criteria was present: nonpalpable lesion, previous surgery with intense scarring, or location difficult for biopsy (close to chest wall). All patients had histopathologically confirmed diagnosis from needle aspiration/excision biopsy and surgical removal. Breast cancer was diagnosed in 17 out of the total 31 cases. The average size of both benign and malignant tumors was less than 1.1 cm.

2.2. MR Imaging

MRI was performed with a 1.5 T system (Magnetom Vision, Siemens, Erlangen, Germany) with two different protocols equipped with a dedicated surface coil to enable simultaneous imaging of both breasts. The patients were placed in a prone position. First, transversal images were acquired with a STIR (short TI inversion recovery) sequence (TR = 5600 ms, TE = 60 ms, FA = 90^∘, IT = 150 ms, matrix size 256 × 256 pixels, slice thickness 4 mm). Then a dynamic T-weighted 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 mm or 2 mm.

The dynamic study consisted of 6 measurements with an interval of 83 s. The first frame was acquired before injection of paramagnetic contrast agent (gadopentetate dimeglumine, 0.1 mmol/kg body weight, Magnevist, Schering, Berlin, Germany) immediately followed by the 5 other measurements. The initial localization of suspicious breast lesions was performed by computing difference images, that is subtracting the image data of the first from the fourth acquisition. As a preprocessing step to clustering, each raw gray level time-series 𝑆(𝜏),𝜏∈{1,…,6} was transformed into a signal time-series of relative signal enhancement 𝑥(𝜏) for each voxel, the precontrast scan at 𝜏=1 serving as reference. Thus, we ensure that the proposed method is less sensitive to changing between different MR scanners and/or protocols.

3. Motion Compensation Technique

The employed motion compensation algorithm is based on the Horn and Schunck method [12] and represents a variational method for computing the displacement field, the so-called optical flow, in an image sequence. It is based on two typical assumptions for variational optical flow methods, the brightness constancy, and smoothness assumption.

In this context, the MR image sequence 𝑓₀ is a differentiable function of brightness values on a four-dimensional spatiotemporal image domain Ω: 𝑓₀∶Ω_⊂ℝ3×ℝ₊→ℝ₊.(1)

From this image sequence, we want to compute a dense vector field 𝐮=(𝑢₁,𝑢₂,𝑢₃)^⊤∶Ω→ℝ^{2,3} that describes the motion between the precontrast image at time point 𝑡 and a postcontrast image at time point 𝑡+𝑘, either for all three dimensions (ℝ³) or only in one transversal slice (ℝ²).

The initial image sequence 𝑓₀ is preprocessed and convolved with a Gaussian 𝐾_𝜎 of a standard deviation 𝜎 to remove noise in the image: 𝑓=𝐾_𝜎∗𝑓₀.(2)

The brightness constancy assumption dictates that under the motion 𝐮, the image brightness values of the precontrast image at time 𝑡 and the postcontrast image at time 𝑡+𝑘 remain constant in every pixel: 𝑓(𝐱+𝐮(𝐱),𝑡+𝑘)=𝑓(𝐱,𝑡),∀𝐱∈Ω.(3) Naturally, this condition by itself is not sufficient to describe the motion field 𝐮 properly, since for a brightness value in an image voxel in the precontrast image, there are generally many voxels in the postcontrast image with the same brightness value, or, in the presence of noise, possibly even none at all. Therefore, we include the smoothness assumption, dictating that neighboring voxels should move in the same direction, which is expressed as the gradient magnitude of the flow field components is supposed 0 to be as follows: ||∇𝑢_{1,2,3}(𝐱)||=0,∀𝐱∈Ω.(4) This constraint by itself would force the motion field to be a rigid translation, which is not the case in MR images. However, if we use both the brightness constancy assumption and the smoothness assumption as weak constraints in an energy formulation, the motion field 𝐮 that minimizes this energy matches the postcontrast image to the precontrast image and is spatially smooth.

The variational method is based on the minimization of the continuous energy functional which penalizes all deviations from model assumptions: 𝐸(𝐮)=_Ω(𝑓(𝐱+𝐮(𝐱),𝑡+𝑘)−𝑓(𝐱,𝑡))²_{Dataterm}+𝛼||∇𝑢₁(𝐱)||²+||∇𝑢₂(𝐱)||²+||∇𝑢₃(𝐱)||²d𝐱_{Smoothnessterm}.(5)

The weight term 𝛼>0 represents the regularization parameter where larger values correspond to smoother 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 [13]. This method works within a single variational framework.

Given the computed motion from the precontrast image to a postcontrast image, the postcontrast image is being registered backwards before its difference image with the precontrast image is being computed for tumor classification: 𝑓_{post-registered}(𝐱)=𝑓_post(𝐱+𝐮(𝐱)).(6) Since this registration explicitly yields the brightness constancy assumption, if the motion has been estimated correctly, the registered postcontrast image is equal to the precontrast image.

Breast MR images are mostly characterized by brightness, since bright regions are created by fatty tissue or contrast agent enhancing tumor tissue, while dark regions describe glandular tissue and background. Tumors are mainly located in the glandular tissue and proliferate either into the fatty tissue (invasive) or along the boundary inside the glandular tissue (noninvasive).

Two major concerns have to be addressed when applying the optical flow approach to breast MRI: (1) The constancy assumption does not hold for objects appearing from one image to the next such as lesions the contrast agent enhancement of which is much stronger than in the surrounding tissue and (2) The lack of constant grid size in all directions since voxel size is smaller than the slice thickness. The first concern is alleviated by masking suspicious areas by a radiologist and by detecting the sharp gradients in the motion field in the unmasked image. Figure 1 shows the masking of the entire upper image and visualizes the motion in the inner slice by a color code. This color code describes the motion direction based on the hue and the motion magnitude based on the brightness and thus identifies suspicious regions by detecting the sharp gradients. The mask 𝑚∶Ω→{0,1} can be easily incorporated in the energy formulation and it forces the data term to disappear in suspicious regions: 𝐸(𝐮)=_Ω𝑚(𝐱)(𝑓(𝐱+𝐮(𝐱),𝑡+𝑘)−𝑓(𝐱,𝑡))²+𝛼||∇𝑢₁(𝐱)||²+||∇𝑢₂(𝐱)||²+||∇𝑢₃(𝐱)||²d𝐱.(7) In addition, there can be gaps between the slices where no nuclei are being excited in order to avoid overlapping of the slices. Figure 2 shows an example of a tumor considerably shifting its position on adjacent transversal slices.

(a)

(b)

(c)

(d)

Figure 1 

Motion detection on a transverse image. (a) Masking the data term: the green lines separate the boundary between masked and unmasked areas. (b) Color code describing motion from the interior of the image. (c) Motion in two directions determined without mask and (d) based on the mask from (a). The values for the standard deviation of the Gaussian presmoothing kernel and for the smoothness term are 𝜎=3 and 𝛼=500.

(a)

(b)

(c)

(d)

Figure 2 

Tumor in adjacent transverse slices of a 512 × 512 × 32 image.

To overcome the second concern, it is important to decide whether the motion in transverse direction is having a significant impact. The present research has shown that there is no significant difference in visual quality if motion is computed in two or three directions. In the following notation, we will consider the motion in three directions, the one in two directions is analogous.

Based on the work of Horn and Schunck, a vast variety of research in optical flow estimation has been conducted, most of it on more robust, edge preserving regularity terms, for example [14]. In our work however, we favor the original quadratic formulation, since we explicitly need the filling-in effect of a non-robust regularizer to fill in the information in masked regions. To overcome the problem of having a non-convex energy in (5), we use the coarse-to-fine warping scheme detailed in [14], which linearizes the data term as in [12] and computes incremental solutions on different image scales.

For this, we approximate the image 𝑓 at the distorted point (𝐱+𝐮(𝐱),𝑡+𝑘) by a first-order Taylor approximation: 𝑓(𝐱+𝐮(𝐱),𝑡+𝑘)−𝑓(𝐱,𝑡)≈𝑓(𝐱,𝑡+𝑘)+∇𝑓(𝐱,𝑡+𝑘)^⊤𝐮(𝐱)−𝑓(𝐱,𝑡).(8) Alternatively, one can develop the Taylor series at time 𝑡, getting 𝑓(𝐱+𝐮(𝐱),𝑡+𝑘)−𝑓(𝐱,𝑡)≈𝑓(𝐱,𝑡)+∇𝑓(𝐱,𝑡)^⊤𝐮(𝐱)+𝜕𝜕𝑡𝑓(𝐱,𝑡)𝑘−𝑓(𝐱,𝑡).(9) Since the term 𝑓(𝐱,𝑡) cancels and the temporal derivative is again approximated by the difference between the two images at point 𝐱, the only difference between (8) and (9) is the time at which the spatial derivative ∇𝑓 is being computed. In optical flow computation, one usually uses the arithmetic mean (1/2)(∇𝑓(𝐱,𝑡+𝑘)+∇𝑓(𝐱,𝑡)). For the purpose of registration, the scalar factor 𝑘 can be neglected since it is arbitrary whether one computes a motion for a 𝑘≠1 and then scales the motion later on with 𝑘 when registering the image, or simply sets 𝑘 to 1. Incorporating this into the energy formulation and leaving out the indices for better readability, the linearized energy functional then reads as 𝐸_lin(𝐮)=_Ω𝑚𝜕𝜕𝑡𝑓+∇𝑓^⊤𝐮²+𝛼||∇𝑢₁||²+||∇𝑢₂||²+||∇𝑢₃||²𝑑𝐱.(10) This functional is convex in 𝐮, and a minimizer can be found by solving its Euler-Lagrange equations: 0=𝑚𝜕𝑓𝜕𝑥𝜕𝑓𝜕𝑥𝑢₁+𝜕𝑓𝜕𝑥𝜕𝑓𝜕𝑦𝑢₂+𝜕𝑓𝜕𝑥𝜕𝑓𝜕𝑧𝑢₃+𝜕𝑓𝜕𝑥𝜕𝑓𝜕𝑡−𝛼Δ𝑢₁∀𝐱∈Ω,0=𝑚𝜕𝑓𝜕𝑥𝜕𝑓𝜕𝑦𝑢₁+𝜕𝑓𝜕𝑦𝜕𝑓𝜕𝑦𝑢₂+𝜕𝑓𝜕𝑦𝜕𝑓𝜕𝑧𝑢₃+𝜕𝑓𝜕𝑦𝜕𝑓𝜕𝑡−𝛼Δ𝑢₂∀𝐱∈Ω,0=𝑚𝜕𝑓𝜕𝑥𝜕𝑓𝜕𝑧𝑢₁+𝜕𝑓𝜕𝑦𝜕𝑓𝜕𝑧𝑢₂+𝜕𝑓𝜕𝑧𝜕𝑓𝜕𝑧𝑢₃+𝜕𝑓𝜕𝑧𝜕𝑓𝜕𝑡−𝛼Δ𝑢₃∀𝐱∈Ω.(11) Since the linearization in (8) or (9) is only valid for small motions in the subpixel range, a typical strategy to overcome the problem of large motions is to downsample the MR images to a coarse resolution, compute an approximate motion on the coarse resolution, interpolate this motion to the next finer resolution, register the second image with the approximate motion, compute the incremental motion from the first image to the registered second image, add the incremental motion to the approximate motion, and repeat this iteration up to the original resolution.

4. Segmentation

Before we proceed with feature extraction, each MR image has to be segmented into two regions, the region of interest (ROI), that is, the voxels belonging to the tumor, and the background. We are using an interactive region growing algorithm, creating a binary mask the tumor voxels of which are true, and all other voxels are false. The image used for the region growing algorithm was the difference image of the second postcontrast image and the native precontrast image. The center of the lesion was interactively marked on one slice of the subtraction images and then a region growing algorithm included all adjacent contrast-enhancing voxels also those from neighboring slices. Thus a 3D form of the lesion was determined. An interactive ROI was necessary whenever the lesion was connected with diffuse contrast enhancement, as it is the case in mastopathic tissue.

Figure 3 shows a transverse image of a tumor in the right breast and its binary segmentation, created with region growing.

(a)

(b)

Figure 3 

Example of an MR images segmentation showing the transverse image of a tumor in the right breast (a) and the binary segmentation of the tumor (b).

5. Feature Extraction

The complexity of the spatio-temporal tumor representation requires specific morphology and/or kinetic descriptors. We analyzed geometric and dynamical features as shape descriptors, provided a temporal enhancement modeling for kinetic feature extraction and the scaling index method for the simultaneous morphological and dynamics representation.

5.1. Contour Features

To represent the shape of a tumor, we compute the following features: the maximal, the minimum, the mean, and the standard deviation of a radius as well as the azimuth, entropy, and compactness.

First, we determine the center of mass of the tumor as 𝑣∶=1𝑛^𝑛_𝑖=1𝑣_𝑖,(12) where 𝑛 is the number of voxels belonging to the tumor, and 𝑣_𝑖 is location of the 𝑖th voxel.

We define for each contour voxel 𝑐_𝑖, the radius 𝑟_𝑖 and the azimuth 𝜔_𝑖 (i.e, the angle between the vector from the center of mass to the voxel 𝑐_𝑖 and the sagittal plane) and give in the following the definitions of the geometrical features to be used in context with the motion compensation technique.

From the set of the values 𝑟₁,…,𝑟_𝑚 representing the radii, we determine the minimum value 𝑟_min and the maximum value 𝑟_max as well as

the radius: 𝑟_𝑖∶=‖‖𝑐_𝑖−𝑣‖‖₂,(13) the azimuth: 𝜔_𝑖∶=arcsine⎛⎜⎜⎜⎝𝑐_𝑖𝑥−𝑣_𝑥𝑐_𝑖𝑥−𝑣_𝑥²+𝑐_𝑖𝑦−𝑣_𝑦²⎞⎟⎟⎟⎠,(14) the mean value: 𝑟∶=1𝑚^𝑚_𝑖=1𝑟_𝑖(15) the standard deviation: 𝜎_𝑟∶=⎷1𝑚^𝑚_𝑖=1𝑟_𝑖−𝑟²,(16) and the entropy: ℎ_𝑟∶=−¹⁰⁰_𝑖=1𝑝_𝑖⋅log₂𝑝_𝑖.(17) The subscripts 𝑥 and 𝑦 denote the position of the voxel in sagittal and coronal direction respectively. 𝜔_𝑖 was also extended to the range from −𝜋 to 𝜋 by taking into account the sign of (𝑐_𝑖𝑦−𝑣_𝑦).

The entropy ℎ_𝑟 was computed from the normalized distribution of the values into 100 “buckets,” where 𝑝_𝑖 is defined as follows:

For 0 ≤ i ≤ 99: 𝑝_𝑖∶=||𝑟_𝑗∣𝑖≤𝑟_𝑗−𝑟_min/𝑟_max−𝑟_min⋅100<𝑖+1||𝑚.(18)

From the radius, 𝑟_min, 𝑟_max, 𝑟, 𝜎_𝑟, and ℎ_𝑟 were used as morphological features of the tumor. From the azimuth, only the entropy ℎ_𝜔 (computed for 𝜔 as in (17) and (18)) was used as a feature, since the values 𝜔_min and 𝜔_max are always around 𝜋 and −𝜋, and the values 𝜔, and 𝜎_𝜔 are not invariant under rotation of the tumor image.

An additional measurement describing the compactness of the tumor, which was also used as a feature, is the number of contour voxels, divided by the number of all voxels belonging to the tumor.

5.2. Geometric Moments

The spatial and morphological variations of a tumor can be easily captured by shape descriptors. Here, we will employ geometric moments [15] as shape descriptors for our lesions. For each lesion, we determine 6 three-dimensional moments that are invariant under rotation, translation, and scaling and are able to completely characterize the shape of the tumor.

Geometric moments are defined in the discrete three-dimensional case as 𝑚_pqr=_𝑥_𝑦𝑥^𝑝𝑦^𝑞𝑧^𝑟𝑓(𝑥,𝑦,𝑧),(19) with 𝑓(𝑥,𝑦,𝑧) being the image intensity function (gray values). Finally, we compute three-dimensional moment invariants [16] that are not computational intensive and provide results stable to noise and distortion.

5.3. Dynamical Features

For each lesion, the time/intensity curve is computed based on the average signal intensity of the tumor before contrast agent administration (SI) and after contrast agent administration (SI_𝐶), the relative enhancement as ΔSI∶=SI_𝐶−SISI⋅100%.(20)

To capture the slope of the curve of relative signal intensity enhancement (RSIE) versus time in the late postcontrast time, we compute the linear approximation of the RSIE curve for the last three time scans. The slope of this curve represents an important parameter for lesion classification.

5.4. Simultaneous Morphology and Dynamics Representations

The scaling index method [17] describes mathematically an estimation of the local scaling properties of a point represented in an 𝑛-dimensional embedding space. For every single point in the distribution, the number 𝑁 of points is determined, which are included in an 𝑛-dimensional sphere of a radius 𝑟 centered at 𝐱_𝐢𝑁𝐱_𝐢,𝑟=_𝑗Θ𝑟−||𝐱_𝐣−𝐱_𝐢||,(21) where Θ represents the Heaviside function. The function 𝑁(𝐱_𝐢,𝑟), is usually determined within a specified range 𝑟∈[𝑟₁,𝑟₂], the so-called scaling range.

If the scaling region is large, then 𝑁(𝐱_𝐢,𝑟) becomes 𝑁(𝐱_𝐢,𝑟)∝𝐫^𝛼 where 𝛼 is the scaling index. A first-order approximation yields 𝛼_𝑖=log𝑁𝐱_𝐢,𝑟₂−log𝑁𝐱_𝐢,𝑟₁log𝑟₂−log𝑟₁(22) with 𝑟₁ and 𝑟₂ being the lower and upper limit of the scaling range.

This technique can easily distinguish between point-like (𝛼=0), string-like (𝛼=1), and sheet-like (𝛼=2) structures, in images and be thus employed for texture extraction. In the context of breast MRI, such a point consists of the sagittal, coronal, and transverse positions of a tumor voxel and its third-time scan gray value, and the scaling index serves as an approximation of the dimension of local point distributions.

The scaling index serves as a descriptor of the simultaneous enhancement kinetics and morphology. In particular, it can capture the blooming and the hook sign [18] and identify cancers with benign-like kinetics, discriminate normal tissue showing cancer-like enhancement curves, and thus improve the performance of computer-assisted detection of small enhancing lesions.

5.5. Morphologic Blooming

Another feature relying on morphology as well as on dynamics is the so-called blooming: the tumor contour becoming blurred and proliferating over the time of contrast agent enhancement is an indication for malignancy [19]. Beside blooming, [PTL+ 06] also attributes “centripetal” enhancement (strong, contrast agent enhancement near the tumor contour, weak enhancement in the tumor center with advancing time), and inhomogeneous contrast agent enhancement in general to a high probability of malignancy of the tumor.

To capture the blooming of a tumor, the difference of the relative signal intensity enhancement (RSIE) of every contour voxel and the average RSIE of its nontumor neighbor voxels was computed, and their mean, standard deviation and entropy from the 3rd to the 5th postcontrast image were used as features. Unlike for computation of the radial transform, for these features not only voxels in the contour of each transverse slice were considered contour voxels, but also voxels of the tumor with a non-tumor neighbor voxel on an adjacent transverse slice. The minimum and maximum values were neglected as features, since the gray values in the images are sensitive to noise.

The general irregularity of contrast agent enhancement in the tumor interior was computed as the standard deviation and entropy of the RSIE 𝑠_𝑖,𝑗 described above for the time scans 1, 3, and 5 (𝑖∈{1,3,5}), and for all tumor voxels 𝑗. For each of the three time scans, the standard deviation and entropy were used as a feature.

In order to compensate general contrast and brightness differences between the images, for all features described in this section, the images were normalized to have the same mean as the precontrast image, and the standard deviation was divided by the standard deviation of the precontrast image.

6. Results

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

Table 1 describes the motion compensation parameters used in the subsequent evaluations.

Table 1 

Motion compensation parameters for two or three directions. 𝜎 represents the standard deviation for presmoothing and 𝛼 the regularization parameter.

The effect of the motion compensation based on motion compensation parameters such as the amount of presmoothing and the regularization parameter was analyzed based on different combinations of feature groups.

We analyzed the performance of each proposed feature set and combinations of those for tumor classification based on receiver-operating characteristic (ROC) and compare the area-under-the-curve (AUC) values.

As a classification method to evaluate the effect of motion compensation on breast MR images, we chose the Fisher's linear discriminant analysis. Discriminant analysis represents an important area of multivariate statistics and finds a wide application in medical imaging problems. The most known approaches are linear discriminant analysis (LDA), quadratic discriminant analysis (QDA), and Fisher's canonical discriminant analysis.

Let us assume that 𝐱 describes a K-dimensional feature vector, that there are 𝐽 classes and there are 𝑁_𝑗 samples available in group 𝑗. The mean in group 𝑗 is given by 𝜇_𝑗, and the covariance matrix is given by Σ_𝑗.

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. 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 𝐽−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 𝐚 that maximizes the discrimination index given as 𝐚^𝐓𝐁𝐚𝐚^𝐓𝐖𝐚,(23) with 𝐁 denoting the interclass sum-of-squares matrix and 𝐖 the intraclass sum-of-squares matrix.

6.1. Effectiveness of Contour Features

The contour features described in Section 5.1 show for almost all motion compensation parameters high ROC values as shown in Table 2. Without motion compensation, the entropy as well as radius mean and maximum yield the best results. The radius minimum followed by the compactness show a significant improvement for motion compensation in 3D directions.

Table 2 

Areas under the ROC curves for contour features using FLDA. The rows represent the motion compensation as given in Table 1. Numbers in boldface show the best results.

6.2. Effectiveness of Geometric Moments

The geometric moments describe a representation of average shape parameters. Table 3 shows that motion compensations in two directions yields better results than in three directions for almost all geometric moments. However, the geometric moments seen not to be an adequate mean to extract features for small lesions.

Table 3 

Areas under the ROC curves for morphological features (geometric moments 𝐼₁ to 𝐼₆) using FLDA. The rows represent the motion compensation as given in Table 1. Numbers in boldface show the best results.

6.3. Effectiveness of Kinetic Features

The slope is derived from the first-order approximation of relative signal intensity enhancement from the last three scans. Table 4 shows that both 2D as well as 3D motion compensation yield almost equally good results.

Table 4 

Areas under the ROC curves for dynamic features (slope) using FLDA. The rows represent the motion compensation as given by Table 1.

6.4. Effectiveness of Spatiotemporal Features

For the two radii 𝑟₁ and 𝑟₂, we used radii in the size of tumor structures, 𝑟₁=3mm and 𝑟₂=6mm. The maximum, mean, and standard deviation and entropy of the set of scaling-indices, computed from tumor points as in (15)–(17), were used as features of the tumor. The minimum was neglected, since for almost every tumor it was 0, due to isolated points.

Table 5 shows the ROC values for the scaling index method for different motion compensation parameters shown in Table 1. The scaling index mean value yields the highest results without motion compensation and for both 2D and 3D motion compensation.

Table 5 

Areas under the ROC curves for scaling index (SI) method using FLDA. The rows represent the motion compensation as given by Table 1. Numbers in boldface show the best results.

The performance results for the spatio-temporal features related to both contour and tumor relative signal intensity enhancement are shown in Table 6 for the third, fourth, and fifth scans. For both the tumor and contour, the entropy showed the best results in the ROC analysis: the third scan for the tumor entropy and the fourth for the contour entropy for both uncompensated and compensated motion.

Table 6 

Areas under the ROC curves for spatio-temporal features using FLDA. The rows represent the motion compensation as given in Table 1.

6.5. Effectiveness of All Combined Features

In a final step, we considered all determined features as a vector, reduced its dimension with PCA and used FLDA as a classification method. The AUCvalues for the optimal number of PCs and depending on the motion compensation parameters are shown in Table 7. We observe again that motion compensation in 2D yields promising results for a CAD for small lesions.

Table 7 

Areas under the ROC curves for all combined features and for the optimal number of PCs using FLDA. The rows represent the motion compensation as given by Table 1.

7. Conclusion

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 mammographic 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. Several novel lesion descriptors such as morphologic, kinetic, and spatiotemporal are applied and evaluated in context with benign and malignant lesion discrimination. 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 [20], the proposed motion compensation technique achieved good results for weak motion artifacts.

The performed ROC-analysis shows that an integrated motion compensation step in a CAD system represents a valuable tool for supporting radiological diagnosis in dynamic breast MR imaging.

Acknowledgments

The research was supported by NIH Grant 5K25CA106799-04.