Elasto-mammography: Theory, Algorithm, and Phantom Study

A new imaging modality framework, called elasto-mammography, is proposed to generate the elastograms of breast tissues based on conventional X-ray mammography. The displacement information is extracted from mammography projections before and after breast compression. Incorporating the displacement measurement, an elastography reconstruction algorithm is specifically developed to estimate the elastic moduli of heterogeneous breast tissues. Case studies with numerical breast phantoms are conducted to demonstrate the capability of the proposed elasto-mammography. Effects of noise with measurement, geometric mismatch, and elastic contrast ratio are evaluated in the numerical simulations. It is shown that the proposed methodology is stable and robust for characterization of the elastic moduli of breast tissues from the projective displacement measurement.


INTRODUCTION
Breast cancer is one of the major threats to public health in the world. Approximately 10% of women will develop breast cancer during the course of their lives in USA and Europe. The specific causes of breast cancer are yet unknown. Therefore early detection of breast tumor is the key to successful treatment.
X-ray mammography is the primary method for early detection of breast cancers [1]. According to the reports of US Food and Drug Administration, mammography can find 85 to 90 percent of breast cancers in women over 50, and can detect a lump up to two years before it can be sensed by manual palpation. While effective for detecting breast abnormality, mammography is not quite specific for differentiating benign and malignant masses, especially when the breast tissue is radiodense. A significant number of suspicious masses identified by mammography for surgical breast biopsy are in fact not malignant [2]. False-positive mammograms induce anxiety, distress, and intrusive thoughts.
It has been well recognized that the tissue stiffness plays an important role in diagnosis of breast cancers, as tumors are stiffer than the surrounding breast tissues [3,4], and malignant tumors are much stiffer than benign ones [5].
In other words, in vivo identification of the elastic moduli of normal and abnormal breast tissues, which describe the stiffness, should improve the accuracy of breast cancer diagnosis. There have been elastography studies based on either ultrasound or MRI breast imaging [6][7][8][9][10][11][12]. Ophir et al. [6,7] and Souchon et al. [8] proposed an ultrasound elastography modality for quantitative imaging of the elastic modulus distributions in biological tissues. Muthupillai et al. [9] and Manduca et al. [10] developed an algorithm to reconstruct the shear modulus distribution using acoustic strain wave propagation measured with MRI technique. Plewes et al. [11] and Samani et al. [12] provided a finite-element iteration method to reconstruct the distribution of elastic moduli in a breast containing suspicious tumors, based on the MRI deformation measurement under compression loading.
The objective of this study is to develop a new imaging modality, called elasto-mammography, for quantification of the elastic moduli of normal and cancerous breast tissues. In contrast to the previous breast elastography developments, elasto-mammography does not require additional biomedical imaging measurements and extra expense; that is, it combines the conventional low-dose X-ray mammography directly with our previously proposed tomography-based  elastography framework [13]. Specifically, by adopting certain anatomically well-motivated assumptions, the geometry of tumors is estimated from the mammography projections, as well as from the displacements at key points. The elastography reconstruction is further conducted with our highly efficient algorithm for the elastic parameters of tissues. This work is organized as follows. In Section 2, we recall our optimization-based algorithm for elastography reconstructions. We further present elasto-mammography simulations using numerical breast phantoms containing tumors in Section 3. Section 4 investigates the influences of various errors with the measurements, including noise with displacements, geometric mismatch, and elastic contrast ratios. Conclusions are drawn in the last section.

