**2. Methodology**

The proposed methodology follows the CADx workflow presented in the previous section. However, asymmetry measurements are used to aid in the diagnosis. To obtain such measurements, two additional stages are incorporated into the workflow: registration and pixelwise subtraction. Additionally, a series of image transformations are incorporated to enhance different characteristics of the breast in the mammograms. This work is based on and follows previous efforts [30–32].

**Figure 2** shows how the bilateral asymmetry information was incorporated into the CADx system. Briefly, soft tissue is first segmented, the image of the left breast is then registered to its right counterpart and a bilateral subtraction of the co-registered images is performed;

**Figure 2.** Workflow of the proposed methodology.

images are then filtered and features are extracted; a multivariate model is selected using a train set; and finally, the model is evaluated on a validation set. A detailed explanation of each stage is presented in the following sections.

#### **2.1. Materials**

Ferrari et al. [26] characterized asymmetry as variations in oriented textural patterns, obtained using directional filtering with Gabor wavelets at different orientations and scales. Using a

Rodriguez-Rojas et al. [21] presented a CADx system targeted to detect high-risk cancer patients. To do so, automated breast tissue segmentations were performed on 200 Mexican subjects labeled as either low- or high-risk according to their BI-RADS score. Then, 50 features were extracted, and bilateral differences between mammograms were defined by subtracting corresponding features in both mammograms. Finally, a genetic algorithm selected a predictive combination of features. Using this methodology, they were able to classify low-risk and high-risk cases with an area under the receiver operating characteristic (ROC) curve (AUC) of 0.88 on a 150-fold cross-validation set. The features included in the model were associated

In summary, and as presented, most asymmetry detection methods are either feature-based, rely on simple bilateral subtraction techniques [14, 27], or depend on an ROI provided by a radiologist [24, 25]. Thus, in order to efficiently measure asymmetry, a better and automatic registration must be performed [28]. To do so, alignment has been improved by using the nipple as a reference point [29] and by co-registering both breasts using a robust point matching approach [22]. Nevertheless, none of those works include a fully automated bilateral registration. In this chapter, a methodology that incorporates an automatic asymmetry analysis with both a feature-based and a pixel-wise bilateral subtraction into a CADx system is presented.

The proposed methodology follows the CADx workflow presented in the previous section. However, asymmetry measurements are used to aid in the diagnosis. To obtain such measurements, two additional stages are incorporated into the workflow: registration and pixelwise subtraction. Additionally, a series of image transformations are incorporated to enhance different characteristics of the breast in the mammograms. This work is based on and follows

**Figure 2** shows how the bilateral asymmetry information was incorporated into the CADx system. Briefly, soft tissue is first segmented, the image of the left breast is then registered to its right counterpart and a bilateral subtraction of the co-registered images is performed;

database with 80 images resulted in a classification accuracy of up to 0.744.

with the differences in signal distribution and tissue shape.

**2. Methodology**

112 New Perspectives in Breast Imaging

previous efforts [30–32].

**Figure 2.** Workflow of the proposed methodology.

A total of 1796 digitalized film mammograms from 449 different subjects were used. From those, 45 were classified as healthy subjects (HS) (mean age of 59.3 and standard deviation (SD) of 9.8 years), 197 as subjects with malignant calcifications (CS) (mean age of 58 and SD of 10.9 years), and 207 as subjects with malignant masses (MS) (mean age of 64.1 and SD of 10.1 years). Each subject had the four standard mammograms taken, namely, left and right craniocaudal (CC), and left and right mediolateral oblique (MLO) projections.

In order to avoid problems associated with intra-scanner variability [17, 22, 33], all mammograms in this study were obtained from the Howtek dataset of the Digital Database for Screening Mammography public database [34], in which all mammograms were digitalized using a Howtek 960 scanner using a sampling rate of 43.5 micrometers per pixel and a 12-bit depth.

#### **2.2. Segmentation**

Segmentation, also called categorization by computer vision definitions, allows delimiting one or several parts of a given image assigning one class label (e.g. bone, muscle, fat, skin, calcification, and mass). This process is defined by the division or segmentation of the image into several homogeneous regions disjointed from their surroundings. A commonly used automatic segmentation of the breast tissue is based on the estimation of the background noise. For this study, an initial segmentation mask was created by estimating the background noise in the image and discarding all pixels below five standard deviations of the noise level. Then, holes were removed by applying closing morphological operations with a 3 × 3 supporting region, as described by Eq. (1):

