Abstract
In some simulations like virtual surgery, an accurate surface deformation method is needed. Many deformation methods focus on the whole model swing and twist. Few methods focus on surface deformation. For the surface deformation method, two necessary characteristics are needed: the accuracy and real-time performance. Some traditional methods, such as position-based dynamics (PBD) and mass-spring method (MSM), focus more on the real-time performance. Others like the finite element method (FEM) focus more on the accuracy. To balance these two characteristics, we propose a hybrid mesh deformation method for accurate surface deformation based on FEM and PBD. Firstly, we construct a hybrid mesh, which is composed of a coarse volume mesh and a fine surface mesh. Secondly, we implement FEM on coarse volume mesh and PBD on fine surface mesh, and the deformation of fine surface mesh is constrained by the displacement of the coarse volume mesh. Thirdly, we introduced a small incision process for our proposed method. Finally, we implemented our method on a simple deformation simulation and a small incision simulation. The result shows an accurate surface deformation performance by implementing our method. The incision effect shows the compatibility of our proposed method. In conclusion, our proposed method acquires a better trade-off between accuracy and real-time performance.
1. Introduction
The mesh deformation algorithm plays a very important role in computer simulation technique, such as computer animation and soft tissue deformation simulation in virtual reality. In some simulations like virtual surgery, the surface deformation plays an important part in realistic perform effect. Generally speaking, two major problems need to be solved for surface deformation algorithm: low time consumption and high accuracy. Although many methods focus on one of these two problems, there is a huge demand for high accuracy deformation under real-time performance. Therefore, a suitable deformation algorithm needs a trade-off between the accuracy requirement and the real-time performance.
The finite element method (FEM) is usually used for the numerical discrete model of constitutive equation in elastic mechanics [1]. The standard method is used to create physical elements according to each geometric element and then to calculate the displacement of relevant element nodes according to the applied load [2]. The boundary conditions need to be well treated. The FEM can build different kinds of physical elements according to corresponding geometric units, such as triangle element in 2D and tetrahedral element in 3D [3]. The advantage of FEM is that it is close to the real world since it utilizes the actual physical parameters to get output with high accuracy of deformation. It is also able to simulate the complex mesh model [4]. However, FEM costs more computational time. The real-time performance needs to be improved [5]. At present, FEM is widely used in mesh deformation simulations. Paulus et al. proposed a new re-mesh method and rapid FEM for soft tissue simulation [6]. Courtecuisse et al. demonstrated a real-time tissue simulation based on the implicit numerical integration method in nonlinear FEM [7]. Haouchine et al. proposed a method to enhance the internal structure of the liver mesh model in real time [8]. Tang and Wan proposed a constraint FEM for the interactive simulation for virtual surgery [9]. Kugelstadt et al. proposed a fast co-rotated FEM using operator splitting for whole mesh model swinging and twisting [10]. It acquires a good trade-off between accuracy and low time consumption but is not compatible with mesh dissection (model incision).
Position-based dynamics (PBD) has been widely used in recent years. Classical model methods such as FEM and MSM calculate the deformation according to Newton’s second law. Compared with these methods, PBD directly calculates the position change based on the solution of the quasi-static problem. Therefore, PBD has better real-time performance and stability [11]. However, PBD is generally not as accurate as classical model methods. In recent years, PBD has been widely used in deformation simulations. Kubiak et al. proposed a surgical line simulation method based on the PBD method of Müller et al. [12, 13]. In order to assist the robot-assisted rehearsal and planning of partial nephrectomy, Camara et al. proposed a real-time simulation platform that allows surgeons to quickly construct a biomechanical model [14]. Pan et al. proposed an interactive method of the hybrid soft tissue model based on extended PBD [15]. Liu et al. proposed a PBD method for tetrahedral element mesh model to simulate the deformation on soft tissue [16].
For the mass-spring method (MSM), it has many advantages such as easy operation, simple modeling, and low time consumption. However, MSM needs suitable constraints. Unsuitable constraints will lead to fault and divergence of deformation results. In MSM, it is important to determine the spring parameters which include stiffness and damping coefficient. These parameters determine the accuracy of MSM. Due to the random selection of these parameters, the accuracy of deformation results is not high. Two ways of parameter selection are data measurement and formula analysis. Data measurement estimates the parameters by measuring the deformation data of real biological tissues. The formula analytical method is based on the constitutive equation of the material to deduce the parameters [17]. The MSM also has a wide range of applications. Skornitzke et al. proposed a mitral valve MSM model [18]. Gao et al. simulated the nonlinear anisotropy of biological materials based on MSM [19]. Zhang et al. proposed a three-stage method based on MSM to simulate soft tissue [20].
In this paper, we propose a hybrid mesh deformation method for accurate surface deformation based on FEM and PBD. Our proposed method combined the advantages of FEM and PBD to satisfy the requirement of accuracy and real time. First, we construct the hybrid mesh which composes a fine surface mesh and a coarse volume mesh. Then, the FEM and PBD is utilized to simulate deformation on the volume mesh and surface mesh, respectively. Finally, the deformation of surface mesh is constrained on the volume mesh. We also introduced a small incision process which is compatible with our proposed method.
The reminder of this paper is organized as follows: Section 2 briefly explains purpose of our proposed method. Sections 3 and 4 describe the mathematical schemes including FEM and PBD. Section 5 introduces the construction of the hybrid mesh. Section 6 presents the algorithm flow of our proposed method. Section 7 introduces the small incision processing for our proposed method. Section 8 introduces the experiment environments, shows the rendering effects and results of deformation, and demonstrates the compatibility of small incision processing. In Section 9, we conclude our work.
2. Overview
In previous work, FEM is used for high accuracy in fine mesh deformation but usually has high time consumption. PBD has the characteristics of low time consumption and low accuracy. In our method, FEM is used to calculate the deformation of the coarse tetrahedral mesh, and the PBD is used to perform the deformation of the fine surface triangle mesh. The nodes of the volume mesh which connect to the surface mesh are control points (black points in Figure 1). These points constrain the surface mesh and make its deformation within a more accurate range. Meanwhile, since the FEM is applied to the coarse tetrahedral mesh and the surface mesh uses the PBD, the time consumption is not too high. Thus, the real-time requirements can be met.
3. Finite Element Method
Our method uses a set of tetrahedral elements to separate the 3D domain. Based on this spatial discretion, the node displacement interpolation method is used to approximate the continuous displacement field. Specifically, according to equation (1), using the node of elements, the shape function is interpolated to form the displacement field function of the mesh element.
is the shape function matrix of elements, and is the displacement of the mesh element node. We choose linear interpolation as the shape function of the tetrahedron discretion. Through the shape function, the relationship between the strain and the displacement of the element node is established:
is the strain matrix of elements. According to the principle of virtual work, the relationship between nodal displacement and external nodal force of elements fe is inwhere Ke is the stiffness matrix of elements (equation (4)) and D is the elastic matrix:
We combine each element stiffness matrix and element external force according to the global index of the mesh nodes to obtain the global stiffness matrix K and global external force f for the whole mesh. The relationship between global node displacement and global external force vector f can be obtained in
When it is necessary, solve equation (5) to calculate the deformation displacement u. K is determined in preprocessing which costs no time in deformation calculation.
4. Position-Based Dynamics
In PBD, it is necessary to predict the position and velocity of each node and update the nodal displacement to proper position by constraint function. The establishment of the constraint function affects the stability and efficiency of the model. Thus, solving the constraint function becomes the most important part of PBD. Let the constraint function be C(P) = 0. For the nodal predicted position P, the nodal position projected by the constraint functions is P + ∆P. So, constraint function C(P + ∆P) iswhere ∆P is the correction of P after the constraint projection and ∇PC(P) is the gradient of the constraint function. In order to maintain the conservation of momentum and angular momentum, ∆P and ∇PC(P) have the same direction. Thus, a Lagrangian multiplier λ can be used to express ∆P and ∇PC(P) relationship:
In our method, we utilize two constraints to perform the deformation of fine surface mesh: stretch constraint and bending constraint. Two constraints are as follows.
4.1. Stretch Constraint
As Figure 2 shows, the stretch constraint simulates the elastic force between two nodes.
The stretch constraint represents equation (8), where l is the initial length of P1P2.
According to equation (7), we obtained ∆P1 and ∆P2 in equation (9), where = 1/mi (i = 1, 2). m1 and m2 are the mass of P1 and P2:
4.2. Bending Constraint
As Figure 3 shows, the constraint simulates the bending between two triangle elements.
The bending constraint function iswhere n1 and n2 are the normal of ∆P1P3P2 and ∆P1P2P4 and φ is the initial angle between ∆P1P3P2 and ∆P1P2P4. According to equation (7), we obtained ∆Pi (i = 1, 2, 3, 4) inwhere = 1/mi, mi is the mass of Pi, and qi (i = 1, 2, 3, 4) is
5. Construction of the Hybrid Mesh
The hybrid mesh includes a fine triangle surface mesh and a matching internal coarse tetrahedral volume mesh. Now, given a fine triangle surface mesh, the steps to generate the hybrid mesh are as follows:(i)Step 1: input the fine triangle surface mesh MS(ii)Step 2: simplify MS to a coarse triangle surface mesh (iii)Step 3: according to Delaunay algorithm offered in open-source software Tetgen, use to generate a coarse tetrahedral volume mesh MV(iv)Step 4: combine MS and MV to obtain the hybrid mesh
The Illustration of these steps is shown in Figure 4. Since Step 3 uses the open-source algorithm, we do not introduce it here [21]. Now, Step 2 is introduced as follows. Given MS, it is necessary to select some nodes from MS as the nodes of the coarse surface triangle mesh . The selection of these nodes is based on the principle of minimum vertices cover; i.e., some nodes in MS are selected to a node set V, so that each edge in MS contains at least one node in V and the number of nodes in V is minimized. We use the method in [22] to obtain the V.
After V is obtained, the way to simplify MS to is introduced in the following. As shown in Figure 5, for one node pi in MS and not in V, a node pj which is in V and adjacent to pi is chosen. After that, pi is merged with pj and the edge pipj is removed. For all adjacent vertices pk of pi, if pk is also adjacent to pj, the edges pipk are also removed. Otherwise, a new edge pkpj will be created after pipk removed.
Then, the selection of pj will become a problem. In order to find a proper pj, we first calculate the average distance davg between pi and pk:where dik is the distance between pi and pk. And then, we select pj that satisfies
The node pj is selected in this condition in order to guarantee a better mesh quality after simplification. When all the nodes pi are simplified in this way, a coarse surface triangle mesh is obtained. Step 2 is now complete.
6. Deformation Algorithm Based on the Hybrid Mesh
The algorithm includes processing of external force, computing the deformation of the coarse volume mesh (MV) using FEM, deformation of MV, position prediction of each node on the fine surface mesh (MS), and constraint projection on MS. Figure 6 shows the illustration of the deformation algorithm on the 2D sectional view of the hybrid mesh. The black points indicate the shared nodes of MV and MS. These nodes are obtained from the MS simplification in Section 5. We call these nodes control nodes/points. Control nodes share the same position information with nodes on the boundary of MV. The white points represent normal nodes in MS. MV is represented by dotted lines, which means that it is not rendered in graphics. The solid line at the bottom of MV indicates that it is subject to displacement constraints.
6.1. Processing of External Force
Figure 6(a) shows that a vertical upward force is applied to the hybrid mesh. Since MS has many nodes, the default load point is on the node. If not, snap the load point to the closest node. After snapping, as Figure 6(b) shows, the position of load point may not be on the control node. In this condition, as Figure 6(c) shows, the equivalent force on the surrounding nodes of the external force is calculated for MV. The relationship between equivalent force and external force satisfies linear interpolation as the shape function of the triangular element.
6.2. Deformation of MV
After processing of external force, deformation of MV is solved according to equation (5). Control nodes share the same position with the nodes on the boundary of MV, so the position of control nodes is also obtained. Solving equation (5) will not take too long to run since MV is coarse.
6.3. Position Prediction of MS
For position prediction of MS, like PBD, the node position changes under the external load, without considering the influence of internal forces between the nodes in MS. This is node position prediction. The final node position of MS will be determined after constraint projection. Take the MS in Figure 6(d) as an example, where only one node is affected by the external load. Then, its position changes. The other vertices are not affected by external loads, so their positions have not changed. In this time, the nodes on MS are ready for constraint projection.
6.4. Constraint Projection of MS
As Figure 6(e) shows, since control nodes share the same position on MS and the boundary of MV, the displacement of control nodes in MS is determined when deformation of MV completes. In addition, constraint projection is also conducted for position correction, and the deformation position of normal nodes in MS is determined. Since only MS participates in graphic rendering and MV does not, this is the final state of hybrid mesh at the current moment.
6.5. Algorithm Flow
Step 1. Initialization: compute the global stiffness matrix K of MV, and initialize the position pi and velocity of each node of MS and the external force f. Step 2. In the current ∆t, get f, compute the deformation of MV according to equation (11), obtain the position of nodes on the boundary of MV, and update the velocity of the nodes on MS according to where = 1/mi (mi is the node mass of MS). Step 3. Predict the position of each node on MS according to Step 4. Project constraints on MS: this step has two cases. For control nodes, = . For other nodes on MS, solve position correction ∆pi on MS according to equations (9) and (11) and obtain the according to Step 5. According to equation (18), update the pi and of each node on MS: Step 6. Start the next ∆t and return to Step 2.
7. Small Incision Processing
In some deformation simulations, operations such as incision on mesh are necessary. Generally, the incision on the surface of a model is small, so the incision displacement of MV can be ignored. Moreover, MV is not rendered in graphics and thus we focus on the incision of MS. To this end, we take the incision treatment as a linear segment which is intersected by MS and a sweep surface generated by the virtual cutting tool. When the operation of incision happens, calculate the intersection points on MS. There are two types of the intersection points: edge points and facet points. When the operation ends, subdivide the triangular elements and update topology change. Besides, the update of constraint on MS is also needed.
As Figure 7 shows, there are three element types for mesh topology change. We use the method of [23] to subdivide the triangular mesh.
(a)
(b)
(c)
To update the constraint functions, the original constraints of the elements in the incision area are deleted. And, the new constraints of the newly subdivided triangular elements are generated. The new constraints generation is introduced in Figure 8.
(a)
(b)
For the stretch constraint, as Figure 8(a) shows, the intersection P is an edge point. Calculate the barycentric coordinate of P (λ1, λ2, λ3) (in Figure 8(a), λ3 = 0). The initial length of AP lAP iswhere lAB is initial length of AB in the original constraint. For stretch constraint, as Figure 8(b) shows, the intersection P is the facet point. The barycentric coordinate of P is (λ1, λ2, λ3). Take the calculation of initial length AP lAP as an example. According to law of cosine, angle α iswhere lAB, lBC, and lAC are the initial lengths in original constraints AB, BC, and AC, respectively. The initial length of DC lDC is
According to the law of cosine, the initial length of AD lAD is
The initial length of AP lAP is
For the bending constraint, the subdivision of triangle elements does not change the initial angle of two original triangles. Thus, the initial angle of triangular subelements remains unchanged. Since a control point is in the intersection segment, the control point P is duplicated due to the element subdivision. These two newly duplicated points are no longer control points.
8. Experiments and Results
8.1. Experimental Environment
We conduct our experiments on a computer using Intel Core i3 (2.7 GHz), 4 GB RAM, HD530 integrated graphics, and Windows 10. The software platform consists of VC++2019, open-source SOFA framework, and OpenGL graphics library.
8.2. Deformation Simulation
In order to present the mesh deformation effect, we implemented our proposed method on a simple deformation simulation.
In simulation, we use a pin-head to apply the external force f on the mesh. As Figure 9(a) shows, at the beginning, the mouse cursor controls the pin-head. When the cursor (pin-head) does not touch the mesh, the position of the pin-head is consistent with the position of the cursor. Once the cursor gets close to the node of mesh on the surface, as Figure 9(b) shows, the pin-head snaps to the closest node if the left mouse button is pressed down. In this situation, the position of pin-head is consistent with this node, not the cursor, and it becomes the external force-applied node. Then, as Figure 9(c) shows, move the cursor to generate external force f.
(a)
(b)
(c)
x is the current position of the cursor. x0 is the original cursor position when the pin-head snaps to the mesh node. The coefficient k can be set by user’s preference. The f is provided to the deformation algorithm. If the left mouse button is released, the external force is withdrawn and the position of pin-head is consistent with the cursor again. The cursor is not rendered in graphics. The flow chart is shown in Figure 10.
Figure 11 shows the deformation effect on two shape models (cube and sphere) and two real mesh models (liver and steak) with tension and pressure. The two real mesh models are from turbosquid.com. We set the material properties as follows: Young’s modulus: 5.5 × 103 Pa; Poisson’s ratio: 0.3. By observing the four models, the surface deformation effects are realistic and the graphic rendering is good.
8.3. Surface Deformation Performance
8.3.1. Accuracy Verification
In order to show the accuracy of our proposed method, we compare the deformation result of liver and steak mesh simulated by our method with FEM [10], PBD [14], and MSM [17]. The mesh used in FEM is a fine tetrahedral volume mesh, which is directly generated by MS through the Delaunay algorithm in Tetgen. The reason we use fine volume mesh in FEM is that the deformation experimental results will be used as the reference in experiment for its high accuracy. The mesh used in PBD and MSM is the same surface fine mesh model as the MS in our hybrid mesh. The information of cube, sphere, liver, and steak mesh models used in different methods is listed in Table 1.
The specific process is as follows: First, we directly applied vertical upward tension and vertical downward pressure on the liver and steak mesh. The magnitudes are 0.5 N, 1.0 N, and 2.0 N for small deformation and 5.0 N, 10.0 N, and 20.0 N for large deformation. Second, we choose 10 points near the force-applied points of the mesh as sample points. Finally, the positions of the sample points after deformation of two shape models (cube and sphere) and two real mesh models (liver and steak) are recorded, and the deformation displacement d are calculated by usingwhere represents the position of the sample points i after the deformation and represents the initial position of the sample points i. The deformation displacement radar diagrams of cube, sphere, liver, and steak models are shown in Figures 12–15 respectively.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
Meanwhile, we calculate the RMSE (root mean square error) of the deformation displacement of our method, PBD, and MSM. Note that the RMSEs are obtained based on the FEM as reference:where r represents the RMSE value, n represents the number of the deformation sample points, represents the position of the deformation sample points i after the deformation of our method, PBD, and MSM, respectively, and represents the position of the deformation sample points i after the deformation of FEM as reference. For RMSE calculation, we choose 50 points nearest to the load-applied point (include the load-applied point) as sample points. Table 2 shows RMSE of deformation displacement of our method, PBD, and MSM on liver and steak mesh models. Table 3 shows RMSE of deformation displacement of our method, PBD, and MSM on cube and sphere shapes.
From the diagrams shown in Figures 12 and 13, for both small deformation (0.5 N, 1.0 N, 2.0 N) and large deformation (5.0 N, 10.0 N, 20.0 N), the deformation displacement curve of our method is closer to that of FEM than PBD and MSM. According to the definition of RMSE, the closer the RMSE value is to zero, the closer the deformation displacements is to the value of FEM. The RMSE value of our method is closer to zero than the value of PBD and MSM on both small deformation and large deformation. The above results verify the accuracy of our method.
8.3.2. Time Verification
For time verification, we compared the consumed time of our method with that of PBD, MSM, and FEM in calculating the deformation results. As shown in Tables 4 and 5, the time consumed by our method, PBD, MSM, and FEM under different forces is recorded on two shape models (cube and sphere) and two real mesh models.
From the result shown in Tables 3 and 4, on the one hand, compared with PBD and MSM, our method costs more computation time due to the addition of the time consumed for MV deformation, but this time increment is not too much. On the other hand, compared with FEM, our method reduces large computation cost, since we use coarse volume mesh instead of fine volume mesh in FEM for deformation calculation. Although the study in [18] has reduced some computational cost of FEM by optimizing the numerical schemes of dealing with the stretch part of deformation energy, it needs more elements to present the accurate surface deformation effect than ours. Our proposed method needs fewer parameters of volume element and surface element to present the effect. The real-time performance of our proposed method is guaranteed.
8.3.3. Discussion of the Trade-Off between Accuracy and Time Consumption
In accuracy verification, by comparing the deformation displacement and RMSE of our method with the ones of PBD, MSM, and FEM, the results indicates that the deformation of our proposed method is similar to the one of FEM, which we take as reference to evaluate the accuracy. The deformation result of our method is more accurate, due to the constraint on MS from the displacement of MV. For the computation time, by comparing the consumed time of our methods with that of PBD, MSM, and FEM, the result indicates that our method has lower time consumption than FEM. Although our method has more time consumption than PBD and MSM, the time increment is not too much since computing the deformation for MV does not take too much time. Comprehensively considering the deformation result and time consumption, our proposed method achieves the trade-off between the high accuracy and low time consumption.
8.4. Small Incision Effect
We also developed a small incision simulation using our hybrid mesh model. We use a pin-head as a virtual cutting tool. If the pin-head is close to the mesh surface and the left mouse button is pressed down, the incision operation begins. As the pin-head moves, the intersection happens on the mesh surface. When the left mouse button is released, the incision operation is finished. Figure 16 shows the small incision simulation effects on the liver and steak mesh model. For the incision to be visible, we applied a tiny force on a node near the incision on the surface mesh. From Figure 16, we can see that the incisions of these two shape models and two real mesh models are smooth. The incision processing is compatible with our proposed method.
(a)
(b)
(c)
(d)
9. Conclusions
In this paper, we propose a deformation method based on the hybrid mesh model which achieves the trade-off between the high accuracy and low time consumption. First, we proposed the construction of the hybrid mesh from fine surface triangular mesh. Its most important step is the simplification process of the fine surface mesh. Second, we combined the FEM for coarse volume mesh deformation and PBD deformation method for fine surface mesh deformation and proposed the algorithm flow based on the hybrid mesh model. Third, we implemented our method on a simple deformation simulation. The experimental result of the deformation algorithm is presented. Compared with the PBD, MSM, and FEM, our method verifies the high accuracy and low time consumption. Finally, we developed an incision simulation which is compatible with our proposed method.
In the future, our proposed method can be improved further in discontinuity processing. The cutting on the hybrid mesh model is a major problem. It involves the subdivision of MV, the way to calculate the deformation of the subdivided MV, and how the surface node is constrained on the subdivided MV. When a large cut happens, the graphic rendering of incision on MV is a problem. It is a challenging task to maintain the trade-off between high accuracy and low time consumption in deformation simulation with subdivision of MS and MV.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the NSFC Joint Fund under Grants GG2090090072, U1332130, and U1713206; 111 Projects under Grant B07033; 973 Project under Grant 2014CB931804; and Key Research and Development Plan of Anhui Province under Grant 1704a0902051.