Efficient Sensitivity Based Reconstruction Technique to Accomplish Breast Hyperelastic Elastography

Hyperelastic models have been acknowledged as constitutive equations which reliably model the nonlinear behaviors observed from soft tissues under various loading conditions. Among them, the Mooney-Rivlin, Yeoh, and polynomial models have been proved capable of accurately modeling responses of breast tissues to applied compressions. Hyperelastic elastography technique takes advantage of the disparities between hyperelastic parameters of varied tissues and the change in hyperelastic parameters in pathological processes. The precise reconstruction of hyperelastic parameters of a completely unknown pathology in the breast in a noninvasive and nondestructive way using the ultrasound elastography has been scrutinized in this paper. In the ultrasound elastography, tissue displacement field is extracted from radio frequency signals or images recorded using the ultrasound medical imaging system; hence the exact displacement field might not be obtained. Our results indicate that the parameters estimated by manipulating the iterative sensitivity-matrix based method converge to tissue's real hyperelastic parameters providing appropriate parameters are assigned to the hypothetical hyperelastic and regularization parameters. Iterative methods have therefore been proposed to compute proper hypothetical hyperelastic and regularization parameters. Accurate estimates of hyperelastic parameters of obscure breast pathology have been achieved even from imprecise measurements of displacements induced in the tissue by the ramp excitation.


Introduction
In the past two decades, substantial efforts have been made to bring the novel promising imaging approach "elastography" proposed by Ophir et al. [1] into clinical use. In contrast to the common medical imaging methods, elastography techniques could noninvasively provide qualitative or quantitative information on mechanical characteristics of biological tissues [2,3]. Distinct perceptible changes in the mechanical attributes of a biological soft tissue in varying pathologies yield the noninvasive detection and classification of its lesions using elastography techniques. Among them, ultrasound elastography approaches, which benefit from the superiority of conventional ultrasound (US) as being a safe, easy-to-use, inexpensive, nondestructive, noninvasive, widely available, and versatile medical imaging system, have received considerable attention.
In general, the practicable US elastography methods are classified on the basis of measured tissue-correlated physical quantities into [4]: (1) Strain Imaging: which is capable of qualitatively imaging Young's modulus, E, of the tissue by taking the response of tissue to a quasistatic load into consideration. (2) Shear Wave Imaging (SWI): which is capable of measuring shear wave speed, c , in the tissue (or its Young's modulus, E,) by taking the dynamic response of tissue to a mechanical vibration or acoustic radiation force into account.
Various praxes, particularly biomechanical characterization and imaging, have attributed the significant rise in the proposal of techniques for estimating linear and nonlinear elastic parameters of materials, i.e., for solving inverse 2 BioMed Research International problems in elasticity. Invaluable surveys of the relevant techniques have been provided by (i) Bonnet and Constantinescu [5], (ii) Hajhashemkhani et al. [6], (iii) Guchhait and Banerjee [7].
The aforementioned reviews have respectively discussed: (i) The formulations and solution methods, specifically the finite element and boundary element based methods, which are pertinent to general numerical methods for solving unspecialized elasticity problems on complex configurations, namely, least-squares functionals, adjoint solutions, and error in constitutive equation (ECE) based cost functions.
(ii) Momentous approaches for characterizing mechanical properties, i.e., hyperelastic parameters, of biological soft tissues and rubber-like materials with regard to the testing methods utilized, e.g., indentation tests, pipette aspiration experiments, propagation of elastic waves in the medium, inflation experiments, equibiaxial tensions, and in vivo experiments.
(iii) Prime elastic parameter estimation techniques involving linear elastic regime and geometrically or materially nonlinear forward model, for example, least-square based cost functionals, sensitivity based tactics, and integral approaches.
Besides the clinical applications, the introduction of numerous innovative strategies to estimate the nonlinear elastic parameters of materials or to reduce their computational costs, for instance, (i) virtual fields method (VFM) [8,9], (ii) subdomain inverse finite element technique [10], (iii) pointwise identification approach [11], (iv) gradient-based quasi-Newton minimization, adjoint, and continuation strategies [12], (v) the minimization-based reconstruction technique enhanced through material parameter grouping and user-supplied gradients of the objective function together with a nonlinear adjoint method [13], to name but a few, infers the importance of the subject of study. As regards our research focus, in brief, a least-square based cost functional was primarily applied by Iding et al. in 1974 [14] to estimate hyperelastic parameters of homogeneous materials. With respect to the geometrically nonlinear response of soft tissues at large deformations, i.e., the nonlinearity of the strain-displacement relation, an integral approach has been proposed by Skovoroda et al. [15] to reconstruct the elasticity distribution of soft tissue. Sensitivity based approaches have been applied formerly by Gendy and Saleeb [16] and recently by Hajhashemkhani and Hematiyan [6,17,18] to estimate hyperelastic parameters of rubber-like materials and soft tissues. The slope-variation method and Nelder-Mead algorithm have been applied by O'Hagan and Samani [19,20] and Naini et al. [21] to evaluate hyperelastic parameters of abnormal breast and deflated lung tissues.
With the aim of diagnosing breast cancers through reconstructing the spatial distributions of linear and nonlinear elastic parameters in patients with benign and malignant tumors [22,23], the inverse nonlinear elasticity problem has been altered into a minimization problem by Gokhale et al. [24]. The gradient-based quasi-Newton optimization method has been applied to minimize the cost function in consideration of the spatial distribution of material properties. The adjoint elasticity equations and continuation (in the material properties) scheme have been employed to calculate the gradient in reasonable time.
The application of US elastography technique to reconstruct nonlinear elastic constants of normal and abnormal breast tissues has been intended in this paper. The response of breast tissue, with and without the lesion, to a ramp stimulus (with low rate of increase in the applied load to ignore the inertia effect) has been simulated with the help of the finite element (FEM) software, Abaqus FEA. Two iterative methods founded on the stress-strain relation and sensitivity matrix have been, respectively, applied to estimate proper hypothetical and precise real hyperelastic parameters for the unknown tissues from limited displacement quantities in the tissues.
The approximate estimation of tissue displacement fields from the simulated US radio frequency (RF) signals, using the cross-correlation algorithm, impelled us to determine proper regularization parameter to converge to the real hyperelastic parameters of the tissue. An analogous iterative tactic has therefore been employed to compute proper regularization parameters. The decrease in the errors of elastic parameters estimated for the tumor, by comparing with the real elastic parameter estimated for the tumor, is the essence of the proposed iterative methods which lead to appropriate hypothetical and accurate real hyperelastic constants and suitable regularization parameter.

