 Research
 Open Access
 Published:
Reconstruction of elasticity: a stochastic modelbased approach in ultrasound elastography
BioMedical Engineering OnLine volume 12, Article number: 79 (2013)
Abstract
Background
The convectional strainbased algorithm has been widely utilized in clinical practice. It can only provide the information of relative information of tissue stiffness. However, the exact information of tissue stiffness should be valuable for clinical diagnosis and treatment.
Methods
In this study we propose a reconstruction strategy to recover the mechanical properties of the tissue. After the discrepancies between the biomechanical model and data are modeled as the process noise, and the biomechanical model constraint is transformed into a state space representation the reconstruction of elasticity can be accomplished through one filtering identification process, which is to recursively estimate the material properties and kinematic functions from ultrasound data according to the minimum mean square error (MMSE) criteria. In the implementation of this modelbased algorithm, the linear isotropic elasticity is adopted as the biomechanical constraint. The estimation of kinematic functions (i.e., the full displacement and velocity field), and the distribution of Young’s modulus are computed simultaneously through an extended Kalman filter (EKF).
Results
In the following experiments the accuracy and robustness of this filtering framework is first evaluated on synthetic data in controlled conditions, and the performance of this framework is then evaluated in the real data collected from elastography phantom and patients using the ultrasound system. Quantitative analysis verifies that strain fields estimated by our filtering strategy are more closer to the ground truth. The distribution of Young’s modulus is also well estimated. Further, the effects of measurement noise and process noise have been investigated as well.
Conclusions
The advantage of this modelbased algorithm over the conventional strainbased algorithm is its potential of providing the distribution of elasticity under a proper biomechanical model constraint. We address the modeldata discrepancy and measurement noise by introducing process noise and measurement noise in our framework, and then the absolute values of Young’s modulus are estimated through the EFK in the MMSE sense. However, the initial conditions, and the mesh strategy will affect the performance, i.e., the convergence rate, and computational cost, etc.
Introduction
In the daily clinical exercises, palpation is frequently used as one part of physical examination to determine suspicious lesion’s size, shape, firmness, or location [1], but this process is largely decided by the doctor’s experiences. Therefore, ultrasound elastography [2], which aims to replace subjective palpation process, is developed to provide quantitative measures of material properties of tissue. As an evolving imaging modality, ultrasound elastography has already been broadly evaluated in the clinical practices and shown great potentials in providing valuable diagnostic information, such as improving the diagnostic accuracy of breast and prostate cancer [3–9], assessing plaque vulnerability [10–13] and guiding minimally invasive therapy [14–16]. Similar elastography techniques have been also developed using magnetic resonance imaging [17–19] and mammographic imaging [20]. Ultrasound elastography can be broadly grouped into quasistatic ultrasound elastography and dynamic ultrasound elastography according to the mechanical excitation applied on the biological tissue [21]. During the quasistatic ultrasound elastography, which is the focus of this paper, the mechanical excitation is usually applied through the ultrasound transducer and this external mechanical stimuli can be considered as a steadystate quasistatic loading [2]. After the deformation caused by the external loading is captured by the ultrasound machine, the elastographic image can be generated using different reconstruction techniques [22]. During the dynamic ultrasound elastography [3, 23], harmonic excitations on the order of 10–1000 Hz, or transient dynamic excitations are applied to the biological tissue in order to generate the shear wave. After the shear wave of the biological tissue is captured by the ultrasound machine, the following images of shear modulus can be generated using similar reconstruction techniques [24].
In spite of varieties of ultrasound elastography, we focus on the estimation of mechanical elasticity in quasistatic ultrasound elastography, which is one popular elastographic technique because of its efficiency and conveniency. The imaging process of quasistatic ultrasound elastography can be described by following four steps: a radiofrequency (RF) echo frame is first acquired from the tissue; second, a small deformation is induced into the tissue through the probe; third, a second RF frame is acquired from the deformed tissue; fourth, the displacement field, strain field and the distribution of tissue elasticity are reconstructed. Numerous techniques have been proposed to solve the reconstruction problem in step four of quasistatic ultrasound elastography in the last decades [25]. According to the way to obtain the tissue elasticity, the existing reconstruction algorithms for step four in quasistatic ultrasound elastography can be classified into two groups: the direct approach and iterative approach [26, 27]. However, the estimation of stable and meaningful mechanical elasticity in ultrasound elastography still remains as a challenging researching topic at present because of its illposed nature [28]. In the recent review [22], the modelbased approach is considered as one promising researching direction in solving the reconstruction problem in step four, and how to define the meaningful model constraint during the quasistatic ultrasound elastography is the key issue in the modelbased approach.
According to the constitute law of elasticity, accurate quantification of material properties, such as Young’s modulus, is feasible only when the distribution of all components of the internal stress is known. In quasistatic ultrasound elastography, however, the distribution of internal stress is hardly available. Hence, in the direct approach, the elastographic image is reconstructed always by converting the strain value as relative Young’s modulus using a very simple model constraint (e.g., Hook’s law) [2]. Mostly, the distribution of relative Young’s modulus (e.g., inverse strain values), is directly interpreted as the elastographic image based on the assumption of stress uniformity. But the assumption of uniformed internal stress distribution is not always valid in clinical cases. Furthermore, poor contrasttransfer efficiency will incur because of the assumption of uniformed stress distribution. In [2], an artifact, called “target hardening”, caused by this simplified assumption has been discussed. In the following work [29], an analytic solution of the elasticity equations derived for an infinite medium subjected to a uniaxial compression induced by a finite size compressor was used to verify the fact that a gradual decay of stress with increasing distance along the axis and off axis would cause “target hardening”. A different approach later tried to compress the target hardening by reconstructing the elastogram with measured strain and predicted stress field [30]. However, this approach still contained artifacts in an inhomogeneous medium: bidirectional shadows appear at the boundaries of the inclusion and the quality of elastographic image is even worse in complex situations. Several other groups including Doyley et al. [31], have sought for better quality of elastographic images by exploring the feasibility of solving the inverse problem based on a theoretical representation of forward elasticity model without requiring the internal stress distribution, but it turned out to be difficult in recovering the elastic modulus by inverting the forward elasticity model. In [32] and [33], a similar technique for quantifying the tissue elasticity was proposed to invert a finite difference representation of forward elasticity model with the pressure terms disappearing in the equations elimination, but artifacts were still observed in the reconstructed images. In another approach [34], the pressure term was included, but it was found that the method failed to employ in practice since hydrostatic pressure is nearly impossible to be measured from within tissues. In a word, the direct approach is trying to simplify the inverse problem (i.e., reconstruction of elastographic images), by inverting the forward model directly using strict controlled conditions or simplified assumptions, but these conditions or assumptions are usually invalid in clinical situation. For example, the assumption of internal stress uniformity within tissue is seldom valid in clinical cases. Overall, the direct approach’s capability is limited to provide the relative elastic modulus, but the ultimate aim of elastography to reconstruct the mechanical elasticity of tissues. What’s more, the differentiating process to compute the strain map from noisy displacement also negatively degrades the accuracy of direct approach because of its onetimetransformation process though the direct method still attracts great attentions from many groups at present because of its efficiency.
To overcome the drawbacks of direct methods, an alternative approach to replace the direct inversion technique is the iterative approach, which can be considered as a repeated direct inversion process, but with much more reasonable constraints, until an optimal estimation of the distribution of meaningful elastic modulus is obtained. Two key features of the iterative approach are: an iterative procedure and an external constraint. In [35], a Newton Raphson iterative scheme and a Tikhanov regularization method was applied to stabilize the reconstruction in the presence of noise. A similar algorithm was adopted in [36] and [26] for solving the inverse problem to recover the elastic modulus using GaussNewton algorithm, but this algorithm would cost heavy computation because of forming a Jacobian matrix. In order to reduce the computational cost from the iterative procedure, the gradient through adjoint elasticity equations was calculated and then the direction of gradient descent was used to find the optimal estimation [37]. In order to better recover smooth kinematic functions, one iterative twodimensional phasematching method was developed to reconstruct a twodimensional displacement vector field during acquisition of two successive RFecho data frames [38]. Another work employed a splitandmerge strategy, which used the strain images to provide a priori knowledge of relative stiffness distribution to reduce the computation time and avoid singularity in solving the inverse problem, to estimate the elastic modulus over the image sequence [39]. In the direct approach, simplified constitutive models can be directly inverted, but easily destroyed by the noises from measurement and complex tissue deformation. Therefore, the mechanical models with the ability of selfcorrecting to minimize the cost function that measures the goodness of fit between computed and measured data are always adopted in the iterative approach. For example, a finite element (FE) model with appropriate mathematical description of the object was adopted as the constraint to reconstruct the modulus of elasticity in [23]; in prostate elastography or intravascular elastography where the compression direction is radial, one FEbased approach was developed to study the reconstruction of tissue elasticity for ultrasound elastography in the polar coordinates [40]; since the material properties of lesion is important to the FE model, one iterative approach that used a FE model was developed to reconstruct the Young’s modulus from the measured forcedisplacement slope [41].
The iterative approach with a biomechanical model constraint was even named as a modelbased method in [28]. In [28], a comparative evaluation of the modelbased method and direct approach was conducted, and a conclusion that better contrast transfer efficiency could be achieved particularly for the highcontrast inclusions. However, the modelbased method is fraught with three major defects: modeldata discrepancy, measurement noise and no unique result. The model mismatch and the measurement noise will undoubtedly compromise the quality and accuracy of the ensuing elastograms. The reason for nonunique solution is that there is no guarantee of producing unique modulus elastograms as the illposed nature of inverse elasticity problems. Nevertheless, with the assistance of iterative procedure and meaningful biomechanical model, the iterative approach should be able to perform better than direct methods for the recovery of absolute values of elastic modulus if three major defects can be properly handled.
To overcome the defects of the modebased approach mentioned above, we have introduced a stochastic modelbased algorithm in our previous paper [42]. Different from earlier approaches, this algorithm treats the displacement and Young’s modulus as random variables that need to be estimated together from ultrasonic RF signals: an EKF runs iteratively until it converges to the steady state of deformed tissue, and the optimal estimation of the full displacement field and Young’s modulus are computed in a MMSE sense. Furthermore, the modeldata discrepancy and measurement noise are treated as white noises after the biomechanical model is transformed into the stochastic state space. However, in [42], we only present qualitative comparison in the preliminary results. In this paper, three groups of data, synthetic data, phantom data and patient data, are used to show the ability of our framework in recovering the absolute values of elastic modulus. The results and quantitative comparisons prove that this algorithm is a practicable way to suppress artifact noise. In the end, one fact should be pointed out that it is straightforward to apply our strategy to other forms of quasistatic elastography, but only two dimensional quasistatic ultrasound elastography are studied in this paper because of limited space.
The rest of this paper is organized as follows: section ‘Methodology’ describes the details of our elasticity reconstruction method: the linear elasticity and its finite element implementation are introduced, followed by the discussion about the integration of the stochastic state space strategy and the dynamic equation for the quasistatic ultrasound elastography. In section ‘Experiments and results’, the results obtained from the experiments using simulated data, real phantom data and invivo clinical data are presented and discussed respectively. In each set of data, our results are compared to the strainbased maps. As the ground truths are available in simulation and phantom studies, the accuracy and robustness of our framework are evaluated by comparing the results generated by our strategy to the ground truths. In the clinical experiment, the elastographic images are compared to CT images with lesions indicated by the doctor. Finally, the discussion of our framework, concluding remarks and future research endeavors are outlined in section ‘Discussion’.
Methodology
Linear elasticity
In order to construct a meaningful, yet computationally feasible computing scheme, material properties of the tissue should be properly modeled. In general, the biological tissue is a nonrigid object that deforms over time and has very complicated material properties in terms of the underlying constitutive laws [43]. However, in this paper, we only study the quasistatic ultrasound elastography, which is to utilize a steady state quasistatic excitation on the tissues using the ultrasound transducer. Since the tissue is compressed slightly during quasistatic elastography, constitute law of linear elasticity is still valid [2]. Hence, for computational simplicity, in our current two dimensional implementation, the linear isotropic continuum material is adopted to describe the mechanical behavior of the tissue during the quasistatic elastography, and the relationship of the stress σ and the strain ε obeys the Hooke’s law [44]:
where S is the strainstress matrix.
For the two dimensional case discussed in this paper, let u(x, y) and v(x, y) be the displacement along the x and yaxis of a point, the infinitesimal strain tensor ε of the point is:
In twodimensional situation, under the plane strain condition [44], strainstress matrix S can be derived as:
S is a materialdependent matrix, E stands for the Young’s modulus, and ν stands for the Poisson’s ratio. Here, the Young’s modulus E and the Poisson’s ratio ν are two materialspecific parameters. This fact that the internal stress caused by the deformation is a function of the displacement vector and the material parameters is quite clear from these equations. Ideally, the problem could be tackled in three dimension to avoid the throughplane motion effect, but in this paper we only study the twodimensional problem because the images of quasistatic ultrasound elastography are two dimensional data collected by the linear array of ultrasonic probe. Further, a linear material model is used to illustrate the basic ideas and rationales of a joint estimation strategy to recover the distribution of real elasticity for quasistatic ultrasound elastography, but more realistic models can be easily adopted into current framework in the future works.
Stochastic finite element method
In the past decades, the deterministic finite element method has been able to provide an effective and convenient platform for biomechanical studies in the past decades [45, 46]. However, it does not have the capability to process the uncertainties of material properties, and kinematic observations, i.e., measurements of the displacement. Especially the ultrasound imagingderived data are usually corrupted by noises of various nature. It is thus necessary to develop a strategy, the stochastic finite element method (SFEM), which has been used widely for structural dynamic analysis in probability analysis frameworks [47, 48], to process the data from the ultrasound elastography. The main difference between SFEM and finite element method (FEM) is that material properties are deterministic variables in FEM, but are described by random fields, possibly with known prior statistics in SFEM.
In order to build our implementation of SFEM for the quasistatic ultrasound elastography, a Delaunay triangulated finite element mesh is constructed at the first image frame before compression. An isoparametric formulation defined in a natural coordinate system is used, where, for trinodal linear element, the basis functions are linear functions of the nodal coordinates [44]. With mesh constructed, assuming that the material parameters E (Young’s modulus) and ν (Poisson’s ratio) are temporally constant but varying spatially, we can have the following dynamic governing equation [44]:
where M, C and K are the mass, damping and system stiffness matrices, R is the load vector, and U is the displacement vector. In this work, U consists of u(x, y) and v(x, y). The calculation of M, K and C matrices are the standard FEM procedure, which can be read in [44]. Because the tissue density can be generally considered to be uniform over the region of interest, M is a known function of the material density and is temporally and spatially constant. K is a function of the constitutive law of material, which is linear elasticity here, and is related to the materialspecific Young’s modulus and Poisson’s ratio, which are considered as constants temporarily in this study. In our framework, these two local material parameters are treated as random variables with known a priori statistics, and will be estimated along with the motion parameters. The damping matrix C is frequency dependent, and we assume small proportional Rayleigh damping with C = α M + β K in our implementation. In Rayleigh damping α is the mass proportional Rayleigh damping coefficient and β is the stiffness proportional Rayleigh damping coefficient [44]. In practice, it is difficult to determine the damping parameters because they are frequency dependent. We assume Raleigh damping for the very low damping biological tissue during quasistatic elastography, and fix the two weighting constants to be 1%.
State space representation
In order to employ our filtering strategy to integrate biomechanical model, the dynamics equation (Equation (4)) needs to be transformed into a statespace representation of the continuoustime linear stochastic system. First the kinematic vector x(k) and the material parameter vector θ are defined as:
Where the kinematic vector x(t) is consisted of displacement U(t) and velocity \stackrel{\u0307}{U}\left(t\right) information, and the material parameter vector θ is consisted of Young’s modulus E and Poisson’s ratio ν. In general, tumor inclusion or tissue blocked from its blood nutrients is stiffer than normal tissue, which mostly reflects in the variation of E. Since benign and cancerous tumors usually have distinguishing elastic properties (i.e., different Young’s modulus value), the value of Poisson’s ratio can be fixed in the following implementation. For example, as one of incompressible materials, the Poisson’s ratio of tissue can be set close to 0.5.
The state space representation of Equation (4) becomes
where n(t) is the process noise which is an additive, zeromean, white noise (E[n(t)] = 0; E[n(t)n(s)^{′}] = Q_{ n }(t)δ_{ ts }, where Q_{ n } is the process noise covariance). Process noise cannot be ignored if modeldata discrepancy exists. The system matrices A_{ c }, B_{ c } and input forces W(t) are:
An associated measurement equation, which describes the observations y(t) extracted from the ultrasonic RF data, can be expressed in the form:
where e(t) is the measurement noise which is additive, zero mean, and white (E[e(t)] = 0;E[e(t)e(s)^{′}] = R_{ e }(t)δ_{ ts }, where R_{ e } is the measurement noise covariance), independent of n(t). Although the noise in the displacement measurement extracted from ultrasonic RF data might not be white Gaussian noise, the property of the noise from ultrasonic RF data is rarely available due to the complicated collecting process of ultrasonic RF data. Hence our best guess is to assume it is white Gaussian noise, which is one common treatment in engineering field [49–51]. Furthermore, we found that the performance of this assumption is satisfied in the following experiments. H is the measurement matrix which should be specified by the relation between state vector x(t) and measurement vector y(t).
In order to run the dynamic equation (Equation (4)) in computer, it should be transformed into a discrete statespace equation, typically seen in control and estimation literature [50, 52]. We discretize Equation (7) and (8) over a constant time interval T. Since the interval T is always a known constant, we can replace kT with k in following equations
where
The associated discrete measurement equation is:
where y(k) is the measurement vector contained the displacement extracted from ultrasound RF data. In most quasistatic elastography, only the axial components of displacement vector are extracted. Therefore, the corresponding place of D will be set to 1 or 0 according to the available measurement data. n(k) and e(k) are process noise and the measurement noise respectively.
In order to perform the estimation of the material parameter distribution, the state vector x is augmented by the material parameter vector θ to form the new state vector z:
Accordingly, the new augmented state equation is:
with
Where v(k) is the process noise, introduced by the model mismatch. The new augmented measurement equation is derived:
with
Using previous assumptions of noises, the augmented noises are also Gaussian with distribution:
Extended kalman filter (EKF) for joint state estimation
The reconstruction of quasistatic elastographic image now can be defined as calculating the optimal estimates of the nonlinear statespace estimation problem represented by the Equation (14) and (17) above. This leads to solve this nonlinear estimation problem using the EKF framework, which is to linearize the augmented state (Equation (14)) at each time step and perform a recursive procedure with natural block structure to estimate the material parameter. EKF operates like Kalman filter, adopting the same form of feedback control in estimation: the filter makes an estimation of the state at one moment and then obtains feedback in form of measurements. The operations of EKF can be divided into two steps: time update equations and measurement update equations. Time update equations are responsible for projecting forward the current state and error covariance estimates to obtain the priori estimates for the next time step, while the measurement update equations are responsible for the feedback with the posterior data.
Initializing the filter with \hat{z}{}^{}\left(0\right)={\hat{z}}_{0} and P\left(0\right)=\sum _{0}, the EKF runs sequentially as follows:

