# A Hybrid Model Method for Accurate Surface Deformation and Incision Based on FEM and PBD

Authors: Shijie Tan; Hongjun Zhou; Jinjin Zheng
Journal: Scientific Programming
 (2021)

---

## 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.


---

## Body

## 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: Section2 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 Figure1). 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.Figure 1
Overview of our proposed method for accurate surface deformation.


## 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.(1)uex=Nexue.Nex is the shape function matrix of elements, and ue 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:(2)εex=Bexue.Bex is the strain matrix of elements. According to the principle of virtual work, the relationship between nodal displacement ue and external nodal force of elements fe is in(3)Keue=fe,where Ke is the stiffness matrix of elements (equation (4)) and D is the elastic matrix:(4)Ke=∫BeTDBedx.We combine each element stiffness matrix and element external force according to the global index of the mesh nodes to obtain the global stiffness matrixK and global external force f for the whole mesh. The relationship between global node displacement u and global external force vector f can be obtained in(5)Ku=f.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 beC(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) is(6)CP+ΔP≈CP+∇PCPΔP,where ∆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:(7)ΔP=λ∇PCP=−CP∇PCP2∇PCP.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 Figure2 shows, the stretch constraint simulates the elastic force between two nodes.Figure 2
Stretch constraint.The stretch constraint represents equation (8), where l is the initial length of P1P2.(8)CstretchP1,P2=P1−P2−l.According to equation (7), we obtained ∆P1 and ∆P2 in equation (9), where wi  = 1/mi (i = 1, 2). m1 and m2 are the mass of P1 and P2:(9)ΔP1=−w1w1+w2P1−P2−lP1−P2P1−P2,ΔP2=+w2w1+w2P1−P2−lP1−P2P1−P2.


### 4.2. Bending Constraint

As Figure3 shows, the constraint simulates the bending between two triangle elements.Figure 3
Bending constraint.The bending constraint function is(10)CbendingP1,P2,P3,P4=arccosn1⋅n2−φ,where 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) in(11)ΔPi=−4wi∑​jwj1−n1⋅n22arccosn1⋅n2−φ∑​jqj2qi,where wi = 1/mi, mi is the mass of Pi, and qi (i = 1, 2, 3, 4) is(12)q1=−q2−q3−q4,q2=−P3×n2+n1×P3n1⋅n2P2×P3−P4×n1+n2×P4n1⋅n2P2×P4,q3=P2×n2+n1×P2n1⋅n2P2×P3,q4=P2×n1+n2×P2n1⋅n2P2×P4.


## 4.1. Stretch Constraint

As Figure2 shows, the stretch constraint simulates the elastic force between two nodes.Figure 2
Stretch constraint.The stretch constraint represents equation (8), where l is the initial length of P1P2.(8)CstretchP1,P2=P1−P2−l.According to equation (7), we obtained ∆P1 and ∆P2 in equation (9), where wi  = 1/mi (i = 1, 2). m1 and m2 are the mass of P1 and P2:(9)ΔP1=−w1w1+w2P1−P2−lP1−P2P1−P2,ΔP2=+w2w1+w2P1−P2−lP1−P2P1−P2.


## 4.2. Bending Constraint

As Figure3 shows, the constraint simulates the bending between two triangle elements.Figure 3
Bending constraint.The bending constraint function is(10)CbendingP1,P2,P3,P4=arccosn1⋅n2−φ,where 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) in(11)ΔPi=−4wi∑​jwj1−n1⋅n22arccosn1⋅n2−φ∑​jqj2qi,where wi = 1/mi, mi is the mass of Pi, and qi (i = 1, 2, 3, 4) is(12)q1=−q2−q3−q4,q2=−P3×n2+n1×P3n1⋅n2P2×P3−P4×n1+n2×P4n1⋅n2P2×P4,q3=P2×n2+n1×P2n1⋅n2P2×P3,q4=P2×n1+n2×P2n1⋅n2P2×P4.


