# Elastography Method for Reconstruction of Nonlinear Breast Tissue Properties

Authors: Z. G. Wang; Y. Liu; G. Wang; L. Z. Sun
Journal: International Journal of Biomedical Imaging
 (2009)

---

## Abstract

Elastography is developed as a quantitative approach to imaging linear elastic properties of tissues to detect suspicious tumors. In this paper a nonlinear elastography method is introduced for reconstruction of complex breast tissue properties. The elastic parameters are estimated by optimally minimizing the difference between the computed forces and experimental measures. A nonlinear adjoint method is derived to calculate the gradient of the objective function, which significantly enhances the numerical efficiency and stability. Simulations are conducted on a three-dimensional heterogeneous breast phantom extracting from real imaging including fatty tissue, glandular tissue, and tumors. An
exponential-form of nonlinear material model is applied. The effect of noise is taken into account.
Results demonstrate that the proposed nonlinear method opens the door toward nonlinear elastography
and provides guidelines for future development and clinical application in breast cancer study.


---

## Body

## 1. Introduction

Breast cancer is one of the  major threats to public health all over the world. Currently, X-ray mammography is the primary method for early detection and characterization of breast tumors [1]. While more efficient in detecting malignancies as age increases or the breast becomes fatty, mammography fails to detect small cancers in dense breasts. Further, mammography may not be specific in terms of tumor benignity and malignancy [2–4].To solve these problems associated with mammography, a number of technologies have been explored. Detection and characterization of breast tumors can be enhanced by recognizing the difference of elastic modulus (stiffness) among normal soft tissues and malignant and benign tumors. Elastic properties of breast tissues may become an indicator of histological diagnosis [5]. An imaging technology called elastography was developed as an approach to imaging tissue elastic modulus in a quantitative manner for detection of breast tumors in 1990s. The general basis of elastography is to induce motion within tissues under investigation by either an external or internal mechanical stimulation. Conventional medical imaging modalities are then used to measure the spatial deformation, from which the mechanical properties can subsequently be reconstructed. Most simulations are either ultrasound or magnetic resonance (MR) image-based [6–10] that take the dynamic or quasistatic interior displacement field, completely or partially, as input for identification of the elasticity properties. Current elastography reconstruction framework is based on the assumption of linear elasticity theory. It is shown, however, that the deformation of most biological soft tissues are not linearly elastic [5, 11]. Consideration of nonlinear models is essential for elastography in clinical applications.The present work aims at developing a nonlinear elastography model for breast tissues. We introduce a nonlinear adjoint gradient method that significantly improves the numerical efficiency and enhances the stability of elastography reconstructions. In Section2, we present an algorithm based on nonlinear adjoint gradient method. We apply the iterative Newton method to solve unknown displacements and forces. In order to find elastic modulus distribution that minimizes the objective function based on measured and calculated forces, the adjoint gradient method is employed to provide user-supplied gradients in nonlinear elastography. Numerical simulation is described in Section 3 including establishing a three-dimensional (3D) heterogeneous FEM breast phantom and applying exponential-form nonlinear material model. Four types of compressive loadings are applied in the forward problem. The iterative reconstruction based on force measurements on surface is also detailed. In Section 4, the result of the inverse problem is obtained and the effect of noise is investigated. A user-defined penalty function is introduced to reduce the impact of noise on reconstruction.


## 2. Methodology

### 2.1. Finite-Strain Deformation Equations

We use continuum description for the breast tissues.  LetΩ0 be a biological object. From the displacement tensor u(X) based on Lagragian coordinate system X, the deformation gradient is F=I+∂u/∂X and the Green strain is E=(FT·F-I)/2 where “·” denotes the contraction operation between two tensors.  The breast tissues are assumed hyperelastic so that the second Piola-Kirchhoff stress is S=∂W(E;p)/∂E, where W is the strain energy and p denotes material parameters in the model. We use “;” to separate material parameters from deformation variables. The governing equation and boundary conditions for u are
(1)∇·(S·FT)+b(X,u(X))=0X∈Ω0,N(X)·(S·FT)=t(X,u(X))X∈Γt0,u(X)=u̅(X)X∈Γu0.
The boundary of Ω0 consists of Γt0, with N the outer normal, on which external force t is applied, and Γu0 where displacement u̅ is prescribed. Here, we consider general problems that the body force b and surface force t are deformation dependent.  Following the standard finite-element (FE) method [12], the displacement u is discretized as nodal displacement vector {u}={u1,u2}T, where u2 corresponds to u̅ prescribed on Γu0, and u1 is to be solved from nonlinear equations:
(2){f1in(u1,u2;p)f2in(u1,u2;p)}-{f1out(u1,u2)f2out}={00}.
The internal nodal force fin corresponds to stress S; that is, it changes with u1 and material parameters  p, as u2 is given. The external nodal force f1out is due to the prescribed surface force t on Γt0 and body force b in Ω0. It varies with the displacement in large deformation. The nodal force f2out is the unknown constraint force on Γu0.  A classic Quasi-Newton method is employed to solve (2) for u1 as shown in Appendix A.