METHODOLOGY OF ELASTOGRAPHY RECONSTRUCTION
In this study, the mechanical properties of normal breast tissue and tumors are assumed to be linearly elastic and isotropic; that is, they are described with elastic Lamé parameters λ and μ. Typical clinical mammography applies two or M (M ≥ 2) individual compressions on the breast so that the maximum amount of tissues can be imaged and examined from different view angles. Information about the displacements is collected from projective images, and is used in the proposed elasto-mammography for identifying the Lamé parameters of the tissues, as will be described in the next section. The three-dimensional (3D) reconstruction algorithm for Lamé parameters is optimization based and follows our previously developed general framework [13]. The displacement and force quantities with the Ith experiment are denoted with superscript (I). Denoting for the Ith loading the measured displacement field in the biomedical medium of interest (Ω) as U (I) (x), and the calculated displacement field associated with the trial distribution of Lamé parameters (λ(x), μ(x)) as u (I) (x), the elasto-mammography seeks Lamé parameters such that the following objective functional Φ(λ(x), μ(x)) is optimally minimized: where the second-order tensor χ (I) simply takes diagonal matrix form, that is, (χ (I) (x)) i j = δ i j ω (I) i (x) (I = 1, 2, . . . , M; i, j = 1, 2, 3). The weight function ω (I) i (x) equals zero if the ith displacement component is not measured at point x. To include the surface displacement as measurement, ω (I) i (x) is considered as a generalized function on the boundary of Ω.
The elasto-mammography reconstruction follows an iterative optimization procedure, as schematically shown in Figure 1. We employ a large-scale limited-memory BFGS (L-BFGS) optimization method [14], which requires usersupplied gradients of the objective functional, that is, ∂Φ/∂λ and ∂Φ/∂μ. Continuum formulas for the gradients have been derived by Oberai et al. [15] for isotropic elastography and by Liu et al. [13] for general anisotropic cases. Here, we give the finite-element presentations for the objective functional Φ and its gradients.
Following standard finite-element procedures (e.g., [16]), the displacement field u (I) (x) is discretized as vector u (I) , which satisfies the equilibrium equation where vector F (I) represents the nodal force. Once finiteelement mesh is generated and discretization method is selected, the stiffness matrix K depends only on the Lamé parameters (λ(x), μ(x)). Consistently, the measured displacement field U (I) (x) is discretized as vector U (I) . Therefore, the objective functional (1) reads in which the matrix X (I) corresponds to the weight function χ (I) (x) and has the same dimension as the stiffness matrix K. It has been shown [13,15] that the gradients of Φ can be calculated conveniently via where the adjoint displacement w (I) is the solution of It is noted that u (I) and w (I) share the same Cholesky factorization (e.g., [16]) for the stiffness matrix K, thus the computational expense for solving w (I) (5) is minimal once u (I) is solved (2). In the proposed elasto-mammography technique, anatomic structures of the normal breast tissue and tumor are prescanned. Therefore the breast can be modeled as a piecewise homogenous medium, with uniform Lamé parameters (λ tissue , μ tissue ) for the normal breast tissue region and uniform parameters (λ tumor , μ tumor ) for the tumor region. Consequently, there are four gradients to be calculated: in which the inner summations are taken over all the elements in the tissue region for the gradients ∂Φ/∂λ tissue and ∂Φ/∂μ tissue , and over all the elements in the tumor region for the gradients ∂Φ/∂λ tumor and ∂Φ/∂μ tumor . In (6), (K e ) N is the element stiffness matrix of the Nth element, and (u (I) e ) N and (w (I) e ) N are the element nodal displacements with the Ith loading. It is also noted that ∂(K e ) N /∂λ and ∂(K e ) N /∂μ are constant matrices for the Nth element.

NUMERICAL SIMULATIONS
In this section, simulations are performed with numerical breast phantoms to identify the elastic parameters for normal tissue and tumor(s). The 3D breast phantoms contain one and two tumors, respectively. To simulate mammography compression, two types of loadings are applied, respectively, on the phantoms from different loading angles. Surface forces and part of the boundary displacements are extracted from the forward computation results, in compliance with the capability of projective imaging, and are used as input for the reconstruction. In the following text, the units are "cm" for length and displacements, "kPa" for elastic moduli, and "kN" for nodal forces.