Hyperelastic Constitutive Models.
The results of multitude experimental scrutinizations of the behaviors of biological soft tissues, such as the breast and its lesions, have confirmed their nonlinear responses to applied stresses [25][26][27]. The researchers in the continuum mechanics work towards constructing mathematical models (mathematical equations) which could [28][29][30] (i) realistically represent the behavior of the understudy material; (ii) assess the material's response to the applied load; (iii) differentiate materials.
Hyperelastic constitutive models have made us competent at representing the nonlinear elastic responses of soft tissues to (large) strains. The constitutive theory of hyperelastic materials regards both the nonlinearity in the material behavior and considerable changes in the shape of material [31][32][33].
The strain energy density function, aka stored energy function, characterizes the energy absorbed by the homogeneous material in consequence of its deformation. Certain strain energy density functions are utilized to describe hyperelastic materials. As regards the deformation gradient tensor, F, (1) defines the strain energy function, , as a function of F, The invariants of deformation, F, aka strain invariants of deformation, make mapping the area and volume between the deformed and reference configurations possible. The first, second, and third invariants of deformation, 1 , 2 , and 3 , are computed by the use of (2) for unconstrained isotropic elastic materials, The left and right Cauchy-Green deformation tensors, B and C, are computed using The principal invariants of B and C are calculated as follows: As represented for B in (7), for incompressible materials, different sets of principal invariants of B and C are applied, Equation (8) specifies the Cauchy stress tensor of an unconstrained isotropic elastic material in terms of strain invariants, 1 , 2 , and 3 , The comprehensive compilations of the relations associated with the hyperelasticity theory have been provided by Bower [34], Felippa [35], and Holzapfel and Ogden [36]. On account of the unity of 3 (and therefore 0 =0), the Cauchy stress tensor relationship for incompressible materials simplifies to where refers to the arbitrary hydrostatic pressure, i.e., the Lagrange multiplier which compels the incompressibility constraint.
With the aim of accurately modeling the nonlinear elastic behaviors observed from soft tissues, a multitude of strain energy functions have been introduced in the literature. The functions have been defined in terms of (i) strain invariants, for instance, 1 , 2 , and 3 or 1 , 2 , and J, (ii) hyperelastic parameters, .
The functions extend from the well-known long-established Neo-Hookean and Mooney-Rivlin models (originated respectively by Treloar [37] in 1943, and Rivlin and Saunders [38] in 1951) to the recently inaugurated models, as the ones proposed by Nolan et al. [39] in 2014, Chen et al. [40] in 2015, and Carniel and Fancello [41] in 2017. On the basis of the outcomes of several investigations, which have been referred in Section 3.1, three strain energy functions, While a uniaxial stress, , is applied to the medium, the deformation gradient tensor, F, is computed using (13), where every (i=1,2,3) refers to the principal stretch parallel to one of the coordinate axes, ] .
The principal stretches determine the principal invariants, 1 , 2 , and 3 , as represented, If the uniaxial stress, , applied to the medium is considered in line with the first coordinate axis, (15) defines the uniaxial strain, , produced in the medium due to the applied stress, where represents the parallel stretch. Through assuming (1) the incompressibility of the medium (i.e., 3 =1) and (2) the equivalence of deformations in the two other coordinate axes, the Cauchy stress tensor equation for an incompressible medium (9) is simplified to [42,43] Regarding the above explanations, (17), (18), and (19) represent the Cauchy stress-stretch relations of the Mooney-Rivlin, Yeoh, and second-order polynomial models, respectively, In consideration of (15), the following stress-strain relationships have been, respectively, achieved for the Mooney-Rivlin, Yeoh, and polynomial models,

