Implementation and clinical application of a deformation method for fast simulation of biological tissue formed by fibers and fluid
- Ana Gabriella de Oliveira Sardinha^{1},
- Ceres Nunes de Resende Oyama^{2},
- Armando de Mendonça Maroja^{1} and
- Ivan F. Costa^{1}Email author
https://doi.org/10.1186/s13029-016-0054-x
© Sardinha et al. 2016
Received: 10 December 2014
Accepted: 31 March 2016
Published: 15 April 2016
Abstract
Background
The aim of this paper is to provide a general discussion, algorithm, and actual working programs of the deformation method for fast simulation of biological tissue formed by fibers and fluid. In order to demonstrate the benefit of the clinical applications software, we successfully used our computational program to deform a 3D breast image acquired from patients, using a 3D scanner, in a real hospital environment.
Results
The method implements a quasi-static solution for elastic global deformations of objects. Each pair of vertices of the surface is connected and defines an elastic fiber. The set of all the elastic fibers defines a mesh of smaller size than the volumetric meshes, allowing for simulation of complex objects with less computational effort. The behavior similar to the stress tensor is obtained by the volume conservation equation that mixes the 3D coordinates. Step by step, we show the computational implementation of this approach.
Conclusions
As an example, a 2D rectangle formed by only 4 vertices is solved and, for this simple geometry, all intermediate results are shown. On the other hand, actual implementations of these ideas in the form of working computer routines are provided for general 3D objects, including a clinical application.
Keywords
Background
A realistic and fast soft tissue model must be used effectively in various medical applications, such as planning surgery procedures, image-guided surgery, image registration, diagnosis, biomechanical data refinement, and for training physicians [1].
There are many deformable physics-based methods used for surgical simulation. Meier [2] and Badosgan [3] report on some of the following methods: boundary element method, tensor-mass model, point-associated finite-field approach, and the most widely used finite element method and mass-spring model. They highlight the advantages and drawbacks of each method regarding the level of accuracy, computational load, difficulties and needs during implementation, numerical stability, etc.
The new method presented by Costa [1] runs in real time and can simulate biological soft tissues formed by fluid and a dense network of deformable fibers connecting surface vertices. The deformed state of the mesh is computed equating internal forces, due to fluid pressure and fiber tension, with external forces acting in an area associated with each superficial vertex. The fibrous tissue is similar to mass-spring simulation. However, unlike mass-spring simulation, point masses are not necessary. The mass is distributed in the entire object through fluid density. By enforcing the volume conservation, a behavior similar to stress tensor is obtained, reminiscent of the finite element method. As a result, this method provides some interesting outcomes. It is suited to anisotropic elasticity and non-linear stress–strain relationship. The results are accurate independent of mesh discretization. Only a few material parameters are needed. On the other hand, an important limitation of this method is that it is only valid for objects filled with fluids. Another drawback is that it has no dynamic behavior. Therefore, movements such as waves and vibrations on the object surface cannot be simulated. However, the quasi-static approach has the advantage of being numerically stable.
For problems concerning long-range connections, like the behavior produced by fibers, this approach defines a mesh of smaller size than volumetric meshes, allowing simulation of complex objects with less computational effort. Moreover, user interaction is minimized by dismissing the tedious and time-consuming need for mesh generation [4] and using a fully automated node-fiber-node model instead. On the other hand, volumetric meshes with only local connections can produce a sparse matrix, in which case the numeric solution would be asymptotically faster than this method.
The central theme of this paper is to provide a step by step discussion, algorithmic, and actual implementations in the form of working computer routines of the ideas expressed in Costa [1] for fast simulation of biological tissue. In addition, the software for clinical applications is shown, for the first time, by using the program to deform a 3D breast image of a real patient through the use of a 3D scanner in a hospital environment.
Methods
Surface geometry definition
The surface geometry of the objects to be simulated can be derived from data scanned from real data or from intermediate models. In general, this data is given as a set of surface points (vertex) positions.
For a deformable object, a vector \( {\overrightarrow{s}}_i \) represents the vertices’ 3D positions, where i enumerates each vertex between 1 and the total number of surface points N. For the computational code, the x-, y- and z-components of the vector \( {\overrightarrow{s}}_i \) are stored in (X[i], Y[i], Z[i]) and N is stored in the variable NVertex.
The deformation for each vertex can be specified by a displacement vector field \( {\overrightarrow{u}}_i \) and its x-, y- and z-components are stored in unique column matrix (u[i], u[i + NVertex], u[i + 2 * NVertex]).
A triangular mesh comprises a set of N′ triangles (in three dimensions) connected by their common vertices and defines the surface shape of an object in 3D solid modeling. To describe each triangle we enumerate three vertices V[1][i], V[2][i], and V[3][i] which must be connected to form a face, similar to the definition of the WaveFront Object (.obj) File Format [5]. In this case, the index i enumerates each triangle between 1 and the total number of surface triangular faces N’ which in the computational code is stored as NFaces. Note that the V[1][i], V[2][i], and V[3][i] have values between 1 and the total number of surface points N.
Computer routines for faces and vertex areas
The flat nature of triangles makes it simple to determine their normal vector, a three-dimensional vector \( {\overrightarrow{A}}_i \) perpendicular to the i-th triangle’s surface. Vector \( {\overrightarrow{A}}_i \) can be obtained by calculating the cross product between the vectors that form two edges of the triangle divided by two. Thus the modulus of \( {\overrightarrow{A}}_i \) is the area of the triangle. The direction of vector \( {\overrightarrow{A}}_i \) can be chosen to point outside of the object.
In our method, the deformation of 3D objects is performed through the displacement of its vertices. An area vector is assigned to each vertex, in order to determine the force acting over its surface. This area is defined from the area vectors of the triangles (faces) in the vicinity of the vertex.
Surface area of objects which the circumscribed sphere (one that touches the polygon at all vertices) has a unitary radius (R=1). The number of vertices N and the number of triangular faces N’ are also shown
Object | N | N’ | Analytical Equation | Analytical Result | \( {\displaystyle \sum_{i=1}^N\left|{\overrightarrow{S}}_i\right|} \) | \( {\displaystyle \sum_{i=1}^{N^{\prime }}\left|{\overrightarrow{A}}_i\right|} \) |
---|---|---|---|---|---|---|
Octahedron | 6 | 8 | \( 4\sqrt{3}{R}^2 \) | 6.92820 | 6.92820 | 6.92820 |
Icosahedron | 12 | 20 | \( {\scriptscriptstyle \frac{40\sqrt{3}}{5+\sqrt{5}}}{R}^2 \) | 9.57454 | 9.57454 | 9.57454 |
≈ Sphere | 642 | 1280 | ≈ 4πR ^{2} | 12.5664 | 12.5037 | 12.5065 |
Deformation routine
Note that Eq. 5 couples the dislocations (and hence, forces) in perpendicular directions (x, y and z). This coupling creates an effect somewhat similar to the stress and strain tensors in standard elastic theory.
The force due to one fiber connecting vertex i and j obeys Hook’s Law and is proportional to the fiber’s area (S _{ i } +S _{ j } )/2 and inverse to its length \( \left|{\overrightarrow{s}}_i-{\overrightarrow{s}}_j\right| \), where \( {\overrightarrow{s}}_i \) and \( {\overrightarrow{s}}_j \) are the positions of vertices i and j. The force due to all fibers \( {\overrightarrow{F}}_i^{fibers} \) is obtained by connecting each vertex i to other vertices by a set of N-1 elastic fibers. These fibers are shown as green lines in the upper inset of Fig. 2 for a sphere.
These equations can be written as a problem of type A . x = B where B is the column vector defined by the right hand side of Eq. 4 (\( {\overrightarrow{F}}_i^{contact}+\rho {h}_ia{\overrightarrow{S}}_i \)) and one extra element equal to zero that imposes the conservation of volume (right hand side of Eq. 5). The left hand sides of Eqs. 4 and 5 define the square matrix A.
The resulting vector x gives all the displacements and the pressure variation. An example of matrices x, A and B can be seen in Eq. 8.
Border conditions
The degenerate first mode, corresponding to the zero eigenvalue, represents a rigid body translation because, although moving, each vertex is stationary relative to the other. Of course, the presence or absence of this degenerate mode will not influence the purely deformable characteristics of the system. This degeneracy can be removed using boundary condition.
For boundary condition some vertices can be considered fixed to an external support. To achieve this condition, vertex i must not move, or moves a negligible amount compared to the movement of the others vertices. The resulting movement x _{ i } of the vertex is controlled by the size of element A _{ ii }, because each element of the product of matrices A and x is naturally expressed as the sum of N products A _{ ij } x _{ j }. In order to impose the border conditions, i.e., make x _{ i } <<x _{ j }, the element A _{ ii } should be done much bigger than the other elements A _{ ij } in line i of matrix A. An example of this procedure can be seen in section "An example for a rectangle".
Computer routines for general 3D objects
We need to create a matrix A[1..3N+1] [1..3N+1] as the input matrix of equation A . x = B. A large number of elements in matrix A in general vanish. Thus we initialize by writing zeros in all elements of this matrix.
For a deformable object the triangles’ positions and shapes change. Therefore, we need to recalculate the perpendicular of each triangle and vertex at each calculation step.
Concave objects
Note that all fibers effectively exist for a convex object. However for a concave object, some fibers connect pairs of vertices beyond the surface of the object, so these fibers do not actually exist. In order to consider this situation, forces due to fibers need to be removed when they pass through the surface. Then the number N is equal to the number of vertices only for the convex object. For concave object, N must be equal to the number of effective connections.
Two tests must be done in order to detect the fibers that pass through the surface of the object. First, we need to test if the fiber starting at vertex i goes inside or outside the object. Then, we calculated the scalar product for all k triangles neighboring the vertex i: \( {\overrightarrow{A}}_k\cdot \left({\overrightarrow{S}}_i-{\overrightarrow{S}}_j\right) \). If there is any negative result, the fiber connecting the vertices i and j goes outside the object and the matrix element A _{ ij } must be set equal to zero.
A second test is needed because if a fiber intersects any triangle, it goes outside the object. Then we used the fast 3D line segment-triangle intersection test developed by Chirkov [6]. If a fiber that connects the vertices i and j intersects any triangle, the matrix element A _{ ij } also must vanish. In this link [7] there is an executable beta version of our software with this and other functionalities.
The matrix solution
To simulate a non-linear stress–strain relation the value of Y _{ ij } can be defined as a function of the displacement \( {\overrightarrow{u}}_i-{\overrightarrow{u}}_j \) in the last step. Shortly after each increment of deformation, the element geometry and stiffness value is updated to model the nonlinear behavior of the material. The variation of external force on each interaction should be chosen to obtain the desired precision or an acceptable calculation time. Therefore, in this case, matrix A must be recalculated at each step.
But for situations where the area \( {\overrightarrow{S}}_i \) and the distance of each vertex \( \left|{\overrightarrow{s}}_i-{\overrightarrow{s}}_j\right| \) varies in a negligible amount on the left hand side of Eqs. 4 and 5 and the elasticity is linear (surface tension γ_{ik} and Young’s Modulus Y _{ ij } are constants independent of the deformation), matrix A is constant and can be operated once. For example, LU decomposition from the book Numerical Recipes in C [8] was used to solve the problem A.x=B. If matrix A is constant, the routine given in this book allows LU decomposition result to be left in place for successive calls with different right-hand sides B to achieve a greatly reduced calculation time.
Results and discussions
An example for a rectangle
For simplicity, the 2D Young’s Modulus is set to be a constant Y _{ ij } =360 and no acceleration effect is taken into account. The boundary conditions are chosen to be as least restrictive as possible: the lower left corner is fixed, the lower right corner is constrained to move only horizontally. In order to obtain these boundary conditions, we make the matrix elements A _{ ii } large, equal to 10 ^{ 9 }, for these vertices in these directions.
Solving this matrix equation for u we obtain the dislocations that should be added to the vertices position to get the new deformed shape. The result for the position of the right corners is s _{ x } =11 and for the upper corners s _{ y } =25, as shown in Fig. 3. Other examples for 3D objects are shown above.
An example for a sphere
Let us now consider that a sphere with a radius of one meter. The border conditions impose that all upper vertices, defining a spherical cap of height 0.2, cannot move. The parameters are defined as: the surface tension γ _{ ik } and Young’s Modulus Y _{ ij } are adjusted to be equal to 3 with units in kNewtons and meters; the acceleration due to gravity is acting down and increases at each iteration 1 m/s ^{ 2 } in a total of ten steps until the final value g=10 m/s ^{ 2 } is reached; the matrix A is constant; the number of vertices is N=162; and the density of the fluid is ρ = 1000 kg/m ^{3}. To show the effect of contact forces \( {\overrightarrow{F}}^{contact} \) , at the final step a force \( \overrightarrow{F}=-300\left(\widehat{i}+\widehat{j}+\widehat{k}\right) \) is applied to the vertex number 1.
The ideas given by Wright [10] were used to draw the objects and shadows. In addition, the routine to make spheres was taken from Shreiner [11].
An example for a real in vivo breast
The tissues were simulated as isotropic and linear, allowing us to do LU decomposition of the matrix A only once before the first iteration. Thus, for each interaction the right hand of Eq. 4 (vector B) are updated and new values for \( {h}_i,\kern0.49em a\kern0.37em \mathrm{and}\ {\overrightarrow{S}}_i \) must be incorporated at the next calculation step. We choose this procedure because, according to Costa [1], the resultant deformation in this case is similar to the deformation of the slower procedure that repeats LU decomposition at each interaction (updating \( {\overrightarrow{S}}_i \) and \( \left|{\overrightarrow{r}}_i-{\overrightarrow{r}}_j\right| \)).
Nevertheless, for nonlinear approaches the matrix A changes for each step and it must be recalculated for each iteration, imposing the necessity of the slower method. An example of a nonlinear elasticity is shown in Costa [1].
Best accuracy is expected if small steps are used (great number of interactions), since in this case the variation of \( {h}_i,\kern0.49em a\kern0.37em \mathrm{and}\ {\overrightarrow{S}}_i \) is smaller and the outcome get closer to the continuous (analytical) result.
We chose the surface tension γ _{ ik } = 48 MN/m, the Young’s Modulus Y _{ ij } = 48 kPa and the density of the fluid ρ = 1000 kg/m ^{3}. The boundary conditions in the case of breast deformation are the chest region that restrains the tissue movement in any direction.
For a continuous (analytical) approach, no difference on uncompressed and decompressed position are expected, and the result of Eq. 9 must be zero. For a numerical (discrete) approach differences on the initial and final position are expected due to errors introduced in the calculation of matrix B in each step. However, in our algorithm the results of Eq. 9 are small; numerically 3 % for 1 interaction and decreased to about 0.3 % and 0.03 % for 10 and 100 interactions, respectively. And the method evaluation results were similar for all four patients. These results demonstrate a good algorithmic behavior for the complete compression and decompression process of the breasts.
Conclusions
This work presents the implementation of a new general approach for modeling soft tissue compression process. The method was successfully applied to compress a 3D breast model from real patients.
The fast simulation of deformation of biomaterial using this algorithm could provide more realistic images that could serve for educational, clinical application, and research purposes. The latter include investigations of different breast imaging techniques involving compressed and uncompressed breasts.
Declarations
Acknowledgments
The authors gratefully acknowledge the University of Brasilia and the financial support from the Brazilian science agency (CNPq) through research project 474831/2012-4.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
- Costa IF. A novel deformation method for fast simulation of biological tissue formed by fibers and fluid. Med Image Anal. 2012;16:1038–46.View ArticlePubMedGoogle Scholar
- Meier U, López O, Monserrat C, Juan MC, Alcañiz M. Real-time deformable models for surgery simulation: a survey. Comput Methods Prog Biomed. 2005;77:183–97.View ArticleGoogle Scholar
- Basdogan C, Sedef M, Harders M, Wesarg S. VR-based simulators for training in minimally invasive surgery. IEEE Comput Graph Appl. 2007;27:54–66.View ArticlePubMedGoogle Scholar
- Filho PJ, Sousa EAC. Reconstruction and mesh generation in three-dimensional biomechanical structures for finite elements analysis. Braz J Biom Eng. 2009;25(1):15–20.Google Scholar
- McHenry K, Peter B. An overview of 3d data content, file formats and viewers, National Center for Supercomputing Applications. 2008. p. 1205.Google Scholar
- Chirkov N. Fast 3D Line Segment—Triangle Intersection Test. J Graph Tool. 2005;10(3):13–8.View ArticleGoogle Scholar
- Costa IF. Authors’ own website. https://www.sites.google.com/site/unbivan/download. Accessed 12 April 2016.
- Press WH, et al. Numerical Recipes in C: the art scientific computing, 2th edn. Cambridge: Cambridge University Press; 1992. (Chapter 2.3 - LU Decomposition and Its Applications)Google Scholar
- Costa IF. Authors’ own website. https://sites.google.com/site/unbivan/FFM.cpp. Accessed 12 April 2016.
- Wright R, Sweet M. OpenGL SuperBible. 2nd ed. Essex: Pearson Education; 1999 (Chapter 9 – Lighting and Lamps).Google Scholar
- Shreiner D. OpenGL programming guide: the official guide to learning OpenGL, versions 3.0 and 3.1. 7th ed. Boston: Addison-Wesley; 2009. Chapter 2 - State Management and Drawing Geometric Objects.Google Scholar