Forward computations
Let us first consider a 3D phantom consisting of a halfspherical matrix with an embedded spherical inclusion (Figure 2(a)). The soft matrix, 10 cm in diameter and center at (x, y, z) = (0, 0, 0), imitates normal breast tissue. The hard inclusion, 1.5 cm in diameter and center at (2, 1.75, 2.25), simulates a tumor. The second phantom (Figure 2(b)) is similar, but has one more tumor of the same size and center at (−1.8, 0, 2). We denote these phantoms as "Phantom I" and "Phantom II," respectively. The phantoms are discretized with standard 3D tetrahedral elements. Phantom I consists of 1114 nodes and 6070 elements, while Phantom II consists of 1657 nodes and 9340 elements.   The materials are assumed isotropic. The Lamé parameters (λ, μ) are (25, 7.5) for the soft breast tissue and are (125, 25) for the tumor. Young's modulus E and Poisson's ratio ν are related to Lamé parameters via Hence, (E, ν) are (20.769, 0.384 62) for soft tissue and (70.833, 0.416 67) for tumor. Note that the tumor is assumed approximately 3.5 times as stiff as the surrounding tissue. In general, a tumor is much stiffer than the surrounding normal tissues. However, the ratio between the stiffness of cancerous and normal breast tissues found in the literature shows variations from a few times to a few ten times [4]. Skovoroda et al. [5] recognized that this is partially due to the nonlinearity effect in which the apparent stiffness increases with the strain applied. Effects of the contrast ratio on elastomammography will be discussed later.
In the simulations, the displacements are zero on the base surface where z = 0. Two compression loadings are applied on the upper surface of breast phantoms, respectively. For Loading 1, nodal force of 0.005 kN is applied on some of the surface nodes, as plotted in Figure 3 for Phantom II. Loading 2 applies nodal force |F x | = |F y | = 0.004 kN on the other set of surface nodes, as shown in Figure 4. Note that the loading  directions are different by π/4. For convenience, we denote the direction with Loading 1 as "direction 1," and that with Loading 2 as "direction 2."

Data acquisition
Given the Lamé parameters for normal breast tissue and tumor(s), the deformations in response to the external loadings 1 and 2 are obtained by solving finite element (2), respectively. First of all, the external surface of a breast at undeformed state can be reconstructed from images taken with a 3D camera (e.g., [17]). Then, for each external loading, two mammography projections are made in the compression direction; that is, one projection with undeformed state and one with deformed configuration. The shape and location of the tumor(s) can further be estimated from the undeformed projections along different orientations. It is recognized that real tumors may be irregular in shape and difficult to reconstruct accurately with limited number of projections. As a first-order approximation, we assume that tumors are spherical initially, and deform into ellipsoids. The initial size and center of tumors are readily estimated with two undeformed projections made in different directions. For instance, directions 1 and 2 in the present simulations, as plotted in Figures  5(a) and 6(a) for Phantom II. Note that Phantom I is considered as Phantom II with absence of the tumor initially at (−1.8, 0, 2).
We extract displacement information from projection of deformed configurations. Based on the micromechanics theory for deformation of an inclusion in a large medium (e.g., [18]), it is reasonable to estimate that an initially spherical tumor deforms into an ellipsoid. Because of the relatively simple uniaxial compression loadings applied in mammography, it is further approximated that vertexes of an object in an undeformed projection remain vertexes in the corresponding projection after compression deformation. For example, in Figure 5(b) for Loading 1, point A is the top vertex of tissue in undeformed projection. It moves to vertex A after deformation. Points B ∼ I are vertexes of the tumors in undeformed projection in direction 1. They displace to vertexes B ∼ I , respectively. Thus, by measuring the vertex locations in projections before and after deformation, their displacement information can be obtained. For example, displacement components u x and u z in Loading 1 for vertexes A ∼ I are extracted from the projections as in Figure 5. Acquisition of displacement information with Loading 2 makes use of projections, see Figures 6(a) and 6(b), and follows the same procedure. It is noted that the two tumors partly overlap in the projections in direction 2, and vertexes C and G are in shadow. For such a case, the vertex displacements are still attainable according to the grey density information in the projections with loss of some accuracy. The collected displacement data are denoted as U (1) and U (2) for elastomammography reconstruction.
Accurate displacement measurement with high spatial resolution will benefit elastography reconstruction in general. However, pinpoint tracking of large number of material points in an object is still a challenge in medical imaging [19], in particular for simple mammography projections that lack natural landmarks. Therefore, we propose elastomammography that only makes use of displacements of a few special points extracted directly from projections. As described above, the points include top vertex on the upper breast surface (A in Figure 5(b)) and vertexes of the tumors in projections (B ∼ I). Displacements measured at other points, for instance, on the external surface with a 3D camera, should enhance the efficiency and accuracy of elastomammography.