(1)
compute the prior estimation of the state, from t−1 to t:
\begin{array}{l}\hat{z}{}^{}\left(t\right)=f\left(\hat{z}\right(t1),w(t\left)\right),\end{array}(22) 
(2)
project the error covariance from t−1 to t:
\begin{array}{l}P{}^{}\left(t\right)={F}_{t}P(t1){F}_{t}^{T},\end{array}(23) 
(3)
compute the Kalman gain at t:
\begin{array}{l}G\left(t\right)={P}^{}\left(t\right){H}^{T}{\left(H{P}^{}\right(t){H}^{T}+{R}_{0})}^{},\end{array}(24) 
(4)
compute posterior estimation of the state with the measurement at t:
\begin{array}{l}\hat{z}\left(t\right)=\hat{z}{}^{}\left(t\right)+G\left(z\right(t)h(\hat{z}{}^{}\left(t\right)\left)\right),\end{array}(25) 
(5)
update the error covariance at t:
\begin{array}{l}P\left(t\right)=(IG(t\left)\right)\left(z\right(t)h(\hat{z}{}^{}\left(t\right)\left)\right),\end{array}(26) 
(6)
repeat (1)(5) to convergence. The needed quantities for Equations (22)–(26) are defined by:
\begin{array}{l}{F}_{t}=\frac{\partial}{\partial}f\left(z\right(t),w(t\left)\right){}_{z=\hat{z}}=\left[\begin{array}{cc}A\left(\hat{\theta}\right)& {M}_{t}\\ 0& I\end{array}\right]\end{array}(27)\begin{array}{l}H\left(\phantom{\rule{0.3em}{0ex}}\hat{z}\right(t\left)\right)=\frac{\partial}{\partial}h\left(z\right(t\left)\right){}_{z=\hat{z}}=\left[\begin{array}{cc}D& 0\end{array}\right]\end{array}(28)\begin{array}{l}{M}_{t}=\frac{\partial}{\partial \theta}\left(A\right(\theta \left)\hat{x}\right(t)+B(\theta \left)w\right(t\left)\right){}_{\theta =\hat{\theta}}\end{array}(29)
In our implementation, the initial error covariance matrix P(0) consists of four parts:
P_{1}(0), P_{3}(0) are the submatrices inside the initial error covariance matrix related to the trustworthiness of the kinematic functions and material parameters respectively, and P_{2}(0) is the kinematicsmaterial correlation submatrix, which is a zero matrix in this paper.
Summary of the algorithmic flow
As described above, the system vector z, which consists of kinematic vector x and material vector θ, is updated through the EKF filtering framework, but in this implementation we only have one set of data from quasistatic elastography, i.e., one frame of displacement measurement acquired in deformed configuration. In Figure 1, the flow chart of our filtering framework is illustrated. If the difference between the estimate of current time step and the estimate of previous time step is small enough, the operation of EKF will be terminated and the final estimate will be considered as close enough to the optimal estimation. Therefore, through our filtering framework, a quantification of the internal material properties distribution jointly with kinematic state could be recovered from the ultrasonic RF data collected from quasistatic elastography.
Experiments and results
In order to demonstrate the performance of our filtering framework, we design three groups of experiments using synthetic data, phantom data and clinical data respectively. All the codes in the following experiments are implemented in Matlab R2013a (MathWorks Corporation, USA) and run in a desktop computer with a core 2 intel CPU and 12G memory.
Synthetic data
The tissue examined in the quasistatic ultrasound elastography is a three dimensional object, but the elastographic image collected by the ultrasound probe is two dimensional data. Hence, in the first experiment, the scanned area of elastography is treated as one rectangular area, which consists of two distinguished parts (inclusion and background) with different materials parameters (E_{ y }=15KPa and ν_{ y }=0.47 for the inclusion part, and E_{ b }=7KPa and ν_{ b }=0.47 for the background part). In order to generate authentic ground truth, the mechanical deformation of the cylinder is simulated in ABAQUS (DS Simulia Corporation, USA) as the ground truth: a constant pressure P=1KPa is then applied on the top of the cylinder, and then the deformation ratio is calculated using the maximized displacement of the top boundary divided by the height of cylinder, which is 3% in this simulation. Then deformation of the rectangular crosssection of cylinder is extracted out from the simulation: the rectangular crosssection of cylinder is first divided into a triangular mesh consisted of 231 nodes and 400 triangular elements; then, the displacements and velocities of 231 nodes are exported from ABAQUS as the ground truth. Because the deformation in this experiment is small, we set rightangled equilateral triangle elements in the ABAQUS to do the FEM analysis. During the experiment of our modelbased algorithm using the synthetic data, Young’s modulus is treated as a constant inside each triangular element, but the value of Young’s modulus can vary from element to element. Hence, the values of Young’s modulus of those 400 elements are estimated jointly with kinematic functions through our filtering framework.
To evaluate the performance of the framework under different noises, two sets of simulated displacement data with different signaltonoise ratio (SNR) ((1) 55.1d B; (2) 29.9d B) are generated from the ground truth in this way: 1) simulated displacement data are exported from ABAQUS; 2) two levels of white Gaussian noise are added into the simulated displacement data. It is convenient to generate simulated displacement with different SNR in Matlab platform, such as using wgn function or awgn function. The axial strain maps calculated from noisy displacement data are first inverted as the relative Young’s modulus as the output of a conventional strainbased algorithm with the assumption of stress uniformity [2] as the comparison to the results of the following filtering experiments ((j) and (k) of Figure 2). The estimated results of Young’s modulus over triangular mesh using 55.1d B and 29.9d B data are shown in (b) and (c) of Figure 2 without smoothness. Interpolated results of Young’s modulus for 55.1d B and 29.9d B data are illustrated in (f) and (g) of Figure 2 respectively over the triangular mesh. In order to evaluate the effect of different noises on our framework, the estimated strain maps of our filtering framework are illustrated in (h) and (i) of Figure 2, in comparison to the noisy strain maps generated from the ground truth in (d) and (e) of Figure 2. In (j) of Figure 2, the middle line of (f) and (g) of Figure 2 and the deviation of neighboring zone (i.e., the distribution of modulus along the center zone), are illustrated. In (k) of Figure 2, the middle line and the deviation of neighboring zones of our filtering framework and the relative Young’s modulus from the conventional strainbased method are compared together. For two groups of simulated data using in our filtering method, the mean values of P 1 and P 3, parts of error covariance matrix P, are calculated in every iteration and displayed in Figure 3.
Phantom data
After testing our filtering framework in the simulated data, two sets of real ultrasonic RF data collected from different phantoms are used to evaluate the performance of our filtering framework in this section. The first set of data is collected in one Elasticity QA Phantom (Model 049, CIRS Inc. USA) using a PCbased Ultrasonix RP ultrasound machine (Ultrasonix Medical Corporation, Burnaby, BC, Canada). The true distribution of Young’s modulus of this set of data is: 80kPa for the inclusion and 25kPa for the background part. The axial components of displacement is extracted out from this set of data using a phaseshift method [53]. The second set of data is provided online in public, which is also collected from another Elasticity QA Phantom [54]. The true distribution of Young’s modulus of this set of data is: 55kPa for the inclusion part and 33kPa for the background part. The twodimensional displacement filed of the second set of data is also provided online in public by [54]. In the following experiment, the relative Young’s modulus is still calculated by inverting axial strain maps based on the assumption of stress uniformity [2]. In Figure 4, the estimated results of our modelbased algorithm, the strain image, the image of strainbased method, and the strain profiles along the center line of strainbased method and our method are displayed together for comparison A quantitative comparison of image contrast is also shown in Table 1. The numbers inside the brackets of Table 1 are mean values of Young’s modulus of the inclusion part and background part. In this table, the ideal contrast is compared to the actual contrast in different deformation ratios for two sets of phantom data.
Clinical data
Three sets of patients’ ultrasonic RF data are provided online in public [54]. These data were collected from patients undergoing open surgical RF thermal ablation for liver cancer who were enrolled between February 6, 2008 and July 28, 2009. The clinical statuses of these patients are described in detail [54]. The ultrasound experiments were performed with the approval of the Health Science Research Ethics Committee of Johns Hopkins University (JHU), and the participants provided written informed consent before beginning the experiment [54]. According to the previous reports on the mechanical property of liver tissue, the tissue with Young’s modulus above 6KPa can be considered as lesion inside liver [55]. All the patients’ data are processed in the same way as previous experiments in phantom data. Since the axial and lateral displacement measurements have been provided in [54], the extraction of displacement is not necessary in this experiment. The initial conditions for all three sets of data are not exactly the same throughout the experiment though the data acquisition is performed on the same kind of organ tissue using the same equipment. For example, noise matrices Q_{ v } and R_{ e } need to be set empirically. According to our experience in the phantom data, the convergence of our filtering framework is not sensitive to the initial values of elements of system error covariance matrix P because P can be automatically tuned during the EKF procedure as we have shown in Figure 3 in previous experiments.
In Figure 5, patients’ CT images scanned after RF ablation, strain images and our estimated results are illustrated respectively. In Figure 5, patients’ CT images scanned after RF ablation and strain images are provided online in public by [54]. In CT images, the details of the organ structure is visible, but it is difficult to tell the abnormal tissue from the normal tissue without clinical experiences. The position of tumors in CT images are marked by clinical experts after RF ablation. However, in the strain images provided by [54], it seems that the zones of abnormal tissue are much larger than the markers in CT images. However, the zones of abnormal tissue in our results are clearly smaller than those in the images generated by the strainbased method.
Discussion
In this paper, a filtering framework is proposed to recover the distribution of Young’s modulus (i.e., a meaningful elastographic image), from sparse and noisy ultrasonic RF data using a biomechanical constraint. The experiments reported in this paper reveal several interesting aspects of our filtering framework, and potentially useful characteristics to solve the recovery problem of elasticity distribution in quasistatic ultrasound elastography. In general, the quality of estimated results is mainly influenced by following factors: noise, deformation ratio, boundary conditions, input information and computational cost.
Two kinds of noises have been studied above: measurement noise and process noise. In most of previous iterative approaches, the measurement noise is considered, but the process noise is always ignored. The main contribution of this paper is to treat the measurement noise and process noise as Gaussian noises. The process noise is caused by the mismatch between the data and the system model, i.e., modeldata discrepancy. Since there are not much previous works to study the properties of process noise, the process noise is modeled as white noise distribution, which is the contribution of this paper for iterative approaches to quasistatic elastography to our knowledge. The measurement noise is also modeled as another independent white noise. Therefore, the tolerance of noises of our modelbased algorithm is one important factor that should be carefully examined. In the experiment of simulated data, as you can see in (j) of Figure 2, the tolerance of noise of our filtering framework gets worse when the noise becomes larger. Although our filtering method has incorporated the biomechanical constraints, the quality of estimation is still compromised by the noise as shown in (j) of Figure 2, which can be explained by the shortcoming of EKF: when the system is nonlinear, EKF might only be able to give suboptimal estimates, as it relies on linearization of the nonlinear system equation to propagate the mean and covariance of the state. But when the results of our modelbased algorithm are compared to the results of conventional strainbased method, as you can see in (k) of Figure 2, the relative elastogram (i.e., relative Young’s modulus calculated from the conventional strainbased method), is far away from the ground truth and easily corrupted by the noise, while our modelbased method still can recover a meaningful elastogram by introducing a biomechanical model into EKF framework. In Figure 4, the results of strainbased method and our filtering framework for the phantom data are displayed and compared. As shown in (d) and (h) of Figure 4, our filtering framework still can generate meaningful distribution of Young’s modulus, but the conventional strainbased method failed to recover the true value of Young’s modulus from ultrasonic data. A quantitative comparison is also shown in Table 1: the values calculated by the ground truth are close to the values calculated by our estimated results. Hence, we can conclude that our filtering method can generate more meaningful results under the assumption of Gaussian noise, which is to say that our filtering method is more robust than the conventional strainbased method in handling the noise from ultrasound elastography.
One factor to compromise the quality of the results of our modelbased algorithm is the measurement noise. Measurement noise is produced during the process of acquiring RF data from organ tissue. It is therefore necessary to take the measurement noise effect into account during the recovery procedure. In the past, some of Direct approaches ignore the measurement noise, and their results are drastically affected. The measurement noise introduced during the acquiring process of ultrasonic RF data mainly depends on ultrasound machine and the skill level of operators. One way to increase the SNR is to employ a large deformation ratio during the elastography. That is revealed in phantom study: with larger deformation ratio, the result of the second set of data is better. As shown in the Table 1, the quality of estimation is comprised by the noise largely under a small deformation ratio, therefore a large deformation ratio can be helpful for a quick and quantitive identification of abnormal tissue as shown in Table 1. However, it is not always practical to increase the deformation ratio arbitrarily because of decorrelation of RF signal, so the reduction of the noise effect is still an important issue that still needs to be addressed. Moreover, as shown in (d) and (h) of Figure 4, and Table 1, though our filtering framework is still able to recover a reasonable distribution of Young’s modulus from real ultrasonic data, the quality of estimated results is not so good as the results in the simulated data because the statistical property of system and measurement noises are available in controlled simulation and not easily available during the actual process of elastography. Therefore, the statistical property of process noise (i.e., modeldata discrepancy), and measurement noise are worthy to be carefully examined in the future work.
The way to enforce the boundary conditions is also important in our filtering framework. In this work, we have employed three different strategies to enforce the boundary conditions: known surface force distribution (simulated data), penalty method [44] and modifying stiffness matrix K with boundary nodes’ displacement information [44]. With known surface force distribution, it is convenient to enforce this condition through external loading R. However, surface force distribution introduced by the probe is difficultly acquired in clinical cases. In the penalty method, a large constant, the penalty factor, is added to the diagonal elements of stiffness matrix K, but it will change the structure of the statespace equation and easily diverge our modelbased algorithm. One possible way is to use the movement of boundaries of the interested area as the boundary conditions, so the external force conditions can be implicitly enforced by enforcing the movement of boundaries of the interested area. In the phantom and clinical data study, the stiffness matrix K is modified by the boundary nodes’ displacement information, which is an alternative way to enforce the boundary conditions without changing the structure of the statespace equation. In our experiments of phantom and clinical data, promising results have proved the efficiency and accuracy of modifying stiffness matrix K with boundary nodes’ displacement information. However, the displacement of moving boundary is extracted out from noisy ultrasonic data. So the quality of estimated results of our filtering framework can be affected by the noise of moving boundary from ultrasonic data.
Our filtering method is to generate the optimal estimates using the biomechanical model constraint and the available measurements. In the experiment of phantom data, as explained above, a full displacement field is provided in the second set of data, but only the axial displacement is given in the first set of data. Since the measurement from the second set of data can provide more reliable observations, our filtering framework can generate a better result in the second set of data than that in the first set of data as shown in Figure 4. Hence, in order to maximize the performance of our filtering framework, more reliable observations (i.e., input information to our filtering framework), should be secured. We should point out that the strain and elasticity images without anatomical information will prevent people from upstanding its benefits in clinical practice. Hence, we believe that a further study of the registration of anatomical images (such as CT) and functional images (strain and elasticity images) should be carried on. We admit that it is so difficult to obtain the ground truth in the clinical experiment, and thus it will be always a challenging and interesting topic to carry further study of the performance of our framework using different clinical data.
The initial conditions for all three sets of data are not exactly the same throughout the experiment though the data acquisition is performed on the same kind of organ tissue using the same equipment. For example, noise matrices Q_{ v } and R_{ e } need to be set empirically. According to our experience in the phantom data, the convergence of our filtering framework is not sensitive to the initial values of elements of system error covariance matrix P because P can be automatically tuned during the EKF procedure as we have shown in Figure 3 in previous experiments. For system error covariance matrix P, two important parts, P 1 and P 3, should decrease when our filtering method converges. As shown in Figure 3, the mean value of P 1 and P 3 decreases through iterations as expected, which indicates that the stochastic dynamic vector and material parametric vector becomes close to the optimal estimates. However, the convergence rate of mean value of P 1 and P 3 in small noise condition is higher than under higher noise condition, which is to say that the final converged state is still affected by the noise.
The computational cost is always a problem to the iterative approaches. In our filtering framework, the multiplication of large matrices is a timeconsuming task during the iterative process of EKF. In our study, the mesh size is 21*11. The size of state vector is 1324*1, and the size of error covariance matrix is 1324*1324. One loop (Eqs. (22)–(26)) costs about 58 seconds, and 200 or more loops are required to reach the convergence of EKF as show in previous Figure 3. It totally takes 1525 minutes to converge, which is not suitable for realtime data acquisition. Therefore, it is not feasible to use fine mesh in our framework in current hardware configuration. However, after incorporating the biomechanical model as the constraint into our filtering framework, it is clear in Figure 5 that the elastogram map with low sampling rate is much more readily identifiable than the strain map with high sampling rate: the size of lesions are much smaller than that indicated from strainbased images. So considering the quality of our filtering framework for the quasistatic ultrasound elastography, the reduction of computation, such as GPU acceleration of matrix operation, will be a promising researching topic in the future.
Conclusions
In our filtering framework, the absolute values of Young’s modulus are estimated through the EFK in a MMSE sense. We address the modeldata discrepancy and measurement noise by introducing process noise and measurement noise in our framework. Three groups of data: synthetic data, phantom data and patient data, are used to validate the performance of our framework. It is straightforward to apply this framework at other occasions of transformation from noisy dynamic information to material parameter distribution, not limited to mere quasistatic ultrasound elastography.
Authors’ information
Minhua Lu and Heye Zhang contributed equally to this work and share first authorship.
Abbreviations
 MMSE:

Minimum mean square error
 EKF:

Extended Kalman filter
 RF:

Radio frequency
 SFEM:

Stochastic finite element method
 FEM:

Finite element method
 FE:

Finite element
 SNR:

Signaltonoise ratio
 JHU:

Johns Hopkins University.
References
Samani A, Zubovits J, Plewes D: Elastic moduli of normal and pathological human breast tissues. Phys Med Biol 2007, 52: 1565–1576. 10.1088/00319155/52/6/002
Ophir J, Cespedes I, Ponnekanti H, Yazdi Y, Li X: Elastography: A quantitative method for imaging the elasticity of biological tissues. Ultrason Imaging 1991, 13(2):111–134.
Barr R, Stamatia D, Logan BL, William ES, Corinne B, Carmel S: Evaluation of breast lesions using sonographic elasticity imaging: a multi center trial. J Ultrasound Med 2012, 31(2):281–287.
Thittai A, Yamal J, Mobbs L, KraemerChant C, Chekuri S, Garra B, Ophir J: Axialshear strain elastography for breast lesion classification: further results from in vivo data. Ultrasound Med Biol 2011, 37(2):189–197. 10.1016/j.ultrasmedbio.2010.11.001
Souchon R, Rouviere O, Gelet A: Visualisation of HIFU lesions using elastography of the human prostate in vivo: Preliminary results. Ultrasound Med Biol 2003, 29(7):1007–1015. 10.1016/S03015629(03)000656
Hiltawsky K: Kruger C M adn Starke: Freehand ultrasound elastography of breast lesions: clinical results. Ultrasound Med Biol 2001, 27(11):1461–1469. 10.1016/S03015629(01)004343
Cespedes E: Elastography: imaging of biological tissue elasticity. PhD thesis, Houston. 1993.
Garra B, Cespedes E, Ophir J: Elastography of breast lesions: Initial clinical results. Radiology 1997, 202: 79–86.
Samani A, Bishop J, Plewes D: A constrained modulus reconstruction technique for breast cancer assessment. IEEE Trans Med Imaging 2001, 20: 877–885. 10.1109/42.952726
Brusseau E, Fromageau J, Finet G, Delachartre P, Vray D: Axial strain imaging of intravascular data: Results on polyvinyl alcohol cryogel phantoms and carotid artery. Ultrasound Med Biol 2001, 27(12):1631–1642. 10.1016/S03015629(01)004513
Korte C, Pasterkamp G, van der SteenA, Woutman H, Bom N: Characterization of plaque components with intravascular ultrasound elastography in human femoral and coronary arteries in vitro. Circulation 2000, 102(6):617–623. 10.1161/01.CIR.102.6.617
de Korte C: Morphological and mechanical information of coronary arteries obtained with intravascular elastography. Feasibility study in vivo. Eur Heart J 2002, 23(5):405–413. 10.1053/euhj.2001.2806
Doyley M, Mastik F, de Korte C: Advancing intravascular ultrasonic palpation toward clinical applications. Ultrasound Med Biol 2001, 27(11):1471–1480. 10.1016/S03015629(01)004574
Kallel F, Stafford R, Pric R: The feasibility of elastographic visualization of HIFUinduced thermal lesions in soft tissues. Imageguided highintensity focused ultrasound. Ultrasound Med Biol 1999, 25(4):641–647. 10.1016/S03015629(98)001847
Righetti R, Kallel F, Stafford R: Elastographic characterization of HIFUinduced lesions in canine livers. Ultrasound Med Biol 1999, 25(7):1099–1113. 10.1016/S03015629(99)000447
Varghese T, Techavipoo U, Liu W: Elastographic measurement of the area and volume of thermal lesions resulting from radiofrequency ablation: Pathologic correlation. AJR Am Roentgen Ray Soc 2003, 181(3):701–707. 10.2214/ajr.181.3.1810701
Sumi C, Nakayama K: A robust numerical solution to reconstruct a globally relative shear modulus distribution from strain measurements. IEEE Trans Med Imaging 1998, 17: 419–428. 10.1109/42.712131
Plewes D, Bishop J, Samani A, Sciarretta J: Visualization and quantification of breast cancer biomechanical properties with magnetic resonance elastography. Phys Med Bioly 2000, 45: 1591–1610. 10.1088/00319155/45/6/314
McGee K, Lake D, Mariappan Y, Hubmayr ARDandManduca, Ansell K, Ehman R: Calculation of shear stiffness in noise dominated magnetic resonance elastography data based on principal frequency estimation. Phys Med Biol 2011, 56(14):4291–4309. 10.1088/00319155/56/14/006
Miga M: A new approach to elastography using mutual information and finite elements. Phys Med Biol 2003, 48: 467–480. 10.1088/00319155/48/4/304
Varghese T: Quasistatic ultrasound elastography. Ultrasound Clin 2009, 4(3):323–338. 10.1016/j.cult.2009.10.009
Doyley M: Modelbased elastography: a survey of approaches to the inverse elasticity problem. Phys Med Biol 2012, 57(3):R35—R73.
Parker K, Huang S, Musulin R, Lerner R: Tissue response to mechanical vibrations for “sonoelasticity imaging”. Ultrasound Med Biol 1990, 16: 241–246. 10.1016/03015629(90)90003U
Barr R: Sonographic breast elastography: a primer. J Ultrasound Med 2012, 31(5):773–783.
Ophir J, Alam S, Garra B, Kallel F, Konofagou E, Krouskop T, Merritt C, Righetti R, Souchon R, Srinivasan S, Varghese T: Elastography: imaging the elastic properties of soft tissues with ultrasound. J Med Ultrasound 2002, 29: 155–171. 10.1007/BF02480847
Doyley M, Meaney P, Bamber J: Evalation of an iterative reconstruction method for quantitative elastography. Ultrasound Med Biol 2000, 45: 1521–1540.
Oberai AA, Gokhale NH, Doyley MM, Bamber JC: Evalation of the adjoint equation based algorithm for elasticity imaging. Phys Med Biol 2004, 49: 2955–2974. 10.1088/00319155/49/13/013
Doyley M, Srinivasan S, Pendergrass S, Wu Z: Comparative elvaluation of strainbased and modelbased modulus elastography. Ultrasound Med Biol 2005, 31(6):787–802. 10.1016/j.ultrasmedbio.2005.02.005
Ponnekanti H, Ophir J, Cespendes I: Axial stress distribution between coaxial compressors in elastography: an analytical model. Ultrasound Med Biol 1992, 18(8):667–673. 10.1016/03015629(92)90117S
Ponnekanti H, Ophir J, Huang Y: Fundamental mechanical limitations on the visualization of elasticity contrast in elastography. Ultrasound Med Biol 1995, 21(4):533–543. 10.1016/03015629(94)001362
Doyley M, Bamber J, Shiina T, Leach M: Reconstruction of elastic modulus distribution from envelope detected Bmode data. In. IEEE Ultrasonics Symposium 1996, 1611–1614.
Skovoroda A, O’Donnell M: Tissue elasticity reconstruction based on ultrasonic displacement and strain images. IEEE Trans Ultrason Ferroelectr Freq Control 1995, 42(4):747–765.
Li J, Cui Y, Kadour M, Nobel J: Estimation of shear modulus distribution in soft tissue from strain distribution. IEEE Trans Biomed Eng 1995, 44(2):193–202.
Raghavan K, Yagle A: Forward and inverse problems in elasticity imaging of soft tissues. IEEE Trans Nuclear Sci 1994, 41(4):1639–1648. 10.1109/23.322961
Kallel F, Bertrand M: Tissue elasticity reconstruction using linear perturbation method. IEEE Trans Med Imaging 1996, 15: 299–313. 10.1109/42.500139
Van Houten EE, Paulsen K, Miga M, Kennedy F, Weaver J: An overlanpping subzone technoque for MRbased elastic property reconstruction. Magnetic Reson Med 1999, 42: 779–786. 10.1002/(SICI)15222594(199910)42:4<779::AIDMRM21>3.0.CO;2Z
Oberai A, Gokhale N, Feijoo G: Solution of inverse elasticity problems in elasticity imaging using the adjoint method. Inverse Problem 2003, 19: 293–313.
Sumi C: Fine elasticity imaging utilizing the iterative RFecho phase matching method. IEEE Trans Ultrason Ferroelectr Freq Control 1999, 46: 158–166.
Li J, Cui Y, Kadour M, Nobel J: Elasticity reconstruction from displacement and confidence measures of a multicompressed ultrasound RF sequence. IEEE Trans Ultrason Ferroelectr Freq Control 2008, 55(2):319–326.
Luo J, Ying K, Bai J: Elasticity reconstruction for ultrasound elastography using a radial compression: an inverse approach. Ultrasonics 2006, 44: E195—E198.
Samani A, Plewes D: An inverse problem solution for measuring the elastic modulus of intact ex vivo breast tissue tumours. Phys Med Biol 2007, 52: 1247–1260. 10.1088/00319155/52/5/003
Wang J, Zhang H, Lu M, Liu H, Hu Z: Elastographic image reconstruction: A stochastic state space approach. In. The International Conference on Image Processing 2011, 2017–2020.
Humphrey J: Review Paper: Continuum biomechanics of soft biological tissues. Proc R Soc A 2003, 459: 3–46. 10.1098/rspa.2002.1060
Bathe K: Finite Element Procedures in Engineering Analysis.. Prentice Hall; 1982.
Benayoun S, Ayache N: Dense nonrigid motion estimation in sequences of medical images using differential constraints. Int J Comput Vis 1998, 26: 25–40. 10.1023/A:1007932523504
Creswell L, Moulton M, Wyers S, Pirolo J, Fishman D, Perman W, Myers K, Actis M, Szabo B, Pasque M, RLand Vannier: An experimental method for evaluating constitutive models of myocardium in in vivo hearts. Am J Physiol 1994, 267(2):H853—H863.
Contreras H: The stochastic finite element method. Comput Struct 1980, 12: 314–348.
Liu W, Belytschko T, Man A: Random field finite element. Int J Numerical Methods Eng 1980, 23: 1831–1845.
Simon D: Optimal State Estimation: Kalman, Hinfinity, and Nonlinear Approaches. John Wiley and Sons; 2006.
Glad T, Ljung L: Control Theory. London: Taylor and Francis; 2000.
BarShalom Y, Li X, Kirubarajan T: Estimation with Applications to Tracking and Navigation. New York: Wiley; 2001.
Shi P, Liu H: Stochastic finite element framework for simultaneous estimation of cardiac kinematic functions and material parameters. Med Image Anal 2003, 7: 445–464. 10.1016/S13618415(03)000665
Yuan J, Zhang H, Lu M, Chen S, Liu H: A strainbased ultrasound elastography using phase shift with prior estimates and meshfree shape function. In. IEEE International Symposium on Biomedical Imaging 2011, 532–535.
Rivaz H, Boctor E, Choti M, Hager G: RealTime Regularized Ultrasound Elastography. IEEE Trans Med Imaging 2011, 30(4):928–945.
Mueller S, Sandrin L: Liver stiffness: a novel parameter for the diagnosis of liver disease. Hepatic Med: Evid Res 2010, 2: 49–67.
Acknowledgements
This work was supported in part by the National Basic Research Program 973 (2010CB732504), the Guangdong Innovation Research Team Fund for LowCost Healthcare Technologies in China, the External Cooperation Program of the Chinese Academy of Sciences (GJHZ1212), the Key Lab for Health Informatics of Chinese Academy of Sciences, the Distinguished Youth Talents Fund of Shenzhen (JC201006020025A), the Natural Science Foundation of Zhejiang(R12F030004), and the National Natural Science Foundation of China (81101120, 61271083, 61031003, 81271645).
Author information
Authors and Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
ML and HZ contributed equally to this work and the manuscript writing as they both conceived of the study, designed the whole work, and drafted the manuscript. JW carried out the experiments on phantom and synthetic data. JY helped to extract the displacement information from clinical data. ZH participated in the design of this study and coordination and helped to draft the manuscript. HL helped to draft the manuscript. All authors read and approved the final manuscript.
Minhua Lu, Heye Zhang contributed equally to this work.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Lu, M., Zhang, H., Wang, J. et al. Reconstruction of elasticity: a stochastic modelbased approach in ultrasound elastography. BioMed Eng OnLine 12, 79 (2013). https://doi.org/10.1186/1475925X1279
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/1475925X1279
Keywords
 Measurement Noise
 Extend Kalman Filter
 Process Noise
 Phantom Data
 Ultrasound Elastography