### 2.2. Nonlinear Elastography Algorithm

Experimental measurements for elastography include displacement and force. We consider that the biological objectΩ0 is discretized into FE mesh, and the measurements are discretized consistently into nodal displacement and nodal force. We catalog the measurements as the following. (i) If the displacement at a node is known (prescribed or measured after deformation), it will be included into u2 which is considered “prescribed” in FE (2). The corresponding nodal force (known or unknown) will be in the constraint force, denoted as F2M. Note that F2M corresponds to f2out in (2).  (ii) All the other nodal displacements will be in u1 and the corresponding nodal force will be in f1out. For category (ii), f1out must be considered “prescribed” to fulfill the requirement of the well-postness of a solid-mechanics problem.For elastography problem, the obtainedu2 and f1out are known in (2), and the constraint force F2M is considered as “measurement.” For given material parameters p, the unknown displacement u1 and constraint force f2out (which depends on p) will be solved from FE (2). Elastography procedure thus searches for p so that the overall difference between measured F2M and computed f2out is minimum, that is, minimize objective function
(3)Φ(P)=(f2out-F2M)TΛ(f2out-F2M),
where the diagonal weight matrix Λ=diag(a1,a2,…,aj,…). The component aj=1 when the jth component of F2M is measured or aj=0 otherwise. The present algorithm is mathematically equivalent to our previous work [10] where the displacement is used as “measurement.” For breast tissue whose tangent stiffness significantly increases with strain, this “force version” shows better numerical efficiency.Efficient and robust optimization-based elastography schemes request user-supplied gradient∂Φ/∂p.  The previous elastography studies used direct or approximate finite-difference method [13, 14] for the gradient calculation.  The computational expense of these methods increases proportionally with the number of material parameters, and becomes unaffordable for problems involving finite-strain deformation. Recently, an adjoint method was introduced to compute the gradient analytically [10, 15–17]. The corresponding nonlinear finite element formulas are shown in Appendix B.  After u1 is solved for given p, the gradient is readily obtained as
(4)∂Φ∂p={w1w2}T{∂f1in/∂p∂f2in/∂p},
where the virtual adjoint displacements w1 and w2 are solved from linear equations
(5)w2=2Λ(f2out-F2M),(Keff)Tw1=-(∂f2in∂u1)Tw2,
with the tangent stiffness matrix Keff defined in (A2). The most significant features of the adjoint method are its analytical formulation, high accuracy, and computational efficiency [18]. Since Keff and its LU factorization have been calculated when solving the FE (2), as in (A1), the additional expense for w1 and w2 in (5) is minimal. Furthermore, it only needs to solve one linear (5) regardless of the number of unknown parameters in p.As shown in Figure1, a well-tested L-BFGS subroutine is applied for the present optimization-based nonlinear elastography.Figure 1
Overall flowchart for nonlinear reconstruction of material parametersλ,μ, and γ of breast tissues.


## 2.1. Finite-Strain Deformation Equations

We use continuum description for the breast tissues.  LetΩ0 be a biological object. From the displacement tensor u(X) based on Lagragian coordinate system X, the deformation gradient is F=I+∂u/∂X and the Green strain is E=(FT·F-I)/2 where “·” denotes the contraction operation between two tensors.  The breast tissues are assumed hyperelastic so that the second Piola-Kirchhoff stress is S=∂W(E;p)/∂E, where W is the strain energy and p denotes material parameters in the model. We use “;” to separate material parameters from deformation variables. The governing equation and boundary conditions for u are
(1)∇·(S·FT)+b(X,u(X))=0X∈Ω0,N(X)·(S·FT)=t(X,u(X))X∈Γt0,u(X)=u̅(X)X∈Γu0.
The boundary of Ω0 consists of Γt0, with N the outer normal, on which external force t is applied, and Γu0 where displacement u̅ is prescribed. Here, we consider general problems that the body force b and surface force t are deformation dependent.  Following the standard finite-element (FE) method [12], the displacement u is discretized as nodal displacement vector {u}={u1,u2}T, where u2 corresponds to u̅ prescribed on Γu0, and u1 is to be solved from nonlinear equations:
(2){f1in(u1,u2;p)f2in(u1,u2;p)}-{f1out(u1,u2)f2out}={00}.
The internal nodal force fin corresponds to stress S; that is, it changes with u1 and material parameters  p, as u2 is given. The external nodal force f1out is due to the prescribed surface force t on Γt0 and body force b in Ω0. It varies with the displacement in large deformation. The nodal force f2out is the unknown constraint force on Γu0.  A classic Quasi-Newton method is employed to solve (2) for u1 as shown in Appendix A.