Ideal elasto-mammography
With the described data acquisition method, displacements at some key points are extracted from deformed and undeformed projections with the two compression loadings, and are used as measurements U (1) and U (2) for elastomammography reconstruction. Compression nodal forces applied on the surface are also known with the loadings. Given initial estimate, the Lamé parameters for tissue and tumor are reconstructed following our optimization procedure (Figure 1).   The ideal case is considered first; that is, the displacements, geometry, and compression nodal forces are exactly measured, and are used as input for reconstruction. Rows "Ideal I" and "Ideal II" of Table 1 give the reconstruction results for Phantom I and Phantom II, respectively. Convergent loci of the Lamé parameters (λ, μ) are plotted in Figure 7. The loci for Phantom I (Figure 7(a)) and Phantom II (Figure 7(b)) are very similar. It is observed that (λ, μ) of the tissue approach the real value rapidly. After about 20 iteration steps, their relative errors are well within the range of 5%. Then they experience some minor adjustment. In contrast, Lamé parameters of the tumor converge slower, in particular for λ, which starts to fall to the real value after about 40 steps. After about 50 steps, all parameters are accurately identified, with the largest error of about ±0.18% (for λ of tumor). Reconstructions using different initial estimates have been conducted. Very similar convergent profiles are found for the parameters, and highly accurate results are obtained. This indicates efficiency and uniqueness of the proposed elasto-mammography using projective measurements.
The slower convergent speed of Lamé parameters of the tumor, in particular for λ, is explained by the roles they play in the deformation due to the applied loadings, as discussed by Liu et al. [13]. In general, parameters with the most significant influence on the deformation are also those that are most accurately and easily identified. The influence of a parameter depends on size and location of the material region it belongs to, as well as characteristics of the deformation. For the present simulations, λ and μ of the tissue are dominant, while those of tumor are much less influential, due to the small size and deep location of the tumor(s). Slower convergence of λ for tumor indicates that the present loadings do not introduce enough volumetric strain in the tumor.

Effect of noise
The above elasto-mammography reconstructions are conducted using ideal inputs. In practice, several factors will affect the performance of elasto-mammography, the most common one among which is the noise with displacement measurement. To investigate the capability of the proposed elasto-mammography modality and algorithm to handle imperfect real data due to inevitable measurement errors, we conduct reconstruction using noisy input; that is, each component of U (1) and U (2) is added with a randomly selected relative error between −5% and 5%. Z. G. Wang et al.  The results are shown as "Noise I" (for Phantom I) and "Noise II" (for Phantom II) in Table 1, and the convergent loci are plotted in Figure 8. The overall convergent loci are very similar to the "ideal" cases. Lamé parameters (λ, μ) of the tissue need about 20 steps to approach closely to the real values, while those of the tumor need about 50 steps for convergence. The tissue parameters are very accurately identified, with the largest relative error of 4% for λ of Phantom II, and errors well within ±1% for the others. The Lamé parameters (λ, μ) of tumor, however, are not as robust, with relative errors of (−14.6%, 0.22%) and (17.4%, −15.5%) for Phantom I and Phantom II, respectively. In spite of these reconstruction errors, it is still positive that the elasto-mammography results are accurate enough for diagnosis of tumors, noting the significant differences of stiffness between normal tissue, and benign and malignant tumors (e.g., [3][4][5]). The better robustness of the tissue parameters is also explained by the strong roles they play in the deformation, as we have discussed above. Furthermore, as suggested by Liu et al. [13], multiple sets of well-designed loadings should help to bring out the influences of all the material parameters, and thus suppress the effects of noise.

Effect of geometry mismatch
Another concern for elasto-mammography is the geometric depiction of the tumor. As described in the section of data acquisition, we use a simple sphere to approximate a real tumor, and estimate its size and location from two undeformed mammography projections. This inevitably introduces geometric mismatch for practical elasto-mammography. To investigate the effect of geometry mismatch, the two phantoms are redesigned by replacing the spherical tumors with cubic tumors. Note that the edge length of the cube is 3/ √ 5 cm. Forward simulations are conducted under the same Loading 1 and Loading 2 with the new phantoms. Then, mammography projections are made of the new undeformed and deformed configurations. To extract geometric and displacement data from the projections, we still use spherical approach. As schematically shown in Figure 9, a cubic tumor is approximated with a spherical one, whose size and location are determined by the two undeformed projections in direction 1 and direction 2. Then, the estimated spherical tumors are used for elasto-mammography reconstruction of the material parameters. The results are shown as "Noise I" (for Phantom I) and "Noise II" (for Phantom II) in Table 1. Convergent loci are found to be similar to the previous cases, and are not shown.
The tissue parameters again show excellent robustness. The geometric mismatch introduces relative errors less than 0.17%. Due to the relatively small size of the tumor(s), their Lamé parameters (λ, μ) are more sensitive to geometric mismatch, with relative errors (3.52%, 0.40%) for Phantom I and (24.2%, 7.78%) for Phantom II. In comparison to the displacement noise, the geometry mismatch seems to have slightly less overall influence on the reconstruction results. However, this point is based on the current phantoms, and needs further investigation.  In this study, a perfect sphere is assumed to simulate the real tumor. From this simplified model, the information of deformation can be easily obtained from the projections. However, investigation of geometry mismatch is needed since most of real tumors have irregular shapes. A cube (a rectangle or square in projections) is used to represent a real tumor, and a sphere (a circle in projections) approximates it. From Figure 9(a), it can be seen that most of the areas between circle and square overlap. The results demonstrate that geometry mismatch does not have a great influence on this reconstruction. It is therefore suggested that, for an irregular shape of a real tumor in projections, we could choose a circle to approximate it and do reconstruction based on this simplified model.