Estimation of Hyperelastic Parameters.
The differences between the hyperelastic constants of normal and abnormal breast tissues, as thoroughly discussed in [19,20,44] for numerous ex vivo breast tissue samples, make the detection and identification of breast tumors through their hyperelastic parameters feasible. With the purpose of precisely estimating the parameters of the selected hyperelastic models, namely, (i) the Mooney-Rivlin parameters, i.e., 10 and 01 of (10), (ii) the Yeoh parameters, i.e., 10 , 20 , and 30 of (11), (iii) the polynomial parameters, i.e., 10 , 20 , 11 , 02 , and 01 of (12), for normal and pathological breast tissues, two iterative algorithms have been appraised. The algorithms are founded on (i) the stress-strain relationship and sensitivity matrix, which have been formed on the basis of the relation of the selected hyperelastic model, as explained in Sections 2.1 and 2.2, (ii) the noninvasive measurement of displacement and strain fields in the understudy tissue from the recorded US RF signals and images.
Two types of software, MATLAB5 (The MathWorks, Inc., Natick, Massachusetts, USA) and a FEM software, for instance, Abaqus FEA (Dassault Systèmes Simulia Corp., Johnston, RI, USA), should be bilaterally connected to make the automatic iterations of the propounded algorithms possible.
At first, the recommended iterative algorithm which progressively provides approximate estimates of hyperelastic parameters of an undiagnosed tissue is described in this section with regard to (1) the explanations and equations mentioned in Section 2.1, particularly the ones associated with the BioMed Research International 5 relations between stress sets, strain sets, and the parameters of the Mooney-Rivlin, Yeoh, and polynomial hyperelastic models; (2) the availability of precise information on mechanical characteristics, i.e., linear or nonlinear elastic parameters, of normal mediums that surround the tissue.
To the best of our knowledge, at least, the elastic parameters of almost all healthy soft tissues, mainly measured through minute ex vivo experiments, have been reported in the literature. For several soft tissues, in particular the breast, as mentioned previously, we could benefit from the accessibility to precise information on their nonlinear mechanical constants.
The above-mentioned algorithm could be described in the following steps: (a) The unknown tissue, i.e., tumor, and its surrounding mediums are imaged before/while responding to a ramp excitation (with low rate of increase in the applied load to negate the inertia effect) for a period of time by the employment of a clinical US imaging system.
(b) The registered precompression US images, loading specifications, and boundary conditions are delicately regarded to accurately simulate the tumor and its adjacent mediums with the help of the FEM software, Abaqus FEA. Further explanation of the simulation strategies has been provided at the end of the section.
(c) In the simulated specimen, the values of 1 Pa and 0.5 are, respectively, assigned to the elastic modulus and Poisson's ratio of the tumor to simulate an equivalent elastic tumor.
(d) The displacement fields at some consecutive step times are extracted (i) for the tumor, from the recorded US RF signals or images, for instance, by the use of the crosscorrelation algorithm; i.e., the exact displacement fields in the tumor at some instants are computed; (ii) for the simulated elastic tumor, employing the FEM software, Abaqus FEA.
(e) The real elastic modulus of the tumor, real , is computed using ( [17,18,45]) where Y real and D, respectively, represent the axial displacement values of some points of the tumor and elastic tumor at the specified moments.
(f) The tumor strain field could be roughly approximated from the displacement measurements. The estimated elastic modulus for the tumor, real , is used to calculate a set of stress values, , from a set of strain values, (which is formed with regard to the strain field in the tissue), through the linear elasticity relation, (24), known as Hook's law, In view of the obtained results, a set of arbitrary strain values could also be considered, although the selection of strain set based on the available information leads to the significant decrease in the number of iterations.
(g) The parameters of the elected hyperelastic models, namely, the Mooney-Rivlin, Yeoh, and polynomial models, are computed using the formed stress and strain sets and the relation between stress, strain, and the parameters of a hyperelastic model, as represented in Section 2.1, i.e., (20), (21), and (22), for instance, with the help of regression algorithms of the MAT-LAB software.
(h) The elastic parameter for the simulated tumor, sim , is calculated after assigning the estimated hyperelastic parameters to the tumor and measuring the axial displacement values of the appointed points at the chosen step times, Y sim , by substituting Y sim for Y real in (23). The error of the computed elastic parameter, sim , while it is compared with real , i.e., the real elastic modulus calculated for the tumor, is used to evaluate the estimated hyperelastic parameters for the tumor. The assumption of the unavailability of initial information about the tumor and its mechanical characteristics has compelled us to consider the minimum of the errors of the estimates of tumor elastic modulus the criterion, as represented, where is the number of iterations of the algorithm.
(i) Through correctly modifying the strain set specified in step (f) and iterating the algorithm from this step, the estimated hyperelastic parameters could even converge to the real hyperelastic parameters of the tumor. A strain set that is slightly greater or less than the strain values, which are roughly computed from the tissue displacement measurements, would be the most appropriate choice; therefore, in this step, the strain quantities should be, respectively, decreased or increased.
The errors of the elastic parameters estimated for the tumor might not decrease below the defined tolerance value, specifically while the strain set is elected arbitrary. The iterative algorithm based on the sensitivity matrix, as defined by Hajhashemkhani and Hematiyan [17,18], would be the right choice to converge to accurate estimates of tumor hyperelastic parameters. The attained results indicate that the selection of suitable initial hyperelastic parameters, which could be obtained by the use of the proposed iterative method, is imperative to converge to precise estimates of tumor hyperelastic parameters through the sensitivity-matrix based algorithm; otherwise, the algorithm might approach the local minima.
The following steps specify the sensitivity-matrix based algorithm which has been applied to estimate the parameters of the Mooney-Rivlin, Yeoh, and polynomial hyperelastic models: (j) Similar to the previous algorithm, the US images and RF signals recorded from the mediums which have surrounded the obscure tissue are considered. The data are being registered for a time period before/while the mediums are responding to a ramp excitation (with low incremental rate to annul the inertia) by the use of the clinical US imaging system.
(k) By dint of the FEM software, the tumor and its adjacent mediums are meticulously simulated considering the saved precompression US images, loading specifications, and boundary conditions. The simulation tactics have further been explained at the end of the section.
(l) The axial displacement values of the specified points of the tumor at the determinate step times, Y real , are extracted from the recorded US RF signals or images, for instance, by employing the cross-correlation algorithm. (o) The sensitivity matrix is constructed in this step, as follows: (1) The set of estimated hyperelastic parameters, C, is assigned to the tumor. In the first iteration of the algorithm, the set of hyperelastic parameters computed for the tumor by the employment of the previously described method is taken into account. The displacement quantities of the selected points of tumor at every specified moment, D ,C (in total, D), are regarded. (2) Based on the hyperelastic model considered for the tumor, the response of the tumor, after slightly altering each of the hyperelastic parameters, , for instance, about 0.1% of its value, is simulated; similarly, the displacement quantities of the same points of tumor at every determinate moment, D , , are extracted. The difference between the calculated displacement vectors at the corresponding step times, D , is computed for all the specified moments.
(3) The sensitivity matrix, S, [17,18] is constructed as follows: where , , and , respectively, represent the number of hyperelastic parameters of the model, the amount of change in the hyperelastic parameters, and the number of consecutive step times when the responses of the tissue, i.e., the displacement quantities of some points of the tissue, have been computed.
(p) With regard to the Tikhonov regularization method, [17,18] is then used to compute the new hyperelastic parameters, C est , for the tumor, where the parameter represents the regularization parameter. On account of the achieved outcomes, there is no need to define the regularization parameter when the tissue displacement fields have been accurately measured from the recorded US RF signals and images.
(q) The estimated hyperelastic parameters, C est , are applied to calculate the tumor elastic parameter, sim , and its error by comparing it with real , as explained in step (h).
(r) As regards the estimated hyperelastic parameters, steps (o) to (q) are repeated until the error of the elastic parameter estimated for the tumor decreases below the defined tolerance value, as illustrated, where and , respectively, represent the tolerance value and the number of iterations of the algorithm.
By reason of recording low-quality US RF signals or images, applying an imprecise motion tracking method, assigning inappropriate values to parameters of the motion tracking algorithm, for instance, the maximum lag in the cross-correlation algorithm, or other attributes, the displacements of selected points of tumor might not be measured correctly. It is required to determine an appropriate regularization parameter to converge to the tumor hyperelastic parameters while the iterative sensitivity-matrix based method is applied to the imprecise displacement measurements.
The errors of the elastic parameters estimated for the tumor could also be considered to allocate proper value to  the regularization parameter. With reference to the attained outcomes, through assigning a very small value to the regularization parameter, the estimated hyperelastic parameters by manipulating the sensitivity-matrix based algorithm could converge to the real hyperelastic parameters of the tumor. Since it is assumed that no initial knowledge of the tumor and its mechanical features is available, similar to the two previously described algorithms, right value for the regularization parameter could also be determined regarding the errors of the elastic parameters estimated for the tumor. Although obviously the meticulous simulation of the understudy tissue and its neighbors on the basis of the recorded predeformation US images, loading specifications, and boundary conditions by the use of the FEM software, Abaqus FEA, could significantly improve the results, the attained outcomes connote the stability of the estimates against the errors induced by the displacement measurement or other attributes. Furthermore, the simulation could be enhanced by the use of innovative strategies proposed to (i) form three-dimensional (3D) volume data from two-dimensional (2D) US scans; hence Voxel-Based Methods (VBM), Pixel-Based Methods (PBM), and Function-Based Methods (FBM) could be applied [46,47]; (ii) compensate for the incomplete knowledge of the boundary conditions; notably, the problem of unknown conditions on part of the boundary has been solved by Hajhashemkhani et al. [6] using the Gauss-Newton method to minimize the cost function defined on the basis of the measured and calculated displacements.