$$\mathcal{S}(A) = \left( A(\mathbf{x}, \mathcal{y}) \oplus \mathcal{B}(\mathbf{x}, \mathcal{y}) \right) \ominus \mathcal{B}(\mathbf{x}, \mathcal{y}) \tag{1}$$

where ⊕ and ⊖ are the grayscale dilation and erosion morphological operations, respectively. B *(x, y)* is a 3 × 3 structural element. *A (x, y)* is the image being segmented and *S (A)* is the resulting segmentation of the *A (x, y)* image. The largest connected region is used as the segmentation mask while all other high-intensity regions are removed from the images. **Figure 3** shows an example of the results of the segmentation procedure.

#### **2.3. Registration**

Image registration can be defined as the intensity and spatial mapping between two images [35]. Given two input images *F* and *M*, image registration can be expressed as *R*′ = *g*[*T*(*F*)], where *T* is a spatial transformation function, *g* an intensity transformation function, and *R*′ the registered image. The transformation function is not always necessary; a lookup table can be used to pinpoint intensities. A visual example of image registration is presented in **Figure 4**, where an image *M* is being registered to match image *F*.

**Figure 3.** Segmentation of breast tissue. The image on the left is the original CC mammogram and the image on the right shows the superimposed segmentation mask in white (image from Ref. [32]).

Image registration has been widely used in medical applications [28, 36, 37]. However, the soft nature of the breast tissue makes them highly deformable, and rigid registration procedures, in which only rotation, translation, and scaling functions are used, are not sufficient. Therefore, nonrigid registration methods are necessary [38, 39]. There are many approaches to deal with medical imaging registration, the most recent comparison of algorithms based on a retrospective evaluation was published by West et al. [40], but it was constrained to do intrapatient rigid registration. Also recently, Diez et al. [28] and Celaya-Padilla et al. [30] compared registration algorithms with breast images as a source, and both concluded that the B-Splines approach was the most consistent.

Breast image registration based on a B-Splines transformation is defined as follows: given two input images (*F* = target image, *M* = image being registered), *M* is deformed by modifying a mesh of control points following a maximization of a similarity measure based on steepest descent gradient [6, 15]. The deformed image is compared to *F* using a similarity metric. If the images are similar enough, the process stops. Otherwise, the process reiterates.

**Figure 4.** Basic example of an image registration procedure. *F* is the target image, *M* is the image to be registered, and *R*′ is the registered image.

**Figure 5** shows a multi-resolution pyramid approach [41] for the B-Spline implementation. There, the images are first registered using low-resolution images, the B-Spline transformation parameters are moved into the next higher resolution and parameter optimization is run again, and so on. This often avoids issues with local minima in the parameter search space and reduces computational time [15].

For this study, the image to be registered was first horizontally flipped. Then, both the moving image and the target image were resampled into a lower resolution image. Next, the pyramids for the multi-resolution were generated. Afterwards, the registration process detailed in **Figure 5** was carried out. And finally, the original moving image was deformed using the final parameters of the registration. For this implementation, mutual information [39] was used as the similarity metric. In **Figure 6**, the checkerboard of an example result from the B-Spline registration procedure is presented. There, it can be seen that the registered image was successfully aligned with its counterpart.

#### **2.4. Image subtraction**

Image registration has been widely used in medical applications [28, 36, 37]. However, the soft nature of the breast tissue makes them highly deformable, and rigid registration procedures, in which only rotation, translation, and scaling functions are used, are not sufficient. Therefore, nonrigid registration methods are necessary [38, 39]. There are many approaches to deal with medical imaging registration, the most recent comparison of algorithms based on a retrospective evaluation was published by West et al. [40], but it was constrained to do intrapatient rigid registration. Also recently, Diez et al. [28] and Celaya-Padilla et al. [30] compared registration algorithms with breast images as a source, and both concluded that the B-Splines

**Figure 3.** Segmentation of breast tissue. The image on the left is the original CC mammogram and the image on the right

shows the superimposed segmentation mask in white (image from Ref. [32]).

