|Home | About | Journals | Submit | Contact Us | Français|
Material property has great importance in deformable body simulation and medical robotics. The elasticity parameters, such as Young’s modulus of the deformable bodies, are important to make realistic animations. Further in medical applications the (recovered) elasticity parameters can assist surgeons to perform better pre-op surgical planning and enable medical robots to carry out personalized surgical procedures. Previous elasticity parameters estimation methods are limited to recover one elasticity parameter of one deformable body at a time. In this paper, we propose a novel elasticity parameter estimation algorithm that can recover the elasticity parameters of multiple deformable bodies or multiple regions of one deformable body simultaneously from (at least two sets of) images. We validate our algorithm with both synthetic test cases and real patient CT images.
Material properties are important in depicting the characteristics of virtual objects for realistic computer animations of soft bodies. In addition, virtual surgical simulation has also been increasingly used for rapid prototyping of clinical devices, pre-operation planning of medical procedures, virtual training exercises for surgeons and medical personnel, etc. And, tissue elasticity properties are important parameters for developing accurate and predictive surgical simulation. Futhermore, to compute desired and accurate force feedback for haptic display requires knowledge about the deformation of soft tissues and organs, which are characterized by patient-specific elastic parameters for different tissues and organs.
Elastography  was first proposed to determine the elasticity properties by measuring the deformation of the tissue due to the application of the known external forces. The known external forces are the boundary condition needed to recover the exact elasticity parameter. Originally, the deformations were measured using ultrasound imagery , but such techniques produced coarse, two-dimensional representations of the moving tissue. More sophisticated imaging techniques, such as magnetic resonance imaging (MRI) [3, 4]and computed tomography (CT), produce three-dimensional images of the deforming tissues, allowing a more accurate measurement of displacement.
Toward realizing the concept of 3D physiological humans, we propose perhaps one of the first elasticity parameter estimation algorithm for multiple, heterogeneous deformable bodies simultaneously using medical images1. Our approach is based on a multi-dimensional optimization method that iteratively performs deformable body simulation using a finite element method on reconstructed organ models with the continuously refined, estimated elasticity parameters. The geometric models of organs are reconstructed based on low-resolution CT images. Our objective function measures the sum of the distance between the nodes of the organ surface. In contrast to elastography methods [5, 6, 7], the only information we need is the displacement of the nodes of the organ surface. We do not need every pixel-wise displacement vector, thus no extra procedures need to be performed on the patient. Two sets of (medical) images are sufficient to recover the elasticity parameters using our method. Therefore, our method can be widely applicable to different imaging technology. It can be used for animatino of soft bodies (see supplementary video) and possibly for cancer staging using only low-resolution CT images.
In computer graphics, extensive research have been done for deformable body simulation [8, 9, 10, 11]. Researchers have also been studying methods for designing simulation with material properties [12, 13, 14, 15] and for realistic bio-structure simulation .
In medical applications, there are mainly two kinds of soft tissue elasticity properties estimation method : invasive and non-invasive techniques. The invasive methods rely on a device to measure the displacement and force response [18, 19, 20]. These methods take organ samples either out of the human or animal bodies and perform the experiment in-vitro (out side the body) or do the procedure in-situ (inside the body). The collected data are then used to solve the inverse problem, which is to recover the elasticity properties, by constructing a polynomial interpolation  or by using a finite element model [22, 23, 17].
The non-invasive methods mostly base on image analysis techniques to measure the displacement. In the 1980s, several methods were proposed to measure the motion of the soft tissue, such as the one proposed by Wilson and Robinson  using radio frequency M-mode signals and the one proposed by Birnholz and Farrell  using ultrasound B-scans. Researchers have also used medical image analysis on 2D ultrasound and/or MRI images to estimate the elastic parameters of soft tissue [26, 27, 28].
In the area of elastography, researchers [5, 6, 7] proposed algorithms based on the distance between two medical images. By solving the least square problem, the elasticity parameters are recovered. Van Houten et al.  used elastography methods to estimate the Young’s modulus distribution of a 2D area, then later extended to solve three-dimensional elastic parameter distribution using MRI . These methods need high-resolution displacement fields to recover the elasticity parameter , where the displacement field is typically obtained through an external device using a vibration actuation mechanism. For organs that are located deep-seated inside a human body, the vibrator may need to be placed inside the organ , making the procedure much more complex and possibly uncomfortable for the patient.
Other than distance field based methods, there are also other measurement algorithms. The modality-independent elastography (MIE) method  measures the elasticity parameters by maximizing the image similarity based on a number of landmarks. However, this technique does not apply to all the soft tissues, as landmarks cannot always be found in some of the organs such as prostate. Statistical and machine learning algorithms have also been used to classify soft tissues and estimate the parameters using multi-spectral MR images .
Although the existing elastography methods can provide an estimation of the elasticity parameter distribution, they require high-resolution magnetic resonance medical images and a device to measure the external force exerted on the soft tissues, which is not always possible or practical. Lee et al. proposed the first model to estimate the Young’s modulus based on low-resolution CT images and no external force is needed to set the boundary condition . By optimizing the distance between the deformed and the reference surface meshes, the elastic parameter is estimated. However, this method can only recover the elasticity parameter of one organ. In contrast, our work can recover the elasticity paramenters of multiple, heterogeneous soft bodies simultaneously using a multi-dimensional optimization method.
We propose a novel method to automatically estimate elasticity parameters of multiple organs using only images, such as those from computed tomography (CT) imaging. Our approach is based on multi-dimensional optimization and simulation of multiple deformable bodies using finite element methods. For each optimization iteration, a finite element method is used to assess the deformation of the organs. The objective function is based on distance between the initial, reference mesh and the deformed surface mesh. This objective function is evaluated and its gradient is used in the multi-dimensional optimization algorithm to search for the optimized elasticity parameters that minimize the value of the objective function. The flow chart of the optimization process is shown in Figure 1.
A three-dimensional reconstructed model of the organs and the signed distance map of the deformed organs are the input of our algorithm. The three-dimensional model is reconstructed from the segmented CT images using ITK-SNAP . After each optimization iteration, the elasticity parameters are updated and used by the finite element model. The simulation-based parameter estimation algorithm is terminated when the optimization converges, which usually takes only a few iterations.
The forward simulation in our method computes the displacement vector u using the elasticity parameters recovered in the inverse problem. The displacement vector u is then used in the inverse problem to evaluate the objective distance function (Eqn. 4). Our simulation framework uses a linear static finite element method  (Eqn. 1). The weak formulation for elasticity problem is given,
where u is the displacement vector and t is the boundary condition(traction act on the boundary Γ). For quasi-static deformation process the ü = 0. We can rewrite Eqn. 1 as,
with the first part of the equation as the internal force and the second part of the equation as the external force.
We use linear elastic material model. For isotropic linear elasticity, the stress strain relation is defined as,
where σ is the stress tensor, ε is the strain tensor and matrix D is defined by the material elasticity parameters. We use Young’s modulus E and Poisson’s ratio ν for describing material properties. Previous elastography methods use external forces as the boundary condition. Our method, in contrast does not need external forces. The initial boundary condition in our algorithm is the known displacement vector of the surface mesh. This boundary condition is applied only to the nodes of the surface mesh. When the three-dimensional model deformed, the force generated by the deformation will drive the simulation.
The inverse problem is the process of elasticity parameter estimation. Our method is based on the multi-dimensional optimization method. By solving the least square problem iteratively, we recover the elasticity parameter. We then use the distance between the initial, reference surface mesh and the deformed surface mesh to iteratively update the objective function.
As our simulation framework is based on the low-resolution CT images, only the displacement of the boundary of the soft tissue is known. Our objective function (Eqn. 4) is constructed using the sum of the distance between the nodes of initial, reference surface and that of the deformed surface. By minimizing the value of the objective function, we find the optimal elasticity parameters μ.
where m is the mth organ and N is the number of the organs that are in the simulation scene. d(ul + Δul, ur) is the bidirectional Hausdorff distance between deformed surface mesh Sm and the initial, reference mesh Sr. μm = Em in which Em is the Young’s modulus of the mth organ. Our method can be extended to optimize more than one parameters. We could also include the Poisson’s ratio into μm. The μ that minimizes the objective function is the optimized set of elasticity parameter.
We propose to use multi-dimensional optimization method to recover the elasticity parameters. Our three-dimensional model can have a large number of nodes, so a significant amount of memory would be needed to store the exact Hessian matrix for the Newton’s optimization method. Therefore, to solve the least square problem, we choose the Limited Memory Quasi-Newton’s method. Using this method, the approximation of the Hessian matrix is maintained instead of the exact Hessian matrix. For each step of this BFGS method ,
where Hk is the approximated Hessian matrix, xk is the variable to be optimized, Φk is the objective function value, and k denotes the kth optimization iteration and Φ is the gradient of the objective function. To compute the gradient, the partial derivative of Φ with respect to the elasticity parameter, the Young’s modulus of the mth organ.
We have implemented our algorithm and performed three sets of experiments to evaluate its accuracy under different conditions using both two synthetic sets of models with known parameters to validate the approach and a reconstructed set of organs from CT images to illustrate its robustness.
The first experiment is designed to test the accuracy of the algorithm, if the three organs sharing boundary. As the number of the organs increases, the problem becomes even more complicated.
We used 3-concentric spheres to build the test model in experiment I. In order to measure the elastic parameter of sphere 1, for the area between sphere 1 and sphere 2 and the area between sphere 2 and sphere 3, tetrahedralization is done within sphere 1, the area between sphere 1 and sphere 2 and the area between sphere 2 and sphere 3. The sliced view of the three-dimensional model is shown in Figure 3. The following table is generated when the three areas are all deformed by slightly less than 10%.
The result of this experiment is shown in Table 1. In this experiment we increase the number of ”organs” to test the accuracy of our algorithm. The result of this experiment is affected by both the fact that the spheres are sharing boundaries and the number of the spheres. Under these complicated conditions, the error rate of our algorithm is generally less than 5% (no more than 15% for the highest stiffness values) when the deformation is less than 10% within each of the three regions respectively.
The second experiment is designed to test the robustness of our algorithm with more scene consists of multiple, separate organs in contact with each other. The simulation scene includes multiple organs within a male’s pelvis area. The surface meshes of the prostate, bladder, and rectum were reconstructed from the patient’s CT images. These reference surface meshes were used to construct the tetrahedral mesh of the simulated scene. A slice view of the tetrahedral mesh is shown in Figure 2. In the tetrahedral mesh, the rectum was modeled hollow inside, while the prostate and the bladder were modeled as a continuum represented by tetrahedral elements. The prostate and the bladder are two organs that we use to recover the elasticity parameters.
The signed distance field within each organ was computed using the initial, reference surface mesh and the deformed surface mesh. The deformed surface mesh was generated based on the displacement of the nodes on the surface mesh. We used the initial displacement to set the initial forces as the boundary condition. Then the boundary condition was used to generate a displacement field, which was computed by applying the boundary condition to the three-dimensional model during each iteration of the optimization. The model was deformed using the current set of elasticity parameters. For the synthetic test case, we generate the deformed surface by using the set of “ground truth” parameters. We then run our algorithm on the resulting deformed surface to estimate the elasticity parameters and compare these recovered values with the ground-truth values as shown in Table 2.
We use patient specific medical images to reconstruct the organ models. The CT images was segmented using ITK-SNAP . After segmentation, we reconstruct the surface mesh of the prostate, bladder and rectum also using ITK-SNAP . The surface meshes are shown in Figure 6. Then, we used TetGen  to generate the tetrahedral mesh based on the surface mesh of the organs and a bounding box, with the rectum being hollow inside. We used the deformation of the rectum to set the boundary condition. The deformed prostate, bladder surface mesh are used to compute the updated, signed distance map.
We first did a search of the initial relative Young’s Modulus. The initial value affects the optimization result and several initial values are used to ensure convergence and determine the global minima. The experimental result is shown in Table 2. In our experiment, we initially assigned the same Young’s modulus to both organs, the prostate and the bladder. The three dimensional models are deformed based on the boundary condition we provided and the Young’s modulus we assigned to the two organs. The boundary condition was defined according to the two surface meshes, which were constructed using two CT images taken from the same patient in different time. From the table, we can observe that the relative errors of the bladder are generally smaller than those of the prostate. From , we know that the stiffness of the cancer tissue tend to be much more than 10% of the normal tissue, while the error of our algorithm is less than 10% in early all cases. Thus, given the relative errors of our algorithm, it can still detect cancer (see next). This result also shows the limitation of the linear elastic material model. The accuracy is affected by the amount of deformation of the organs.
We use multiple sets of CT images from each of eight patients (totaling 180 sets of images) in our real cancer correlation experimental study. The simulation scene which includes multiple organs within a male’s pelvis area is the same as the second experiment.
The experiment is designed to determine the correlation between the prostate cancer T stage and the elasticity of prostate and the surrounding organs. The T-stage is defined in TNM (Tumor, lymph Nodes, Metastasis) system  which is a common cancer staging system. The result of our experiment is shown in Figure 7. Using the box plot, we can observe that the mean of both the prostate’s and the bladder’s Young’s increase with the cancer staging (the resulting Young’s modulus is the average value of the entire organ). We further analyze the statistical signficance of this correlation between the T-stages and the elasticity of prostate and bladder. The resulting Pearson correlation coefficient for prostate’s Young’s modulus and T-stage is 0.658 and the p-value for two-tailed probability is 3.16 × 10−5, which indicate a strong correlation between the Young’s modulus of the prostate and the cancer T stage. The resulting Pearson correlation coefficient for bladder Young’s modulus and T-stage is 0.481 and the p-value for two tailed probability is 4.60 × 10−3. This result indicates that the Young’s modulus of the bladder increases with the stages of prostate cancer, but they are not strongly correlated. Our findings reconfirm the studies  that prostate cancer increases the probability of bladder cancer. As cancer continues to advance to higher stages, it spreads to neighboring tissues.
We have used a range of deformation that is larger than the normal tissue deformation to stress test our algorithm and to analyze the relationship between the degree of accuracy vs. the amount of deformation. This analysis helps us understand when non-linear models should be used. Additional sensitive analysis with respect to segmentation, simulation resolution (i.e. the size of mesh), and use of nonlinear FEM model can also be performed to provide additional information to the users. We ‘jumpstart’ the iterative optimization process with some range of default values and the algorithm usually converges quickly within less than 10 iterations in practice. Our implementation currently addressed the possibility of multiple solutions by using multiple (3–5) initial values sampled over a wide range (50–300) of possible values (say 50, 150, 250) and use multiple sets of the image data from few different days to compute the average values, after eliminating possible ’outliner value(s)’. With this approach, our algorithm is able to find the elasticity parameters that are very close to the ’ground truth’ values in practice.
In this paper, we presented a novel multi-body elasticity parameter estimation method using the low-resolution CT images. As our method do not require any external forces to be measured and only the deformation of the organ surface is needed, it can be applied to organs that locate deepseated in the human body. There are limitations, however. The amount of deformation of the soft tissue can affect the accuracy of the algorithm. The larger the deformation the higher the relative error is from the estimation. This is due to the fact that we have adopted a linear, static finite element method and linear elastic material model. Linear models are generally considered accurate and sufficient when the deformation is small and within a certain range where linearity assumption is applicable. Our experimental results support this observation. The linear elastic material model is not suitable for the simulation of human organs when they undergo large deformation. In the future, we plan to adopt a more complex, non-linear elastic material model for soft tissue simulation, such as Mooney Rivlin model. The accuracy may likely be higher when the amount of deformation is significant, though we expect the computation cost to increase as well.
The algorithm we proposed in this paper is based on a multi-dimensional optimization method, which can also be used to estimate multiple elasticity parameters of a single organ of multiple, heterogeneous tissue properties for different regions of a (human) body. Because of the importance of Young’s modulus in noninvasive cancer detection, we choose to estimate this parameter for multiple organs simultaneously. However, Poisson’s ratio has also been suggested as a significant indicator for breast cancer. Therefore, in the future we plan to further study the accuracy of multi-dimensional optimization method and hope to use it to estimate multiple elasticity parameters of a single organ for more accurate cancer screening and grading.
This research is supported in part by the Joint NSF and NIH Smart and Connected Health Program, NIH#R01EB020426–01. We would like to thank Drs. Ronald Chen and Edward Chaney for CT images from their Lab, and Dr. Huai-Ping Lee for data from his thesis in our experiments.
1In this paper, we use computed tomography (CT) images. But, the algorithm is also applicable to magnetic resonance imaging (MRI) images.
Shan Yang, University of North Carolina at Chapel Hill.
Ming Lin, University of North Carolina at Chapel Hill.