Effect of contrast ratio
The material stiffness is a key feature that distinguishes benign from malignant tumors [4][5][6]. Contrast ratio, defined as the ratio between Young's modulus of tumor and normal breast tissue, covers a wide range. For benign tumors, contrast ratio typically varies from 2.0 to about 5.0. For malignant tumors, it is considerably higher. Our numerical experiments [13] indicate that accuracy of elastography reconstruction depends not only on type of loading and measurement accuracy, but also on the contrast ratio. In case that the tumor is very hard, the material parameters may be identified qualitatively, but not quantitatively.
To investigate the effect of contrast ratio, we conducted elasto-mammography reconstructions with soft and hard phantoms, whose Lamé parameters (λ, μ) are set to create contrast ratios (CR) of 1.5 and 8.0, respectively. Table 2 gives the real Lamé parameters and reconstruction results. Compared to the previous case with CR about 3.5 ( Table 1, "Ideal I" and "Ideal II"), results for the soft phantoms (CR = 1.5) are even more accurate, in particular for λ of the tumor. For the hard phantoms (CR = 8.0), the tissue parameters are also exact; however, λ of tumor carries relative reconstruction errors of about 13% for both phantoms, which is considerably larger than the soft cases with CR = 3.5 and 1.5. The reason is that deformation of a relatively softer tumor is larger than a hard one, and thus is more sensitive to small variation of its material parameters. Also as discussed above, λ of the tumor seems to have the least influence on the specific deformations considered in the simulations. On the other hand, it is convincing that the proposed elastomammography is efficient in revealing the contrast ratio and telling whether a tumor is malignant or benign. In general, a tumor is suspected of malignancy when the contrast ratio is higher than 6. In the present simulations, when the "real" contrast ratio is 8.0, the elasto-mammography reconstruction yields 8.11, which is fully acceptable for the diagnostic purpose.
The elastography simulations of Liu et al. [13] assumed tumor and normal breast tissue as general anisotropic materials, and applied four sets of loading on a breast phantom to bring out all the elastic parameters. Their isotropic simulation suggested that displacements measured from a single loading are adequate for unique identification of the Lamé parameters (λ, μ) of tissue and tumor. However, their measurement includes displacement on the entire external surface and the tissue-tumor interface of a breast, requiring more complex imaging equipment. In our elastomammography proposal, displacement measurement has been reduced to a few vertexes, and can be readily obtained from simple mammography projections. A tradeoff is that two or more sets of compression loadings may be needed to obtain adequate identification.
Mathematical proof for uniqueness results of elastomammography using projection measurements is yet under further investigation. Our simulations always yield the same material parameters (within the numerical processing errors), regardless of the initial estimate. With ideal measurements, the resulting parameters exactly match the real values specified for the models. When displacement noise and geometry mismatch are taken into consideration, the resulting parameters have reconstruction errors, however, are close enough to their real values for application purpose. In summary, the proposed elasto-mammography method is numerically stable and robust, is relatively simple to perform, and thus has great potential for clinical applications.

CONCLUSIONS
A new method that combines elastography and mammography to reconstruct the elastic field of the breast is reported. Displacement and geometry measured from deformed and undeformed mammography projections are applied as input data to reconstruct the isotropic material parameters for normal breast tissue and tumor. Our numerical simulations demonstrate that unique and accurate results can be obtained using information extracted from only two sets of projections. Displacement noise, geometry mismatch, and material contrast ratio do not adversely affect the results, demonstrating that our method is stable and robust. These findings are sufficiently encouraging to warrant both further development and clinical evaluation of our reconstruction method.