## 2.2. Nonlinear Elastography Algorithm

Experimental measurements for elastography include displacement and force. We consider that the biological objectΩ0 is discretized into FE mesh, and the measurements are discretized consistently into nodal displacement and nodal force. We catalog the measurements as the following. (i) If the displacement at a node is known (prescribed or measured after deformation), it will be included into u2 which is considered “prescribed” in FE (2). The corresponding nodal force (known or unknown) will be in the constraint force, denoted as F2M. Note that F2M corresponds to f2out in (2).  (ii) All the other nodal displacements will be in u1 and the corresponding nodal force will be in f1out. For category (ii), f1out must be considered “prescribed” to fulfill the requirement of the well-postness of a solid-mechanics problem.For elastography problem, the obtainedu2 and f1out are known in (2), and the constraint force F2M is considered as “measurement.” For given material parameters p, the unknown displacement u1 and constraint force f2out (which depends on p) will be solved from FE (2). Elastography procedure thus searches for p so that the overall difference between measured F2M and computed f2out is minimum, that is, minimize objective function
(3)Φ(P)=(f2out-F2M)TΛ(f2out-F2M),
where the diagonal weight matrix Λ=diag(a1,a2,…,aj,…). The component aj=1 when the jth component of F2M is measured or aj=0 otherwise. The present algorithm is mathematically equivalent to our previous work [10] where the displacement is used as “measurement.” For breast tissue whose tangent stiffness significantly increases with strain, this “force version” shows better numerical efficiency.Efficient and robust optimization-based elastography schemes request user-supplied gradient∂Φ/∂p.  The previous elastography studies used direct or approximate finite-difference method [13, 14] for the gradient calculation.  The computational expense of these methods increases proportionally with the number of material parameters, and becomes unaffordable for problems involving finite-strain deformation. Recently, an adjoint method was introduced to compute the gradient analytically [10, 15–17]. The corresponding nonlinear finite element formulas are shown in Appendix B.  After u1 is solved for given p, the gradient is readily obtained as
(4)∂Φ∂p={w1w2}T{∂f1in/∂p∂f2in/∂p},
where the virtual adjoint displacements w1 and w2 are solved from linear equations
(5)w2=2Λ(f2out-F2M),(Keff)Tw1=-(∂f2in∂u1)Tw2,
with the tangent stiffness matrix Keff defined in (A2). The most significant features of the adjoint method are its analytical formulation, high accuracy, and computational efficiency [18]. Since Keff and its LU factorization have been calculated when solving the FE (2), as in (A1), the additional expense for w1 and w2 in (5) is minimal. Furthermore, it only needs to solve one linear (5) regardless of the number of unknown parameters in p.As shown in Figure1, a well-tested L-BFGS subroutine is applied for the present optimization-based nonlinear elastography.Figure 1
Overall flowchart for nonlinear reconstruction of material parametersλ,μ, and γ of breast tissues.


## 3. Numerical Simulations

This section presents phantom simulations to identify the nonlinear elastic properties of the fatty, glandular, and cancerous tissues in a breast. First of all, a 3D breast FEM phantom attracting from the real data is introduced. Fung’s model [19] is applied to describe the deformation of breast. With the applied loading, the forward-problem computation is performed. Furthermore, boundary forces are extracted from the forward computation results and are used as input for reconstruction for a breast phantom.


### 3.1. Breast Phantom