Breast image registration based on a B-Splines transformation is defined as follows: given two input images (*F* = target image, *M* = image being registered), *M* is deformed by modifying a mesh of control points following a maximization of a similarity measure based on steepest descent gradient [6, 15]. The deformed image is compared to *F* using a similarity metric. If the

**Figure 4.** Basic example of an image registration procedure. *F* is the target image, *M* is the image to be registered, and

images are similar enough, the process stops. Otherwise, the process reiterates.

approach was the most consistent.

114 New Perspectives in Breast Imaging

*R*′ is the registered image.

Once the images were co-registered, a pixel-wise absolute difference was computed between the left and right images, as defined by Eq. (2) as follows:

$$I\_{\Delta}(\mathbf{x}, y) = \left| I\_{\prime}(\mathbf{x}, y) - I\_{\prime}(T(\mathbf{x}, y)) \right| \tag{2}$$

where *Ir* (*x, y*) represents the right image, *Il* (**T**(*x, y*)) represents the left image registered to the right image space, and *I*∆(*x, y*) represents the map of absolute differences. **Figure 7** shows an example of the differential image for two given input images.

**Figure 5.** B-Spline registration typical framework.

**Figure 6.** Checkerboard comparison of images *pre* and *post* B-Spline registration. The image in the left shows a comparison between a left and horizontally flipped right breast before registration, and the right image shows the results of the registering process.

**Figure 7.** Image subtraction example. Left: unaltered CC view of left breast, middle: horizontally flipped CC view of right breast, and right: color map of the subtraction image *I* ∆. White and black pixels inside the breast tissue represent small and large intensity differences, respectively.

#### **2.5. Image enhancement**

To study the appearance of the architectural distortions, two enhancing filters were applied to the images: a morphological high-frequency enhancement filter (*H*) designed to enhance fiber-like tissues, and a Laplacian of Gaussian filter (*L*) that enhances high-frequency patterns inside the breast tissue. Additionally, since the texture between normal and abnormal tissues is different [42], two texture maps were created. The first map computed the local standard deviation (*σ*) of the mammograms, and the second map computed the local fractal dimension (*F*). All image processing was implemented in C++ using Insight Segmentation and Registration Toolkit (ITK) libraries for image manipulation following previous efforts [32, 43].

#### **2.6. Feature extraction**

There are several features that may be quantified when aiming to detect early signs of cancer. For this analysis, 43 features were extracted from each image. These features can be grouped in three main categories: shape (i.e. area, perimeter, compactness, elongation, region centroid, region scatter), signal (i.e. mean, median, energy, variance, standard deviation, dynamic range, *z* mean, entropy, skewness, kurtosis, *z* range, fraction greater than *z* deviations, fraction lower than *z* deviations, value at fraction, 5% trimmed mean, 5% trimmed standard deviation, 5% trimmed *z M*ean), and morphology (i.e. total signal, signal centroid, signal scatter, and signal surface). Details of the full feature extraction procedure can be found in Ref. [32].

The enhancement filters and texture maps presented in Section 2.5 were applied to the four screening mammograms (i.e. left and right CC, and left and right MLO) and to the two bilateral subtraction images (CC and MLO), yielding a set of 15 images for both the CC and the MLO views: *I r* , *I l* , *I <sup>Δ</sup>*, *H*<sup>r</sup> , *Hl* , *HΔ*, *Lr* , *Ll* , *LΔ*, *σ<sup>r</sup>* , *σl* , *σΔ*, *Fr* , *Fl* , and *FΔ*, where *I* is the raw image, *H, L, σ*, and *F* are the enhanced images described in Section 2.5, and *r, l*, and Δ, stand for the right, left, and bilateral subtraction images, respectively. Features were then extracted from this set of images. Additionally, to study the feature-based asymmetry analysis, the average and absolute difference of each left-right pair of measurements was also analyzed, resulting in 860 additional features, resulting in a total of 2150 features per subject.

#### **2.7. Feature selection**

**2.5. Image enhancement**

results of the registering process.

116 New Perspectives in Breast Imaging

right breast, and right: color map of the subtraction image *I*

small and large intensity differences, respectively.

To study the appearance of the architectural distortions, two enhancing filters were applied to the images: a morphological high-frequency enhancement filter (*H*) designed to enhance fiber-like tissues, and a Laplacian of Gaussian filter (*L*) that enhances high-frequency patterns inside the breast tissue. Additionally, since the texture between normal and abnormal tissues is different [42], two texture maps were created. The first map computed the local standard deviation (*σ*) of the mammograms, and the second map computed the local fractal