## 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 meshMS(ii)
Step 2: simplifyMS to a coarse triangle surface mesh MS′(iii)
Step 3: according to Delaunay algorithm offered in open-source software Tetgen, useMS′ to generate a coarse tetrahedral volume mesh MV(iv)
Step 4: combineMS and MV to obtain the hybrid meshThe Illustration of these steps is shown in Figure4. 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 MS′. 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.Figure 4
Construction of the hybrid mesh.AfterV is obtained, the way to simplify MS to MS′ 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.Figure 5
Simplification ofMS to MS′.Then, the selection ofpj will become a problem. In order to find a proper pj, we first calculate the average distance davg between pi and pk:(13)davg=1n∑kndik,where dik is the distance between pi and pk. And then, we select pj that satisfies(14)dij=minkdik−davg.The nodepj 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 MS′ 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.Figure 6
Overview of our proposed deformation algorithm. (a–c) Processing of external force. (d) Deformation ofMV and position prediction of MS. (e) constraint projection of MS.


### 6.1. Processing of External Force

Figure6(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 ofMV

After processing of external force, deformation ofMV 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 ofMS

For position prediction ofMS, 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 ofMS

As Figure6(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 vi 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 pi′ on the boundary of MV, and update the velocity vinew of the nodes on MS according to (15)vinew=vi+fΔtwi,wherewi  = 1/mi (mi is the node mass of MS).Step 3. Predict the position of each node pinew on MS according to(16)pinew=pi+vinewΔt.Step 4. Project constraints on MS: this step has two cases. For control nodes, pisolve = pi′. For other nodes on MS, solve position correction ∆pi on MS according to equations (9) and (11) and obtain the pisolve according to(17)pisolve=pi+Δpi.Step 5. According to equation (18), update the pi and vi of each node on MS:(18)vi=pisolve−piΔt,pi=pisolve.Step 6. Start the next ∆t and return to Step 2.


## 6.1. Processing of External Force

Figure6(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 ofMV

After processing of external force, deformation ofMV 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 ofMS

For position prediction ofMS, 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 ofMS

As Figure6(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 vi 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 pi′ on the boundary of MV, and update the velocity vinew of the nodes on MS according to (15)vinew=vi+fΔtwi,wherewi  = 1/mi (mi is the node mass of MS).Step 3. Predict the position of each node pinew on MS according to(16)pinew=pi+vinewΔt.Step 4. Project constraints on MS: this step has two cases. For control nodes, pisolve = pi′. For other nodes on MS, solve position correction ∆pi on MS according to equations (9) and (11) and obtain the pisolve according to(17)pisolve=pi+Δpi.Step 5. According to equation (18), update the pi and vi of each node on MS:(18)vi=pisolve−piΔt,pi=pisolve.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 ofMV 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 Figure7 shows, there are three element types for mesh topology change. We use the method of [23] to subdivide the triangular mesh.Figure 7
Three types of topology change: (a) uncut with an edge point; (b) partial cut; (c) full cut.
(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 Figure8.Figure 8
Generation of new stretch constraint: (a) intersectionP on edge; (b) intersection P on facet.
(a)(b)For the stretch constraint, as Figure8(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 is(19)lAP=λ1λ1+λ2lAB,where 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 α is(20)cosα=lBC2+lAC2−lAB22lBClAC,where lAB, lBC, and lAC are the initial lengths in original constraints AB, BC, and AC, respectively. The initial length of DC lDC is(21)lDC=λ2λ2+λ3lBC.According to the law of cosine, the initial length of ADlAD is(22)lAD=lDC2+lAC2−2lDClACcosα.The initial length of APlAP is(23)lAP=λ2+λ3lAD.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 pointP 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 forcef 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.(24)f=kx−x0.Figure 9
Diagram of deformation simulation: (a) pin-head untouched; (b) pin-head touched; (c) deformation.
(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 10
Flowchart of the simple deformation simulation.Figure11 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.Figure 11
The surface deformation effects of four models in two load cases. Models (from top to bottom): cube, sphere, liver, and steak. Load cases: left, tension; right, pressure.


### 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.Table 1
The information of four models.
ModelsFine surface mesh (MS)Coarse volume mesh (MV)Fine volume mesh#node#elemSize (kB)#node#elemSize (kB)#node#elemSize (kB)Cube12912888106436122045302814874465Sphere14023158115587163160344117821538Liver13793122114543159258335716985526Steak15433378126654171965386520898623#node: number of nodes; #elem: number of elements; size is expressed in kB. Our method:MS + MV; PBD and MSM: MS; FEM: fine volume mesh generated by MS.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 displacementd are calculated by using(25)d=xsi−x0i,where xsi represents the position of the sample points i after the deformation and x0i 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.Figure 12
The surface deformation displacement with respect to the sample nodes of the cube model. (a) Small deformation: from left to right are 0.5 N, 1.0 N, and 2.0 N. Top indicates tension and bottom indicates pressure. (b) Large deformation: from left to right are 5.0 N, 10.0 N, and 20.0 N. Top indicates tension and bottom indicates pressure.
(a)(b)Figure 13
The surface deformation displacement with respect to the sample nodes of the sphere model. (a) Small deformation: from left to right are 0.5 N, 1.0 N, and 2.0 N. Top indicates tension and bottom indicates pressure. (b) Large deformation: from left to right are 5.0 N, 10.0 N, and 20.0 N. Top indicates tension and bottom indicates pressure.
(a)(b)Figure 14
The surface deformation displacement with respect to the sample nodes of the liver model. (a) Small deformation: from left to right are 0.5 N, 1.0 N, and 2.0 N. Top indicates tension and bottom indicates pressure. (b) Large deformation: from left to right are in 5.0 N, 10.0 N, and 20.0 N. Top indicates tension and bottom indicates pressure.
(a)(b)Figure 15
The surface deformation displacement with respect to the sample nodes of the steak model. (a) Small deformation: from left to right are 0.5 N, 1.0 N, and 2.0 N. Top indicates tension and bottom indicates pressure. (b) Large deformation: from left to right are 5.0 N, 10.0 N, and 20.0 N. Top indicates tension and bottom indicates pressure.
(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:(26)r=∑i=1nxsi−xfi2n,where r represents the RMSE value, n represents the number of the deformation sample points, xsi represents the position of the deformation sample points i after the deformation of our method, PBD, and MSM, respectively, and xfi 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.Table 2
RMSE of deformation displacement of our method, PBD, and MSM on two real mesh models.
RMSE (mm)LiverSteakOursPBDMSMOursPBDMSMTension0.5 N0.08500.17840.20670.08780.18010.20801.0 N0.04760.13660.22710.04870.14500.23882.0 N0.09520.27320.45420.09740.28990.47765.0 N0.21140.85801.07080.25950.88811.024510.0 N0.52881.80452.27970.57711.79842.372320.0 N1.02204.07985.51881.03824.08985.4101Pressure0.5 N0.04760.13660.22710.08540.18450.21561.0 N0.05720.17400.29910.05110.16500.26602.0 N0.09520.27320.45420.09720.32990.49205.0 N0.28870.88761.06900.25750.87651.007710.0 N0.53471.79922.30370.58171.80912.251620.0 N1.03403.90925.26671.03423.90905.3660Table 3
RMSE of deformation displacement of our method, PBD, and MSM on two shape models.
RMSE (mm)CubeSphereOursPBDMSMOursPBDMSMTension0.5 N0.02950.10040.13770.03020.10210.14551.0 N0.04620.16130.20610.04930.14320.24742.0 N0.09990.29480.48190.09570.26180.5115.0 N0.22240.86140.98090.22660.84370.936410.0 N0.5691.80812.42560.53551.69412.182520.0 N0.95454.17775.88861.03724.31885.4804Pressure0.5 N0.03050.10150.13870.02810.11630.13271.0 N0.04720.17020.21190.04620.18110.22832.0 N0.09620.29670.49220.10090.30990.50875.0 N0.21610.88420.99740.22070.87751.027910.0 N0.57341.81642.52530.58171.82722.513720.0 N1.07024.02415.60921.11644.04825.7581From the diagrams shown in Figures12 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 Tables4 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.Table 4
Deformation calculation time consumption of our method, PBD, MSM, and FEM on two shape models.
Time (ms)CubeSphereOursPBDMSMFEMOursPBDMSMFEMTension0.5 N3.82.41.911.63.92.62.210.91.0 N3.92.22.012.53.82.12.211.82.0 N3.62.32.212.63.92.12.312.15.0 N3.72.31.913.63.82.32.413.510.0 N3.82.42.213.03.72.52.113.320.0 N4.02.62.213.73.92.32.113.9Pressure0.5 N3.62.32.111.13.72.12.411.11.0 N3.72.22.312.63.82.32.212.12.0 N3.82.32.113.33.92.32.412.75.0 N3.92.52.413.33.72.72.413.110.0 N3.92.52.314.43.82.52.213.320.0 N4.02.62.214.94.02.32.013.2Table 5
Deformation calculation time consumption of our method, PBD, MSM, and FEM on two real mesh models.
Time (ms)LiverSteakOursPBDMSMFEMOursPBDMSMFEMTension0.5 N3.82.22.18.23.62.42.316.81.0 N3.72.32.19.33.62.22.217.22.0 N3.92.22.210.33.82.32.317.95.0 N3.62.12.112.53.52.32.219.610.0 N3.82.32.112.13.72.42.120.720.0 N3.72.42.313.83.92.22.122.8Pressure0.5 N3.72.22.17.93.62.32.318.41.0 N3.62.32.29.13.72.42.218.82.0 N3.82.22.19.93.92.52.319.35.0 N3.92.42.211.43.72.62.320.110.0 N3.62.32.112.23.82.32.121.420.0 N3.82.42.113.13.72.52.223.4From the result shown in Tables3 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 onMS 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. Figure16 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.Figure 16
The small incision effects of two shape models and two real mesh models. (a) Cube, (b) sphere, (c) liver, and (d) steak.
(a)(b)(c)(d)


## 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 forcef 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.(24)f=kx−x0.Figure 9
Diagram of deformation simulation: (a) pin-head untouched; (b) pin-head touched; (c) deformation.
(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 10
Flowchart of the simple deformation simulation.Figure11 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.Figure 11
The surface deformation effects of four models in two load cases. Models (from top to bottom): cube, sphere, liver, and steak. Load cases: left, tension; right, pressure.


## 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.Table 1
The information of four models.
ModelsFine surface mesh (MS)Coarse volume mesh (MV)Fine volume mesh#node#elemSize (kB)#node#elemSize (kB)#node#elemSize (kB)Cube12912888106436122045302814874465Sphere14023158115587163160344117821538Liver13793122114543159258335716985526Steak15433378126654171965386520898623#node: number of nodes; #elem: number of elements; size is expressed in kB. Our method:MS + MV; PBD and MSM: MS; FEM: fine volume mesh generated by MS.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 displacementd are calculated by using(25)d=xsi−x0i,where xsi represents the position of the sample points i after the deformation and x0i 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.Figure 12
The surface deformation displacement with respect to the sample nodes of the cube model. (a) Small deformation: from left to right are 0.5 N, 1.0 N, and 2.0 N. Top indicates tension and bottom indicates pressure. (b) Large deformation: from left to right are 5.0 N, 10.0 N, and 20.0 N. Top indicates tension and bottom indicates pressure.
(a)(b)Figure 13
The surface deformation displacement with respect to the sample nodes of the sphere model. (a) Small deformation: from left to right are 0.5 N, 1.0 N, and 2.0 N. Top indicates tension and bottom indicates pressure. (b) Large deformation: from left to right are 5.0 N, 10.0 N, and 20.0 N. Top indicates tension and bottom indicates pressure.
(a)(b)Figure 14
The surface deformation displacement with respect to the sample nodes of the liver model. (a) Small deformation: from left to right are 0.5 N, 1.0 N, and 2.0 N. Top indicates tension and bottom indicates pressure. (b) Large deformation: from left to right are in 5.0 N, 10.0 N, and 20.0 N. Top indicates tension and bottom indicates pressure.
(a)(b)Figure 15
The surface deformation displacement with respect to the sample nodes of the steak model. (a) Small deformation: from left to right are 0.5 N, 1.0 N, and 2.0 N. Top indicates tension and bottom indicates pressure. (b) Large deformation: from left to right are 5.0 N, 10.0 N, and 20.0 N. Top indicates tension and bottom indicates pressure.
(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:(26)r=∑i=1nxsi−xfi2n,where r represents the RMSE value, n represents the number of the deformation sample points, xsi represents the position of the deformation sample points i after the deformation of our method, PBD, and MSM, respectively, and xfi 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.Table 2
RMSE of deformation displacement of our method, PBD, and MSM on two real mesh models.
RMSE (mm)LiverSteakOursPBDMSMOursPBDMSMTension0.5 N0.08500.17840.20670.08780.18010.20801.0 N0.04760.13660.22710.04870.14500.23882.0 N0.09520.27320.45420.09740.28990.47765.0 N0.21140.85801.07080.25950.88811.024510.0 N0.52881.80452.27970.57711.79842.372320.0 N1.02204.07985.51881.03824.08985.4101Pressure0.5 N0.04760.13660.22710.08540.18450.21561.0 N0.05720.17400.29910.05110.16500.26602.0 N0.09520.27320.45420.09720.32990.49205.0 N0.28870.88761.06900.25750.87651.007710.0 N0.53471.79922.30370.58171.80912.251620.0 N1.03403.90925.26671.03423.90905.3660Table 3
RMSE of deformation displacement of our method, PBD, and MSM on two shape models.
RMSE (mm)CubeSphereOursPBDMSMOursPBDMSMTension0.5 N0.02950.10040.13770.03020.10210.14551.0 N0.04620.16130.20610.04930.14320.24742.0 N0.09990.29480.48190.09570.26180.5115.0 N0.22240.86140.98090.22660.84370.936410.0 N0.5691.80812.42560.53551.69412.182520.0 N0.95454.17775.88861.03724.31885.4804Pressure0.5 N0.03050.10150.13870.02810.11630.13271.0 N0.04720.17020.21190.04620.18110.22832.0 N0.09620.29670.49220.10090.30990.50875.0 N0.21610.88420.99740.22070.87751.027910.0 N0.57341.81642.52530.58171.82722.513720.0 N1.07024.02415.60921.11644.04825.7581From the diagrams shown in Figures12 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 Tables4 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.Table 4
Deformation calculation time consumption of our method, PBD, MSM, and FEM on two shape models.
Time (ms)CubeSphereOursPBDMSMFEMOursPBDMSMFEMTension0.5 N3.82.41.911.63.92.62.210.91.0 N3.92.22.012.53.82.12.211.82.0 N3.62.32.212.63.92.12.312.15.0 N3.72.31.913.63.82.32.413.510.0 N3.82.42.213.03.72.52.113.320.0 N4.02.62.213.73.92.32.113.9Pressure0.5 N3.62.32.111.13.72.12.411.11.0 N3.72.22.312.63.82.32.212.12.0 N3.82.32.113.33.92.32.412.75.0 N3.92.52.413.33.72.72.413.110.0 N3.92.52.314.43.82.52.213.320.0 N4.02.62.214.94.02.32.013.2Table 5
Deformation calculation time consumption of our method, PBD, MSM, and FEM on two real mesh models.
Time (ms)LiverSteakOursPBDMSMFEMOursPBDMSMFEMTension0.5 N3.82.22.18.23.62.42.316.81.0 N3.72.32.19.33.62.22.217.22.0 N3.92.22.210.33.82.32.317.95.0 N3.62.12.112.53.52.32.219.610.0 N3.82.32.112.13.72.42.120.720.0 N3.72.42.313.83.92.22.122.8Pressure0.5 N3.72.22.17.93.62.32.318.41.0 N3.62.32.29.13.72.42.218.82.0 N3.82.22.19.93.92.52.319.35.0 N3.92.42.211.43.72.62.320.110.0 N3.62.32.112.23.82.32.121.420.0 N3.82.42.113.13.72.52.223.4From the result shown in Tables3 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 onMS 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.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.Table 1
The information of four models.
ModelsFine surface mesh (MS)Coarse volume mesh (MV)Fine volume mesh#node#elemSize (kB)#node#elemSize (kB)#node#elemSize (kB)Cube12912888106436122045302814874465Sphere14023158115587163160344117821538Liver13793122114543159258335716985526Steak15433378126654171965386520898623#node: number of nodes; #elem: number of elements; size is expressed in kB. Our method:MS + MV; PBD and MSM: MS; FEM: fine volume mesh generated by MS.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 displacementd are calculated by using(25)d=xsi−x0i,where xsi represents the position of the sample points i after the deformation and x0i 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.Figure 12
The surface deformation displacement with respect to the sample nodes of the cube model. (a) Small deformation: from left to right are 0.5 N, 1.0 N, and 2.0 N. Top indicates tension and bottom indicates pressure. (b) Large deformation: from left to right are 5.0 N, 10.0 N, and 20.0 N. Top indicates tension and bottom indicates pressure.
(a)(b)Figure 13
The surface deformation displacement with respect to the sample nodes of the sphere model. (a) Small deformation: from left to right are 0.5 N, 1.0 N, and 2.0 N. Top indicates tension and bottom indicates pressure. (b) Large deformation: from left to right are 5.0 N, 10.0 N, and 20.0 N. Top indicates tension and bottom indicates pressure.
(a)(b)Figure 14
The surface deformation displacement with respect to the sample nodes of the liver model. (a) Small deformation: from left to right are 0.5 N, 1.0 N, and 2.0 N. Top indicates tension and bottom indicates pressure. (b) Large deformation: from left to right are in 5.0 N, 10.0 N, and 20.0 N. Top indicates tension and bottom indicates pressure.
(a)(b)Figure 15
The surface deformation displacement with respect to the sample nodes of the steak model. (a) Small deformation: from left to right are 0.5 N, 1.0 N, and 2.0 N. Top indicates tension and bottom indicates pressure. (b) Large deformation: from left to right are 5.0 N, 10.0 N, and 20.0 N. Top indicates tension and bottom indicates pressure.
(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:(26)r=∑i=1nxsi−xfi2n,where r represents the RMSE value, n represents the number of the deformation sample points, xsi represents the position of the deformation sample points i after the deformation of our method, PBD, and MSM, respectively, and xfi 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.Table 2
RMSE of deformation displacement of our method, PBD, and MSM on two real mesh models.
RMSE (mm)LiverSteakOursPBDMSMOursPBDMSMTension0.5 N0.08500.17840.20670.08780.18010.20801.0 N0.04760.13660.22710.04870.14500.23882.0 N0.09520.27320.45420.09740.28990.47765.0 N0.21140.85801.07080.25950.88811.024510.0 N0.52881.80452.27970.57711.79842.372320.0 N1.02204.07985.51881.03824.08985.4101Pressure0.5 N0.04760.13660.22710.08540.18450.21561.0 N0.05720.17400.29910.05110.16500.26602.0 N0.09520.27320.45420.09720.32990.49205.0 N0.28870.88761.06900.25750.87651.007710.0 N0.53471.79922.30370.58171.80912.251620.0 N1.03403.90925.26671.03423.90905.3660Table 3
RMSE of deformation displacement of our method, PBD, and MSM on two shape models.
RMSE (mm)CubeSphereOursPBDMSMOursPBDMSMTension0.5 N0.02950.10040.13770.03020.10210.14551.0 N0.04620.16130.20610.04930.14320.24742.0 N0.09990.29480.48190.09570.26180.5115.0 N0.22240.86140.98090.22660.84370.936410.0 N0.5691.80812.42560.53551.69412.182520.0 N0.95454.17775.88861.03724.31885.4804Pressure0.5 N0.03050.10150.13870.02810.11630.13271.0 N0.04720.17020.21190.04620.18110.22832.0 N0.09620.29670.49220.10090.30990.50875.0 N0.21610.88420.99740.22070.87751.027910.0 N0.57341.81642.52530.58171.82722.513720.0 N1.07024.02415.60921.11644.04825.7581From the diagrams shown in Figures12 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 Tables4 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.Table 4
Deformation calculation time consumption of our method, PBD, MSM, and FEM on two shape models.
Time (ms)CubeSphereOursPBDMSMFEMOursPBDMSMFEMTension0.5 N3.82.41.911.63.92.62.210.91.0 N3.92.22.012.53.82.12.211.82.0 N3.62.32.212.63.92.12.312.15.0 N3.72.31.913.63.82.32.413.510.0 N3.82.42.213.03.72.52.113.320.0 N4.02.62.213.73.92.32.113.9Pressure0.5 N3.62.32.111.13.72.12.411.11.0 N3.72.22.312.63.82.32.212.12.0 N3.82.32.113.33.92.32.412.75.0 N3.92.52.413.33.72.72.413.110.0 N3.92.52.314.43.82.52.213.320.0 N4.02.62.214.94.02.32.013.2Table 5
Deformation calculation time consumption of our method, PBD, MSM, and FEM on two real mesh models.
Time (ms)LiverSteakOursPBDMSMFEMOursPBDMSMFEMTension0.5 N3.82.22.18.23.62.42.316.81.0 N3.72.32.19.33.62.22.217.22.0 N3.92.22.210.33.82.32.317.95.0 N3.62.12.112.53.52.32.219.610.0 N3.82.32.112.13.72.42.120.720.0 N3.72.42.313.83.92.22.122.8Pressure0.5 N3.72.22.17.93.62.32.318.41.0 N3.62.32.29.13.72.42.218.82.0 N3.82.22.19.93.92.52.319.35.0 N3.92.42.211.43.72.62.320.110.0 N3.62.32.112.23.82.32.121.420.0 N3.82.42.113.13.72.52.223.4From the result shown in Tables3 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 onMS 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. Figure16 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.Figure 16
The small incision effects of two shape models and two real mesh models. (a) Cube, (b) sphere, (c) liver, and (d) steak.
(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 ofMV, 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.



---

*Source: 8343312-2021-12-15.xml*