Soft Tissue Simulation.
With the aim of demonstrating the efficacy of the proposed method in estimating hyperelastic parameters of an unknown tissue, particularly an unidentified tumor in the breast, we have utilized the FEM software, Abaqus FEA, to simulate 3D breast tissue geometry. Three tissues, the fat, fibroglandular, and an interior spherical tumor, constitute the simulated breast tissue, which has been depicted in Figure 1. Figure 1 also represents the load applied to the medium and the defined boundary condition, i.e., Encastre. As previously explained, the iterative algorithms have been employed to estimate the parameters of three hyperelastic models, namely, the Mooney-Rivlin, Yeoh, and second-order polynomial models. The hyperelastic parameters allocated to the normal fat and fibroglandular breast tissues and benign and malignant breast tumors, i.e., fibroadenoma, invasive lobular carcinoma (ILC), and invasive mucous carcinoma (IMC), have been reported in Tables 1-4. It is assumed that the mechanical parameters, i.e., the linear and nonlinear elastic parameters, are constant throughout each tissue type. Table 1 demonstrates a set of elastic parameters and Mooney-Rivlin hyperelastic parameters of breast tissues which has been widely applied to simulate the breast [48][49][50]. The Yeoh and polynomial hyperelastic parameters, presented in Tables 2-4, have been reported by Samani and Plewes [44] (for the normal tissues) and O'Hagan and Samani [19] (for the benign and malignant tumors). An iterative method has been proposed by O'Hagan and Samani to estimate hyperelastic parameters of ex vivo breast tissue specimens from their responses to the lowfrequency sinusoidal load. Their results confirm that the Yeoh and polynomial models are the hyperelastic models which conform more to the experimental data recorded from the breast tissues [19,20].
In pursuance of the simulation of US RF signals and images by the use of the Field II US Simulation Program,  Figure 1: The load, boundary condition, and mesh applied to the simulated breast tissue which is composed of the normal fat and fibroglandular tissues and a spherical tumor.
as explained in Section 3.3, the positions of scatterers after applying the load to the simulated medium should be calculated. It is feasible to precisely compute their positions through augmenting the number of nodes in the medium. As depicted in Figure 1, 36303 4-node hybrid tetrahedron (C3D4H) elements with 7603 nodes have constituted the mesh of the simulated phantom. The accuracy of the simulation outcomes has been verified by the convergence analyses.