**Figure 7.** Image subtraction example. Left: unaltered CC view of left breast, middle: horizontally flipped CC view of

**Figure 6.** Checkerboard comparison of images *pre* and *post* B-Spline registration. The image in the left shows a comparison between a left and horizontally flipped right breast before registration, and the right image shows the

∆. White and black pixels inside the breast tissue represent

The first step of the feature selection process consisted discarding highly correlated to avoid redundancy. For any pair of features with a Spearman correlation coefficient larger than 0.96, one feature was randomly selected to be kept, and the other removed from the selection. The dataset was normalized using the empirical distribution of the healthy subjects and a *z*-normalization was performed using the rank-based inverse normal transformation [44].

In order to select the most accurate and compact set of features from each dataset, the least absolute shrinkage and selection operator (LASSO) method was used [45]. The shrinkage and selection method minimizes the sum of squared errors and penalizes the regression coefficients, as described by Eq. (3) as follows:

$$\hat{\boldsymbol{\beta}}^{\text{lasso}} = \operatorname\*{argmin}\_{i=1}^{N} \sum\_{j=1}^{N} \left( y\_i - \boldsymbol{\beta}\_0 - \sum\_{j=1}^{p} \boldsymbol{x}\_{ij} \boldsymbol{\beta}\_j \right)^2 \\ \text{subject to:} \\ \sum\_{j=1}^{p} \left| \boldsymbol{\beta} \right| \le t \tag{3}$$

Given a set of input measurements *x*<sup>1</sup> …*xn* and an outcome y, the lasso method fits a linear model where *xi* is the covariate vector for the *i* th case and *yi* is the outcome, *t* is a tuning parameter that determines the amount of regularization, and *N* is the number of cases.

The multivariate search was performed using a class balanced data sample of 100 subjects for training and the remaining subjects as a blind test set. The models were calibrated using a leave-one-out cross-validation strategy, training the models at every split using *N* – 1 subjects and evaluating the model using the remaining subjects [46]. The final reported performance was obtained by applying the final model gathered on the training stage and evaluating it in the blind test set.

#### **3. Results**

A total of 1796 mammograms were successfully segmented. The image sets of nine subjects had to be removed from the experiment due to problems with the registration process, six were from MS, two from CS, and one from HS. All the remaining subjects were included in the subsequent stages of the analysis. The 2150 extracted features were filtered by the correlation process, removing 826 features.

**Table 1** shows the features that were selected for each model: the CS versus HS (*n* = 12), and the MS versus HS (*n* = 16). The former achieved an accuracy of 0.825 with an AUC of 0.882 and the latter an accuracy of 0.698 with an AUC of 0.807. **Figure 8** shows the ROC curves for both the models.


Note: Features are grouped by dataset, symmetric features are denoted with:

$$\begin{aligned} \text{Note:} & \text{Features are grouped by dataset, symmetric features are denoted with:}\\ I\_{\text{\\_array}} &= \frac{I\_r + I\_l}{2}, H\_{\text{\\_avgy}} = \frac{H\_r + H\_l}{2}, \mathcal{L}\_{\text{\\_avgy}} = \frac{L\_r + L\_l}{2}, \sigma\_{\text{\\_avgy}} = \frac{\sigma\_r + \sigma\_l}{2}, F\_{\text{\\_avgy}} = \frac{\sigma\_r + \sigma\_l}{2}, \\ I\_{\text{\\_A}} &= \left| I\_r - I\_l \right|, H\_{\text{\\_A}} = \left| H\_r - H\_l \right|, \mathcal{L}\_{\text{\\_A}} = \left| H\_r - H\_l \right|, \sigma\_{\text{\\_A}} = \left| H\_r - H\_l \right|, F\_{\text{\\_A}} = \left| H\_r - H\_l \right|. \end{aligned}$$

**Table 1.** Features of the proposed models.

**Figure 8.** ROC curves for the classification models. Left: MS versus HS, right: CS versus HS. The dashed line represents the model with the features from only the difference images, the solid line the one with features from only the raw images, and the dotted line the one with the features from all images (from Ref. [32]).