To perform numerical simulations, a 3D numerical heterogeneous breast phantom extracting from real CT images, containing fatty tissue, glandular tissue, and a tumor is established (Figure2). Boundaries of these regions are described with sets of splines. The phantom is discretized with standard 3D tetrahedral elements, consisting of 7303 elements and 1583 nodes (Figure 3).Figure 2
A 3D heterogeneous breast phantom extracted from real CT images, containing fatty tissue, glandular tissue, and a tumor.Figure 3
The phantom is discretized with standard 3D tetrahedral elements, consisting of 7303 elements and 1583 nodes.It is well known that the mechanical behavior of biological soft tissue is nonlinear. Hyperelastic models have commonly been used to represent the stress-strain relation of biological soft tissue [20, 21]. Fung and coworkers developed a set of hyperelastic models for bio-tissues [20]. In this work, we employed a Fung-type isotropic strain-energy function for the breast tissues, as W(E;p)=γ(exp(λ(E:I)2)/2+μE:E)-1 , where Eis the Green strain, p=={γ,λ,μ} denotes the material parameters. The second Piola-Kirchhoff stress Sis thus(6)S=∂W(E;p)∂E,
Specifically, in uniaxial tension/compression, the relation between the axial strain Eandaxial stressSis(7)S=γexp(EYoungE22)EYoungE,
where EYoung=μ(3λ+2μ)(λ+μ)-1  is the Young’s modulus.Equation (7) is the exponentially nonlinear model that can be applied for breast tissues. Two parameters γ and EYoung can be determined by the experiment [22]. Wellman [5] developed a technique to measure the nonlinear elastic parameters of breast tissues using force-displacement data of thin slice tissues undergoing indentation experiment. The tissue samples tested were obtained during surgery and were tested immediately after removal from the body. Six breast tissues (fatty, gland, lobular carcinoma, fibroadenoma, infiltrating ductal carcinoma, and ductal carcinoma in situ) were tested with their nonlinear stress-train curves shown in Figure 4. It is shown that the mechanical properties between normal soft breast tissue and tumors are quite different. Tumors are stiffer than the surrounding breast tissues and malignant tumors are much stiffer than benign ones.Figure 4
Nonlinear stress-strain curves for six breast tissues (fatty tissue, glandular tissue, lobular carcinoma, fibroadenoma, infiltrating ductal carcinoma, and ductal carcinoma in situ). (Figure redrawn based on [5] and (7).In this study, exponential hyperelastic model is used and material parameters are regressed from Wellman [5] experiment. Three materials are selected: fatty tissue, glandular tissue and ductal carcinoma in situ. Elastic parameters are λf=35, μf=12.5, γf=0.4 for fatty tissue, λg=50, μg=25, γg=0.25 for glandular tissue and λd=80, μd=35, γd=1.5 for ductal carcinoma in situ.  λ and μ are dimensionless, γ is in kPa.


### 3.2. Forward Problem

After the breast phantom and material model are established, a forward problem is solved in which material parameters and external loading are given and deformation is solved. Displacementsu=0 are prescribed on the base of phantom. Tilted compression by paddles is applied on upper surface (Figure 5). The paddle close to tumor gives compression and another paddle is fixed during loading. Four types of compression loading with different angles are applied. As shown in Figure 6, blue lines represent paddle locations before loadings, green lines represent those after loadings. Figure 7 shows the comparison of paddle locations in four loadings. Note that the right paddle is fixed for all loadings.Figure 5
Titled compression is given by two paddles. The one close to tumor gives displacement loading on breast phantom and another is fixed during loading.Four types of compression loading by paddles. Blue lines represent paddle locations before loadings, green lines represent after loadings.  Note that the right paddle is fixed for all loadings.
(a)(b)(c)(d)Comparison of paddle locations in four loadings. Note that the right paddle is fixed for all loadings.
(a)(b)The reason that four sets of loadings are applied in forward problem is to provide more information to reconstruct material parameters. Most of inverse problem in elasticity is nonuniqueness. Previous research [10] has demonstrated that one set of measurements of displacements and forces may not provide sufficient information of the reconstruction of modulus distribution.Given material parameters and loadings, the displacements and forces can be calculated. The surface force will be used as input to reconstruction material parameters in the inverse problem. In fact, surface displacement and force are equivalent as input to solve inverse problem. Most of previously research applied displacement in linear elastography. In this nonlinear elastography study, however, it is found the reconstruction is more sensitive to force than displacement. Therefore surface force is measured and compared with calculated force to reconstruct material parameters.


### 3.3. Inverse Problem

Reconstruction for nonlinear elastic moduli in 3D breast phantom take input extracted from the deformation in response to loading modes A~D described in Figure 6. In each loading, the forces on surface nodes are measured as input. Then initial guesses for material parameters are given. In this study, the same initial guesses are applied to three materials: λ0=20, μ0=10, and γ0=1. Therefore, the surface force can be calculated based on initial guesses of material parameters. The error between calculated and measured force is used to set up the objective function (3).  Following the iterative optimization procedures (Figure 1), gradients of the objective function are calculated and material parameters are updated to approach the optimal values.


## 3.1. Breast Phantom

To perform numerical simulations, a 3D numerical heterogeneous breast phantom extracting from real CT images, containing fatty tissue, glandular tissue, and a tumor is established (Figure2). Boundaries of these regions are described with sets of splines. The phantom is discretized with standard 3D tetrahedral elements, consisting of 7303 elements and 1583 nodes (Figure 3).Figure 2
A 3D heterogeneous breast phantom extracted from real CT images, containing fatty tissue, glandular tissue, and a tumor.Figure 3
The phantom is discretized with standard 3D tetrahedral elements, consisting of 7303 elements and 1583 nodes.It is well known that the mechanical behavior of biological soft tissue is nonlinear. Hyperelastic models have commonly been used to represent the stress-strain relation of biological soft tissue [20, 21]. Fung and coworkers developed a set of hyperelastic models for bio-tissues [20]. In this work, we employed a Fung-type isotropic strain-energy function for the breast tissues, as W(E;p)=γ(exp(λ(E:I)2)/2+μE:E)-1 , where Eis the Green strain, p=={γ,λ,μ} denotes the material parameters. The second Piola-Kirchhoff stress Sis thus(6)S=∂W(E;p)∂E,
Specifically, in uniaxial tension/compression, the relation between the axial strain Eandaxial stressSis(7)S=γexp(EYoungE22)EYoungE,
where EYoung=μ(3λ+2μ)(λ+μ)-1  is the Young’s modulus.Equation (7) is the exponentially nonlinear model that can be applied for breast tissues. Two parameters γ and EYoung can be determined by the experiment [22]. Wellman [5] developed a technique to measure the nonlinear elastic parameters of breast tissues using force-displacement data of thin slice tissues undergoing indentation experiment. The tissue samples tested were obtained during surgery and were tested immediately after removal from the body. Six breast tissues (fatty, gland, lobular carcinoma, fibroadenoma, infiltrating ductal carcinoma, and ductal carcinoma in situ) were tested with their nonlinear stress-train curves shown in Figure 4. It is shown that the mechanical properties between normal soft breast tissue and tumors are quite different. Tumors are stiffer than the surrounding breast tissues and malignant tumors are much stiffer than benign ones.Figure 4
Nonlinear stress-strain curves for six breast tissues (fatty tissue, glandular tissue, lobular carcinoma, fibroadenoma, infiltrating ductal carcinoma, and ductal carcinoma in situ). (Figure redrawn based on [5] and (7).In this study, exponential hyperelastic model is used and material parameters are regressed from Wellman [5] experiment. Three materials are selected: fatty tissue, glandular tissue and ductal carcinoma in situ. Elastic parameters are λf=35, μf=12.5, γf=0.4 for fatty tissue, λg=50, μg=25, γg=0.25 for glandular tissue and λd=80, μd=35, γd=1.5 for ductal carcinoma in situ.  λ and μ are dimensionless, γ is in kPa.


## 3.2. Forward Problem

After the breast phantom and material model are established, a forward problem is solved in which material parameters and external loading are given and deformation is solved. Displacementsu=0 are prescribed on the base of phantom. Tilted compression by paddles is applied on upper surface (Figure 5). The paddle close to tumor gives compression and another paddle is fixed during loading. Four types of compression loading with different angles are applied. As shown in Figure 6, blue lines represent paddle locations before loadings, green lines represent those after loadings. Figure 7 shows the comparison of paddle locations in four loadings. Note that the right paddle is fixed for all loadings.Figure 5
Titled compression is given by two paddles. The one close to tumor gives displacement loading on breast phantom and another is fixed during loading.Four types of compression loading by paddles. Blue lines represent paddle locations before loadings, green lines represent after loadings.  Note that the right paddle is fixed for all loadings.
(a)(b)(c)(d)Comparison of paddle locations in four loadings. Note that the right paddle is fixed for all loadings.
(a)(b)The reason that four sets of loadings are applied in forward problem is to provide more information to reconstruct material parameters. Most of inverse problem in elasticity is nonuniqueness. Previous research [10] has demonstrated that one set of measurements of displacements and forces may not provide sufficient information of the reconstruction of modulus distribution.Given material parameters and loadings, the displacements and forces can be calculated. The surface force will be used as input to reconstruction material parameters in the inverse problem. In fact, surface displacement and force are equivalent as input to solve inverse problem. Most of previously research applied displacement in linear elastography. In this nonlinear elastography study, however, it is found the reconstruction is more sensitive to force than displacement. Therefore surface force is measured and compared with calculated force to reconstruct material parameters.


## 3.3. Inverse Problem

Reconstruction for nonlinear elastic moduli in 3D breast phantom take input extracted from the deformation in response to loading modes A~D described in Figure 6. In each loading, the forces on surface nodes are measured as input. Then initial guesses for material parameters are given. In this study, the same initial guesses are applied to three materials: λ0=20, μ0=10, and γ0=1. Therefore, the surface force can be calculated based on initial guesses of material parameters. The error between calculated and measured force is used to set up the objective function (3).  Following the iterative optimization procedures (Figure 1), gradients of the objective function are calculated and material parameters are updated to approach the optimal values.


## 4. Results and Discussion

### 4.1. Ideal Input

For three materials, elastic parameters are:λf=35, μf=12.5, and γf=0.4 for fatty tissue; λg=50, μg=25, and γg=0.25 for glandular tissue, and λd=80, μd=35, and γd=1.5 for ductal carcinoma in situ. Table 1 shows the initial estimate and reconstructed results for nonlinear elastography. The results in the first part are based on the ideal input. It is demonstrated that the reconstructed results are very close to the real values. The maximum error is 1.916% (λ for tumor) since the effect of the tumor on surface force measurement is the smallest. A similar situation occurs in linear elastography reconstruction.Table 1
Initial estimate and reconstructed results for nonlinear elastogrpahy. The results in the first part are based on ideal input (without noise). The ones in the second part are based on input with 5% noise and regularization is not used. The third part is based on input with 5% noise and regularization is applied to reduce the impact of noise.  (λ and μ are dimensionless, γ is in kPa.)
FattyGlandularTumorλfμfγfλgμgγgλdμdγdReal3512.50.450250.2580351.5Guess201012010120101Ideal inputRecon34.998812.50040.399950.051224.99980.249881.533134.94231.50035% Noise, without regularizationRecon22.27538.41471.503856.197021.54140.25090.000141.55622.00105% Noise, with exponential form regularizationRecon37.5629012.483180.459250.001425.00420.244480.000239.15391.5020


### 4.2. Adjoint Methods

Elastography includes forward and inverse problem. In forward problem, material parameters and loadings are given to calculate the deformation; while in inverse problem, external loadings and deformation are known to reconstruct material parameters. Most researchers established certain objective function and minimized it with a proposed iterative algorithm. The challenge is how to calculate the gradient of objective function efficiently and accurately. A straightforward calculation of gradients requires solving stiffness matrix in each iteration, which takes most of the time consumed in the finite element method.In this study an adjoint method is employed to analytically calculate the gradients. Oberai et al. [16] adopted an adjoint method and proposed a numerical scheme for reconstructing the nonuniform shear modulus field for incompressible isotropic materials using one component of displacement field. Liu et al. [10] applied this method for anisotropic materials. In this study, the adjoint method is applied for nonlinear elastography. The advantage of adjoint method is to solve two adjoint displacements during each iteration (w1 and w2 at (A4) in Appendix B), instead of the whole stiffness matrix, that increases the numerical efficiency significantly. In [18] Oberai et al. compared three different iterative methods: (1) a gradient-based method where the adjoint approach is used to calculate the gradient; (2) a gradient-based method where the straightforward approach is used to calculate the gradient, and (3) the Gauss-Newton method. The results show that “leading-order costs for the gradient-based method with the adjoint approach are smaller than the other two methods.” In fact, without the adjoint method, nonlinear elastography can only be discussed in concepts [23] or applied on simple objects using supercomputing power [24]. While in this study, material parameters in a real breast FEM phantom are reconstructed accurately by adopting this adjoint method. The results are encouraging for further clinical applications.


### 4.3. Multiple Sets of Measurements

Above results are based on four sets of force measurements on surface. For 2D isotropic elastography, Barbone and Bamber [25] has shown that one set of measurements for the displacements and forces, especially those taken only on the boundaries, may not provide sufficient information to the reconstruction of modulus distribution, because of the nonuniqueness nature of most inverse problems in elasticity. Barbone and Gokhale [26] further demonstrated the feasibility of using multiple displacement fields to reduce the likelihood of nonuniqueness for 2D isotropic elastography. Liu et al. [10] discussed the multiple sets of measurements in 3D anisotropic media. However research effort for 3D nonlinear elastography has not been found in literatures. In this study four independent tilted compression loadings are designed and surface force is measured to reconstruct internal material parameters. Due to different initial guesses, stable material parameter can be reconstructed, showing that the multiple sets of measurements are feasible to obtain the unique solution of inverse problem in nonlinear elastography.The key points for using multiple sets of measurements are to bring more deformation modes simultaneously into consideration. The loadings should be close to tumors in order to make tumors have larger deformation. In Figures5–7, four sets of loadings close to tumors are designed to obtain more deformation for reconstruction. This explanation serves as a guideline for design of loadings in clinical applications. Reduction of the number of required loadings will increase the clinical efficiency and benefit the patients. In other words, loadings with the richest stress-strain information are most desired. The design of feasible loadings is important for success in clinical application of nonlinear elastography.


### 4.4. Input with Noise

The above results are based on ideal input. However, noise cannot be avoidable in experiments. Its impact on reconstruction is therefore investigated with 5% of noise applied on surface measurements.  The results are shown in the second part of Table1. It is obvious that the reconstructed parameters are away from the real values. A possible reason is that, while the goal is to find real material parameters that minimize the objective function, a global minimum is not well defined with noise. In another words, real material parameters may not give a global minimum of the objective function with noise. Several similar parameters around real ones may give some local minimums. Our algorithm may reach one of them but fails to find real one because they all give local minimums. A typical way to solve this problem is to add a penalty term which provides additional constraint on the solution space. The penalty term may push the solution into the right area of the solution space and minimize the resulting objective.  A new objective function is therefore proposed as
(8)F=Φ+χΠ,
where Φ is the original objective function (3) which represents the difference between measured and calculated values. In addition, χ is the regularization factor and Π is the penalty term.Specific forms of penalty term have been designed for different problems [27]. In this study, an exponential form of penalty term is applied as
(9)Π=∑k=1K(1-exp(-(ξk-ak)2)),
where K=9 is the total number of material parameters, and ξk and ak are the reconstructed and true elastic parameters, respectively. If the true material parameters are unknown, ak can be estimated as close as possible. This is reasonable since several experiments have provided the nonlinear elastic parameter for breast tissues.The results are shown in the third part of Table1, demonstrating that the regularization method significantly improves the reconstruction data. Most of parameters are close enough to the true values. The largest error occurs for χ in fatty tissue with approximately 15%. By adding the penalty term, the reconstructed values are pulled from local minimum into the range close to true values. It is noted that the regularization factor χ varies in different scenarios. For example, χ can be smaller for higher confidence of experiment-data accuracy. On the contrary, if ak is closer to the true value, χ should be larger. For this study, the elastic parameters in fatty and glandular tissues are stable, comparing with the ones in tumor (Figure 4). Therefore, larger regularization factors are used for fatty and glandular tissues while smaller χ can be applied for tumors.


## 4.1. Ideal Input

For three materials, elastic parameters are:λf=35, μf=12.5, and γf=0.4 for fatty tissue; λg=50, μg=25, and γg=0.25 for glandular tissue, and λd=80, μd=35, and γd=1.5 for ductal carcinoma in situ. Table 1 shows the initial estimate and reconstructed results for nonlinear elastography. The results in the first part are based on the ideal input. It is demonstrated that the reconstructed results are very close to the real values. The maximum error is 1.916% (λ for tumor) since the effect of the tumor on surface force measurement is the smallest. A similar situation occurs in linear elastography reconstruction.Table 1
Initial estimate and reconstructed results for nonlinear elastogrpahy. The results in the first part are based on ideal input (without noise). The ones in the second part are based on input with 5% noise and regularization is not used. The third part is based on input with 5% noise and regularization is applied to reduce the impact of noise.  (λ and μ are dimensionless, γ is in kPa.)
FattyGlandularTumorλfμfγfλgμgγgλdμdγdReal3512.50.450250.2580351.5Guess201012010120101Ideal inputRecon34.998812.50040.399950.051224.99980.249881.533134.94231.50035% Noise, without regularizationRecon22.27538.41471.503856.197021.54140.25090.000141.55622.00105% Noise, with exponential form regularizationRecon37.5629012.483180.459250.001425.00420.244480.000239.15391.5020


## 4.2. Adjoint Methods

Elastography includes forward and inverse problem. In forward problem, material parameters and loadings are given to calculate the deformation; while in inverse problem, external loadings and deformation are known to reconstruct material parameters. Most researchers established certain objective function and minimized it with a proposed iterative algorithm. The challenge is how to calculate the gradient of objective function efficiently and accurately. A straightforward calculation of gradients requires solving stiffness matrix in each iteration, which takes most of the time consumed in the finite element method.In this study an adjoint method is employed to analytically calculate the gradients. Oberai et al. [16] adopted an adjoint method and proposed a numerical scheme for reconstructing the nonuniform shear modulus field for incompressible isotropic materials using one component of displacement field. Liu et al. [10] applied this method for anisotropic materials. In this study, the adjoint method is applied for nonlinear elastography. The advantage of adjoint method is to solve two adjoint displacements during each iteration (w1 and w2 at (A4) in Appendix B), instead of the whole stiffness matrix, that increases the numerical efficiency significantly. In [18] Oberai et al. compared three different iterative methods: (1) a gradient-based method where the adjoint approach is used to calculate the gradient; (2) a gradient-based method where the straightforward approach is used to calculate the gradient, and (3) the Gauss-Newton method. The results show that “leading-order costs for the gradient-based method with the adjoint approach are smaller than the other two methods.” In fact, without the adjoint method, nonlinear elastography can only be discussed in concepts [23] or applied on simple objects using supercomputing power [24]. While in this study, material parameters in a real breast FEM phantom are reconstructed accurately by adopting this adjoint method. The results are encouraging for further clinical applications.


## 4.3. Multiple Sets of Measurements

Above results are based on four sets of force measurements on surface. For 2D isotropic elastography, Barbone and Bamber [25] has shown that one set of measurements for the displacements and forces, especially those taken only on the boundaries, may not provide sufficient information to the reconstruction of modulus distribution, because of the nonuniqueness nature of most inverse problems in elasticity. Barbone and Gokhale [26] further demonstrated the feasibility of using multiple displacement fields to reduce the likelihood of nonuniqueness for 2D isotropic elastography. Liu et al. [10] discussed the multiple sets of measurements in 3D anisotropic media. However research effort for 3D nonlinear elastography has not been found in literatures. In this study four independent tilted compression loadings are designed and surface force is measured to reconstruct internal material parameters. Due to different initial guesses, stable material parameter can be reconstructed, showing that the multiple sets of measurements are feasible to obtain the unique solution of inverse problem in nonlinear elastography.The key points for using multiple sets of measurements are to bring more deformation modes simultaneously into consideration. The loadings should be close to tumors in order to make tumors have larger deformation. In Figures5–7, four sets of loadings close to tumors are designed to obtain more deformation for reconstruction. This explanation serves as a guideline for design of loadings in clinical applications. Reduction of the number of required loadings will increase the clinical efficiency and benefit the patients. In other words, loadings with the richest stress-strain information are most desired. The design of feasible loadings is important for success in clinical application of nonlinear elastography.


## 4.4. Input with Noise

The above results are based on ideal input. However, noise cannot be avoidable in experiments. Its impact on reconstruction is therefore investigated with 5% of noise applied on surface measurements.  The results are shown in the second part of Table1. It is obvious that the reconstructed parameters are away from the real values. A possible reason is that, while the goal is to find real material parameters that minimize the objective function, a global minimum is not well defined with noise. In another words, real material parameters may not give a global minimum of the objective function with noise. Several similar parameters around real ones may give some local minimums. Our algorithm may reach one of them but fails to find real one because they all give local minimums. A typical way to solve this problem is to add a penalty term which provides additional constraint on the solution space. The penalty term may push the solution into the right area of the solution space and minimize the resulting objective.  A new objective function is therefore proposed as
(8)F=Φ+χΠ,
where Φ is the original objective function (3) which represents the difference between measured and calculated values. In addition, χ is the regularization factor and Π is the penalty term.Specific forms of penalty term have been designed for different problems [27]. In this study, an exponential form of penalty term is applied as
(9)Π=∑k=1K(1-exp(-(ξk-ak)2)),
where K=9 is the total number of material parameters, and ξk and ak are the reconstructed and true elastic parameters, respectively. If the true material parameters are unknown, ak can be estimated as close as possible. This is reasonable since several experiments have provided the nonlinear elastic parameter for breast tissues.The results are shown in the third part of Table1, demonstrating that the regularization method significantly improves the reconstruction data. Most of parameters are close enough to the true values. The largest error occurs for χ in fatty tissue with approximately 15%. By adding the penalty term, the reconstructed values are pulled from local minimum into the range close to true values. It is noted that the regularization factor χ varies in different scenarios. For example, χ can be smaller for higher confidence of experiment-data accuracy. On the contrary, if ak is closer to the true value, χ should be larger. For this study, the elastic parameters in fatty and glandular tissues are stable, comparing with the ones in tumor (Figure 4). Therefore, larger regularization factors are used for fatty and glandular tissues while smaller χ can be applied for tumors.


## 5. Conclusions

This paper presents a study on nonlinear elastography of biomedical imaging, in which a 3D model is developed for heterogeneous breast tissues extracting from real images including fatty tissue, glandular tissue, and tumors. Based on the large-deformation constitutive law, discretized nonlinear equations are solved for displacement, strain, and stress fields in breast tissues with given tumors under external compression at breast boundaries. A 3D inverse-problem algorithm is developed to reconstruct the material parameters for nonlinear elastic constitutive relation of breast phantoms with tumors. The adjoint gradient method is introduced to improve the numerical efficiency and enhance the stability of elastography reconstruction. Results demonstrate that this work opens the door toward nonlinear elastography, and provides guidelines for future developments and clinical applications in breast cancer study.



---

*Source: 406854-2009-07-09.xml*