Estimation of Hyperelastic Parameters from Precise and
Imprecise Displacement Measurements. The displacement distribution, specifically in the axial direction, in the in vivo tissue could be measured from the data recorded by the use of the standard medical imaging systems such as the US or magnetic resonance imaging (MRI) system in a noninvasive way; therefore, in view of estimating hyperelastic parameters of the tumor, we have employed the axial displacement values of some points of tissue while the tissue is responding to the ramp stimulation. With the purpose of ignoring the inertia effect, the rate of increase in the applied load is considered low.
With regard to the displacement values of several points of tumor at some step times, the estimates of parameters of the Mooney-Rivlin, Yeoh, and second-order polynomial models, achieved by the use of the iterative methods explicated in Section 2.2, have briefly been represented in Tables 5-10. The results have been reached by applying the displacements of maximum fifteen points of tumor in maximum twelve consecutive instants (with the difference of 0.25 s between the step times). Higher values have been considered for the second-order polynomial model with five indeterminate parameters, particularly while the displacement measurements were inexact.
The results associated with estimating the hypothetical and real elastic parameters and hyperelastic parameters of the Mooney-Rivlin model for the normal and abnormal breast tissues by the use of the proposed iterative methods have been represented in Tables 5 and 6. The outcomes of the estimation of hypothetical and real hyperelastic parameters of the Yeoh model for the benign and malignant breast tumors, i.e., the fibroadenoma and ILC, by applying the suggested algorithms have been summarized in Tables 7, 8, and 9. Table 10 demonstrates the capability of the recommended techniques in properly computing the hypothetical and real hyperelastic parameters of the second-order polynomial model based on the response of the malignant IMC tumor to the ramp excitation.
The calculation of displacement distributions in the stimulated breast tissue from the pre-and postdeformation RF signals simulated by the use of the Field II US Simulation Program [51] (A MATLAB5 toolbox for US field simulation) has implied the necessity to consider the errors of displacement measurements. Table 11 indicates the competence of two propounded methods, i.e., the iterative stress-strain relationbased and sensitivity-matrix based algorithms, in properly computing the hypothetical and real hyperelastic parameters from the inexact estimates of tumor displacement fields.      Table 7).

