Abstract
Characterizing heterogeneous elastic property distribution of soft tissues is of great importance in disease detection. In this paper, we investigate an inverse approach to map the heterogeneous material property distribution of soft solids using harmonic motion data. To examine the feasibility of this approach, a number of numerical examples are presented. We observe that the shear modulus distribution is recovered well using harmonic motion measurements. Compared to the static inverse approach, the proposed dynamic inverse method improves the quality of the recovered shear modulus distribution significantly. We also study the influence of the uncertainty in the driving frequency on the reconstruction results and observe that the influence is not very significant in recovering the shape of the inclusion. The proposed inverse algorithm has potential to be a promising tool to diagnose diseases in clinical medicine.
1. Introduction
Mechanical signal or response has a long history of being used for health assessment and disease detection. For instance, more than 2000 years ago, Chinese developed a pulse diagnosis where physicians assessed patients’ health conditions based on their wrist-pulse [1, 2] since alteration of the frequency, amplitude, etc. of the wrist-pulse was highly correlated to the pathologic changes. With the fast development of mechanical techniques, a large number of mechanical based medical devices including ultrasonography have been invented and widely used in healthcare institutions. In particular, the development of imaging modalities [3, 4] provides us with the availability to measure the deformation and motion of tissues and organs inside our body. With the full-field measured data of the region of interest, we are capable of mapping the heterogeneous mechanical property distribution of soft tissues. Since many diseases such as cancerous diseases [5, 6] and cardiovascular diseases [7, 8] induce the alteration of mechanical properties of associated tissues, mapping heterogeneous mechanical behavior might be useful for the disease detection.
With the accurate imaging data, identifying nonhomogeneous material property distribution requires solving an inverse problem. Due to its ill-posed nature, it is very difficult to solve the inverse problem. There are many methods to solve the inverse problem such as direct approaches [9–12] and statistical approaches [13]. A prevalent approach is to regularize the problem and pose the inverse problem to be a constrained minimization problem [11, 14–16]. To improve the computational efficiency of solving the inverse problem using an iterative inverse algorithm, the adjoint method has been proposed. This approach has been widely applied to map both linear [17, 18] and nonlinear [19, 20] elastic property distributions of soft tissues utilizing measured displacement fields in the quasi-static case.
In this paper, we will generalize the regularized inverse approach into dynamic cases. We assume the tissue or biological organ is subjected to harmonic motion, solve the inverse problem in the frequency domain, and map the linear elastic property distribution using harmonic motion data. This paper is organized as follows: In Methods, we will discuss the mathematical details of the proposed approach and a number of numerical examples will be presented in Results. We will discuss the proposed method and associated results in Discussion and paper closure will be in Conclusions.
2. Methods
The wave equation for harmonic motion in the frequency domain is written aswhere and denote the displacement vector and stress tensor, respectively. and represent the mass density and angular frequency, respectively. In addition, the vectors and are the boundary conditions at. In this paper, the solid is assumed to be an incompressible, linear elastic and in the state of plane stress; thus the stress-strain relation iswhere is shear modulus and is the strain. For a known shear modulus distribution, the displacement field can be acquired by solving a forward problem using the finite element method (FEM), leading to the following discretized equations:where K and M are the stiffness matrix and mass matrix, respectively. Meanwhile, and f are the displacement and force vectors, respectively. Since FEM has been widely used for solving equations of motion, for brevity, we do not intend to discuss here.
The inverse problem is solved by an optimization approach where an objective function is minimized in the L2 norm:where and are the nodal measured and computed displacements, respectively. The computed displacement field is obtained by solving the forward problem at the current estimated shear modulus distribution. The shape function W represents the approximation from the continuous displacement field to the associated discretized field. The second term in (4) is the regularization term. In this paper, we employ the total variation diminishing (TVD) regularization term (, where c is a small constant and is set to 10−2 in order to avoid the singularity when computing the derivative of the regularization term with respect to shear moduli). is the regularization factor to control the contribution of the regularization term to the objective function. A smaller leads to strong distortion of reconstruction, while a larger value will oversmooth the final results. In this paper, the optimal regularization factor is visually determined. To be specific, we start with a very large regularization, keep solving the inverse problem with a decreasing regularization factor, and then observe the shear modulus reconstructions. This optimal regularization factor will be determined when the reconstruction of the background starts oscillation. The same strategy has also been utilized in [17–21].
The inverse problem is solved by a quasi-Newton approach, the L-BFGS (limited-Broyden–Fletcher–Goldfarb–Shanno) method, which requires the objective function value and its spatial gradient with respect to shear moduli. The gradient of the objective function can be calculated as follows:where j represents the global node number. denotes the inner product. Differentiating (3) with respect to the nodal shear modulus yieldsSubstituting (6) into (5) leads toThis straightforward approach to evaluate the gradient is computationally intensive; thus the adjoint method is adopted to calculate the gradient in a highly efficient way [21]. More specifically, if we rewrite (7) by taking advantage of the definition of a transpose, the following equation can be obtained:Accordingly, the adjoint equation can be acquired:If we solve the adjoint equation for the vector D, the gradient can be expressed asThus, we merely need to solve two forward problems at every minimization iteration using the adjoint method. The inverse solver will terminate when one of the two following stop criteria is satisfied: the difference between the objective function values at the current and last minimization iterations is smaller than the machine precision; the norm of the gradient of the objective function with respect to shear moduli is smaller than the machine precision.
The in-house inverse algorithm is implemented by Fortran and parallelized by openMP. For the L-BFGS algorithm, we adopt an open source L-BFGS subroutine developed by [22, 23]. In this work, we will primarily test the inverse algorithms; thus measured data obtained by simulation will be employed in numerical examples, assuming the material property distribution is known. We then use the simulated data to solve the inverse problem and compare the reconstructed shear modulus distribution with the target values. To mimic real data, we add up to 10% white Gaussian noise into the simulated data, and the noise level is defined aswhere is the total number of displacement data. and are the measured displacement and noise-free displacement, respectively.
3. Results
In this section, we will present numerical examples where the geometric model is shown in Figure 1(a). More specifically, a circular inclusion with a shear modulus value of 500Pa is embedded in the 1×1cm2 square background with a shear modulus value of 100Pa. The radius of the inclusion is 0.1cm. The square model is discretized by 3600 bilinear elements. In respect of the boundary conditions, we fully fix the bottom edge and apply 1% shear deformation on the top edge. When solving the inverse problem, the initial guess of shear modulus distribution is homogeneous throughout the problem domain and the initial shear modulus value is 10Pa. In addition, we restrict the search domain of the shear modulus of every node to interval 10,3000Pa.
(a)
(b)
(c)
(d)
As seen in Figure 1, when the noise level is low (3%), the shear modulus distribution is recovered with very good quality, since both the value and shape of the shear modulus of the inclusion are close to the target. With the increase of noise level, the reconstruction becomes worse (see Figures 1(b)–1(d)). In particular, when the noise level reaches 10%, the mapped inclusion is distorted and a strong artifact is observed on the background in Figure 1(d). We also investigate the sensitivity of the regularization factors to the reconstructed shear modulus distribution (see Figure 2). We observe that, for a very small regularization factor ( in Figure 2(b)), the reconstructed shear modulus distribution oscillates significantly. For a very large regularization factor ( in Figure 2(d)), the shear modulus distribution of the background is very smooth. Meanwhile, the inclusion can be recovered successfully but becomes much larger than the target. Thus, the optimal regularization factor should be selected between them as shown in Figure 2(c).
(a)
(b)
(c)
(d)
For a relative low driving frequency (Figure 3), we find that the inclusion is also recovered well in the case of low noise level (3% noise). However, compared to the static case (see Figure 1(b)), both the shape and value of the inclusion are slightly worse mapped. Besides, though increasing noise level will deteriorate the mapped shear modulus distribution, the quality of reconstructed results in this dynamic case is slightly better than the static case (compare Figures 1(c) and 1(d) and Figures 3(c) and 3(d)). For higher driving frequencies of 20Hz and 40Hz (Figures 4 and 5), we observe that the inclusion is well recovered for noise levels up to 5%. Even with 10% noise level, the inclusion is recovered well without too many artifacts. Comparing the reconstructed results for varying driving frequencies, we also observe that the driving frequency of 40Hz yields the best reconstruction results.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
We also study the effect of the uncertainty in the driving frequency on the reconstructed results. In Figures 6 and 7, we also add noise to the driving frequency in solving the inverse problem. We discover that the decrease of driving frequency reduces the value of the mapped inclusion (see Figure 6) but also reduces the artifact of the background simultaneously. The shape of the mapped inclusion seems to remain the same level for different driving frequencies. When we raise the driving frequency (Figure 7), it is clear that the value of the shear modulus of the mapped inclusion rises but the background experiences a stronger oscillation than that using the exact driving frequency.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
4. Discussion
This paper presents the regularized inverse approach to map the heterogeneous elastic property distribution of the soft solids using harmonic motion data. It took 5000-10000 iterations for convergence to solve the inverse problem in this paper. We compared the reconstructed results when the driving frequencies are 0Hz (static case), 2Hz, 20Hz, and 40Hz, respectively. We also varied the initial guess of the shear modulus from 10Pa to 1000Pa and acquired very similar reconstructed shear modulus distribution. We notice that using static data is capable of mapping the elastic property distribution well with a low noise level. However, for higher noise level, the dynamic data yields better reconstruction results. Thereby, this is one of the advantages of utilizing dynamic measurement. Another merit of using harmonic motion data is that we can quantitatively determine the shear modulus distribution merely using displacement measurements. This is impossible for the static case owing to the homogeneity of the equilibrium equations. Thereby, we must know nonzero force or traction information, or shear modulus values in a certain subregion for the static case. Otherwise, the shear modulus distribution can only be determined relatively up to a multiplicative factor. This has been well studied in [18]. We also learn that the uncertainty in the driving frequency might not necessarily reduce the quality of reconstructed results. However, the reason behind this remains an open question. Though we merely test the plane stress case in this paper, the proposed approach can be easily generalized to 2D plane strain and 3D cases. Besides, the experimental data should be utilized to test the feasibility of the proposed method. In this paper, we reveal that using harmonic data is capable of yielding higher quality of reconstruction even with very high noise level. Further, the error in driving frequency does not decrease the reconstruction quality significantly. Therefore, this analysis is of great significance in applying the approach to the practical cases.
5. Conclusions
In this paper, we study the feasibility of characterizing the nonhomogeneous elastic modulus distribution utilizing full-field harmonic motion data. We test several numerical cases and observe that this approach is capable of mapping the elastic property distribution well even with high noise levels. We also investigate how the uncertainty in the driving frequency affects the reconstruction of the shear modulus distribution. We realize that the uncertainty in the driving frequency might not reduce the quality of reconstruction results.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
Peng Yu thanks the support from China Scholarship Council (201406150089).