Simulation of US RF Signals and Images.
The Field II US Simulation Program [51] has been applied to simulate the pre-and postdeformation RF signals (on the basis of the response of the breast tissue to the applied load) to appraise the errors of displacement measurements. While the RF signals of breast tumors surrounded by the normal fat and fibroglandular tissues were being simulated, zero acoustic impedance was allocated to the abnormal breast tissues. In the Field II US Simulation Program, the properties selected to model the probe array and simulate the US RF signals have been considered as follows: In addition, the positions of scatterers in each deformation state have been calculated by linearly interpolating the displacements of the adjacent nodes computed through the FEM software, Abaqus FEA. To illustrate, two sets of B-mode images of the breast tissue, i.e., the pre-and postdeformation B-mode images, constructed from the RF signals associated with the times of 7.50 s and 9.00 s after starting to apply the ramp excitation, have been represented in Figure 2.
The gold standard of displacement estimation from the US RF signals, i.e., the cross-correlation algorithm, has been applied to evaluate the errors of displacement measurements. The displacement distributions in the tumor at the selected step times have been calculated from the pre-and postdeformation RF signals. In comparison with the displacement measurement from the US images, more precise estimate of tissue displacement field with higher spatial resolution is obtained by utilizing the RF signals.
In view of declining the computation time, which significantly raises by the use of the RF signals due to their higher sampling rate, two tactics have been employed: (1) Through defining maximum lag in regard to the estimated displacements of adjacent regions, i.e., preceding samples and lines, each search area (the number of lags considered in specifying the maximum of the cross-correlation function for two corresponding windows of the RF signals) was reduced considerably [52]. (2) The computational cost of the cross-correlation algorithm was appreciably declined by eliminating repeated calculations with the help of precomputed sum tables [53].

Conclusions
The correct identification of benign and malignant tumors in the breast through their nonlinear elastic parameters by the use of a noninvasive and nondestructive method has   been intended in this paper. With respect to the principles of elastography technique, two successive iterative algorithms, founded on (1) the relation between stress, strain, and the parameters of a hyperelastic model, (2) the sensitivity matrix, which correlates the changes in the displacement field in the tissue to the variations of hyperelastic parameters, (3) the estimation of displacement and strain fields in the tissue from the recorded US RF signals and images using the motion tracking techniques, (as explained in the Materials and Methods section of the paper), have been utilized to precisely estimate the parameters of three hyperelastic models, (i) the Mooney-Rivlin model, i.e., the parameters 10 and 01 of (10), (ii) the Yeoh model, i.e., the parameters 10 , 20 , and 30 of (11), (iii) the second-order polynomial model, i.e., the parameters 10 , 20 , 11 , 02 , and 01 of (12).
The exact/inexact displacement values of limited points (at restricted instants) of the tumor excited by the ramp stimulus have been applied to calculate the parameters. The dependency of the proposed method to the displacement and strain quantities warrants its competence as a noninvasive stratagem. The displacement and strain fields in the in vivo tissue could be computed from the recorded US RF signals or images by the employment of motion tracking approaches, which have been typically classified into three categories [53], The reliance of a method to precise values of deformation variables except or further than the displacement or strain, as the ones proposed by Omidi et al. [54], Roy and Desai [55], Liu et al. [56], Boonvisut and Ç avuşoglu [57], and Wang et al. [58], to mention a few, prevents the technique from being considered an elasticity imaging approach with the capability to noninvasively depict the nonlinear elastic features of in vivo tissues.
The inadequacy of the displacement-based techniques, for instance, two methods proposed by Mehrabian and Samani [59][60][61] and Hajhashemkhani and Hematiyan [17,18], necessitates the improvement of the tactics or the introduction of novel strategies for the hyperelastic elastography. The main defects of the aforementioned methods (as some of them have been assessed in [62]) are as follows: (1) the dependency of the defined coefficient matrix in the former method to the precise displacement measurements of a great number of adjacent points inside the medium and the reliance of the latter method to the displacement values of some boundary points of the tumor, which might not be accurately extracted from the recorded data, e.g., the US RF signals or images, (2) the necessity to have initial knowledge of the tumor in order to (a) consider proper initial guesses for the hyperelastic parameters to initiate the algorithms and (b) be assured of converging to the main hyperelastic parameters, specifically on account of the defined criteria to stop the algorithms, (3) the requisite to employ appropriate regularization methods and parameters, for example, as indicated by Mehrabian and Samani, the truncated singular value decomposition (SVD), Tikhonov regularization, and Wiener filtering techniques [59][60][61], which have been rectified in the proposed method. The higher error of displacement estimates in the boundary region has been discussed in [63][64][65]. It has been assumed that no primary knowledge of the tumor is accessible except the exact or even approximate measurements of displacement and strain fields in the tumor. With consideration of the estimated hyperelastic parameters, the minimum error of the elastic parameters calculated for the tumor, while they are compared with the main elastic parameter computed for the tumor from the estimated displacement fields in the tumor, is the principal criterion for the choice of the best estimates of hypothetical and real hyperelastic parameters.
Thanks to the defined criterion, the repeated manipulation of the stress-strain relation of the hyperelastic model results in (i) converging to the proper initial estimates of the parameters of the hyperelastic model; therefore, the possibility of inaccurately estimating the hyperelastic parameters of the understudy tissue declines to zero; (ii) significantly reducing the number of iterations of the iterative sensitivity-matrix based method, i.e., the computational cost of the real hyperelastic parameter estimation algorithm.
In pursuance of fulfilling the breast hyperelastic elastography, the propounded technique will be applied to estimate the nonlinear elastic constants of in vivo breast abnormalities from the recorded US RF signals and images in the future research.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Disclosure
The funding had no role in the study design, data collection, or analysis.