An Optimized Mass-spring Model with Shape Restoration Ability Based on Volume Conservation

Abstract

To improve the accuracy and realism of the virtual surgical simulation system, this paper proposes an optimized mass-spring model with shape restoration ability based on volume conservation to simulate soft tissue deformation. The proposed method constructs a soft tissue surface model that adopts a new flexion spring for resisting bending and incorporates it into the mass-spring model (MSM) to restore the original shape. Then, we employ the particle swarm optimization algorithm to achieve the optimal solution of the model parameters. Besides, the volume conservation constraint is applied to the position-based dynamics (PBD) approach to maintain the volume of the deformable object for constructing the soft tissue volumetric model base on tetrahedrons. Finally, we built a simulation system on the PHANTOM OMNI force tactile interaction device to realize the deformation simulation of the virtual liver. Experimental results show that the proposed model has a good shape restoration ability and incompressibility, which can enhance the deformation accuracy and interactive realism.

1. Introduction

Traditional surgical training system regards human corpses, rubber mannequins or animals as objects for training and simulation. However, human corpses cannot be reused after dissection, and the animal body is distinct from the human body in terms of physiological structure. In addition, corpses, especially human corpses, are in an extreme shortage [1]. These make traditional surgical training suffer from severe challenges. With the development of simulation technology, virtual surgical training system has become an alternative to the traditional surgical training system, which plays an important part in clinical medicine [2]. The virtual surgical training system not only renders a virtual environment with a sense of interaction and realism for the clinical practice but also helps doctors to practice repeatedly by making a reasonable surgical scheme, which improve straining efficiency and reduces the cost and risk of surgery. Therefore, virtual surgical training system is of great significance for the development of medical science [3]. 

Modeling soft tissue deformation is one of the most critical techniques in the virtual surgical training system [4]. The common models for soft tissue deformation are the finite element model (FEM) [5-7], mass-spring model (MSM) [8-12] and position-based dynamics(PBD) approach [13-15]. Since the FEM can capture the biological characteristics of soft tissue precisely by Young's modulus and Poisson's ratio, it has a good deformation accuracy. However, the computational cost is expensive when reorganizing the mesh topologies, therefore not suitable for real-time deformation. MSM is simple and easy to implement and has a small computation amount of deformation and reliable real-time performance. However, the parameters of the mass-spring model (MSM) not connected with the biological characteristics of the soft tissue, resulting in poor model stability and low deformation accuracy. The PBD approach adopts the structure of point clouds, which do not require a complex topological relationship between points, and can directly control the position of objects. In addition, the model can describe the volumetric characteristic of the object, has high speed, strong stability, and good controllability. However, its accuracy is lower than that of FEM and MSM.

To address the above-mentioned challenges, this paper proposed an optimized mass-spring model with shape restoration ability based on volume conservation to simulate the deformation of a virtual liver. Our proposed MSM incorporates a new flexion spring to construct the surface model, uses the particle swarm optimization algorithm to optimize the model parameters, employs the PBD approach to conduct volume conservation constraint for characterizing the volumetric characteristic. The contributions of this paper are included:

(i) Incorporate a new flexion spring to the MSM for resisting bending based on the angle between the initial position vector and final position vector after the external force is removed, so as to improve the shape restoration ability for the surface of the MSM, especially when large deformation occurs.

(ii) Utilize the particle swarm optimization algorithm to optimize the elastic coefficient and damping coefficient of the MSM. Choose the FEM [16] as the reference model, and establishes a fitness function between these two parameters and the Young’s modulus and Poisson’s ratio. Tune the parameters by minimizing the fitness function until the deformation of the proposed MSM is close enough to the FEM, thereby improving the accuracy and realism of the deformation.

(iii) Adopt the PBD approach to conduct the volume conservation constraint based on tetrahedron elements for the interior structure of soft tissue in order to ensure the incompressibility, i.e., maintain the volume of the soft tissue.

The rest of this paper is organized as follows: Section 2 illustrate our proposed model. Experimental results and analysis are demonstrated in Section 3. Finally, Section 4 gives a conclusion on our model. 

2. Method

This paper proposed an optimized MSM with shape restoration ability based on volume conservation, which can simulates the volumetric characteristic of soft tissue. A particle swarm optimization algorithm is applied to optimize the model parameters, which can characterize the biological characteristics of soft tissue. The proposed model effectively enhances the shape restoration ability, accuracy and realism of deformation, and simulates the incompressibility of the soft tissue. The modeling procedure is shown in Fig. 1.

Fig. 1. Modeling procedure

2.1 MSM with a Flexion Spring

The traditional MSM utilizes discrete mass points and springs together with dampers as the connection to create a surface model. When an external force is applied on the surface, the elastic force and the damping force are generated by the spring unit to resist the external force. Moreover, the external force will be transmitted to the adjacent point through the spring unit and then to the farther adjacent point, bringing about the soft tissue deformation. However, the traditional MSM has a weak shape restoration ability and cannot accurately simulate the deformation when the large deformation occurred. To increase the shape restoration ability, we incorporate a new flexion spring into the MSM to resist bending.

The common surface model in MSM has three different kinds of topological structures, namely triangular mesh, quadrilateral mesh, and hexagonal mesh, as shown in Fig. 2. The quadrilateral mesh cannot recover its original shape after deformation, the hexagonal mesh is difficult and complicated for simulation, and the triangular mesh has strong stability to model objects. Therefore, the triangular mesh is employed to construct the surface model, where each triangular mesh consists of three mass points and three spring units, and each spring unit is composed of a structure spring, a damper, and a flexion spring, as shown in Fig. 3. The structure spring is used to describe the linear elasticity of the soft tissue, the damper is used to represent the damping characteristic of the soft tissue, and the flexion spring is used to resist bending to improve the shape restoration ability of the soft tissue.

Fig. 2. The topological structure of surface model in MSM (A) Triangular mesh (B) Quadrilateral mesh (C) Hexagonal mesh 

Fig. 3. The structure of spring unit

In this paper, the surface model of the virtual liver is constructed on triangular meshes. Suppose that triangular meshes consist of a set of points N and a set of spring units M ,two arbitrary adjacent points x_i and x_j , i,j = 1,2,...,N and i ≠ j , are connected by the spring unit. Structure spring and flexion spring are following the Hooke's law. When the external force f_ext is applied on the point x_i , according to the Hooke's law, the structure spring force \(f_{i j}^{s}\) at time t is computed as follows:

\(f_{i j}^{s}=k_{i j}^{s}\left(\left\|l_{i j}\right\|-\left\|l_{i j}\right\|^{0}\right) \frac{l_{i j}}{\left\|l_{i j}\right\|}\)       (1)

l_ij = x_j(t) - x_i(t)       (2)

where k^s_ij is the elastic coefficient of the structure spring connecting x_i and x _j , l_ij is the current displacement vector between x_i and x _j at time t, || l_ij || represents the modulus of the vector l_ij, and || l_ij ||⁰ denotes the initial distance between x_i and x _j . Similarly, the damping force \(f_{i j}^{d}\) at time t is given as follows:

\(f_{i j}^{d}=c_{i j}^{d} \Delta l_{i j}\)       (3)

Δl_ij =|| l_ij ||− || l_ij ||⁰       (4)

where \(c_{i j}^{d}\) is the damping coefficient of the damper connecting x_i and x _j , and ∆l_y is the derivative of || l_ij ||− || l_ij ||⁰. 

Since the traditional MSM is difficult to achieve a good shape restoration ability during the deformation simulation, the soft tissue cannot recover its original shape when the external force exerted to model surface is removed. Good shape restoration ability is of great significance for simulation [17]. To obtain a high-realism deformation model, a new type of flexion spring is incorporated into the traditional MSM for resisting bending and flexion spring force, which ensures the connecting points can return their initial positions. The flexion spring force \(f_{i j}^{f}\) at time t is defined as follows:

\(f_{i j}^{f}=k_{i j}^{f} \cdot \frac{\theta_{i j}}{\left|x_{j}(0)-x_{i}(0)\right|} \cdot \frac{d_{i j}}{\left|d_{i j}\right|}\)       (5)

where \(k_{i j}^{f}\) is the elastic coefficient of the flexion spring connecting x_i and x_j , θ_ij is the angle between the vector x_j(0) - x_i(0) − and x_j(t) - x_i(t), | x_j(0) - x_i(0)| is the modulus of the vector x_j(0) - x_i(0), and d_ij is the direction vector of the flexion spring force, which can be defined as follows:

\(d_{i j}=\frac{\left[x_{j}(t)-x_{i}(t)\right] \otimes\left[x_{j}(0)-x_{i}(0)\right] \otimes\left[x_{j}(t)-x_{i}(t)\right]}{\left|\left[x_{j}(t)-x_{i}(t)\right] \otimes\left[x_{j}(0)-x_{i}(0)\right] \otimes\left[x_{j}(t)-x_{i}(t)\right]\right|}\)       (6)

where ⊗ represents the cross product operation. The angle θ_ij is computed as

\(\theta_{i j}=\arccos \left\{\frac{\left[x_{j}(t)-x_{i}(t)\right] \mathrm{o}\left[x_{j}(0)-x_{i}(0)\right]}{\left|x_{j}(t)-x_{i}(t)\right| \cdot\left|x_{j}(0)-x_{i}(0)\right|}\right\}\)       (7)

where o represents the dot product operation.

The function of the flexion spring is described in Fig. 4. We select a triangular mesh containing points x_i , x_j and x_k from the surface model and fix the position of points x_j and x_k . We can observe that, assuming an external force f_ext   is acted on the point x_i at time t , the point can move from the initial position x_i(0) to the position x_i (0) and produce displacement. When the external force is withdrawn, the point will eventually move to the position x_i(t) instead of the position x_i(0)  because its kinetic energy has not been exhausted, which means that the point x_i is unable to return to its initial position after the external force is removed. At this time, it is necessary to incorporate the flexion spring and apply it on the point x_i to guarantee it can return to the initial position x_i(0) from the current position x_i(t) , i.e., let the angle θ be equal to zero.

Fig. 4. Schematic diagram of the flexion spring

The points x_j and x_k are not fixed at time t in the real situation. The point x_j moves to the position x_j(t) from the initial position x_j(0), and the point x_k moves to the position x_k(t) from the initial position x_k(0) under the external force. So the flexion spring ought to be employed to resist bending, ensuring that the point xi can return to the initial position, which enables the angle θij between the initial vector x_i(0)-x_j(0) and the final vector x_i(t)-x_j(t) be equal to zero. Similarly, the angle θ_ik between vectors x_i(0)-x_k(0) and x_i(t)-x_k(t) also can be easily return to zero. Therefore, the soft tissue can successfully return to the original shape when deformation generates.

2.2 Parameter Optimization

A common problem in MSM is that the model parameters cannot characterize the biological characteristics of the soft tissue, resulting in poor model stability and low simulation accuracy. This is because the biological characteristics are expressed with Young's modulus and Poisson's ratio, while the elastic coefficient and damping coefficient in MSM are not associated with them. Therefore, for achieving a behavioral realism experienced by surgeons, these two parameters are often obtained by trial and error. In this paper, we take the FEM, which has high accuracy and realism as the reference model, and employ the particle swarm optimization algorithm to optimize the elastic coefficient and damping coefficient.

The particle swarm optimization algorithm is suitable for the problem of finding the optimal solution of the elastic coefficient and damping coefficient. The flow chart is presented in Fig. 5. In this paper, the objective of the optimization problem is to minimize the difference in position of the corresponding points in our MSM and the reference model FEM, i.e., find the optimal values of these parameters so that the deformation caused by our MSM is close enough to the one caused by the reference model FEM.

Fig. 5. The procedure of parameter optimization using particle swarm optimization algorithm

We discretize the model surface based on the FEM into a set of points N, and these points have the same position as points in our MSM. Assuming all structure springs have the same elastic coefficients, all flexion springs have the same elastic coefficients, all dampers have the same damping coefficients, and two models are subjected to the identical external force. The fitness function is defined to describe the optimization problem in our paper is achieve the optimal solution of the elastic coefficient and damping coefficient as follows:

\(\Theta\left(k^{s}, k^{f}, d\right)=\sum_{i}\left|x_{i}^{M S M}\left(k^{s}, k^{f}, d, f_{e x t}\right)\right|-\left|x_{i}^{F E M}\left(E, v, f_{e x t}\right)\right|\)       (8)

where Θ is the sum of Euclidean distances between all corresponding points in the two models, \(x_{i}^{M S M}\) and \(x_{i}^{F E M}\) are the position of the i^th point in MSM and FEM, respectively. k^s and k^f are the elastic coefficients of the structure spring and the flexion spring, respectively. d is the damping coefficient, E is the Young’s modulus, v is the Poisson’s ratio, and f_ext is the external force.

According to the particle swarm optimization algorithm [18], the update equation of the velocity and position of the particles is presented as follows:

v_i = w·v_i + c₁·rand() ·(pbest_i-x_i) + c₂ · rand() · (gbest - x_i)       (9)

x_i = x_i + v_i       (10)

where v_i  and x_i are the velocity and position of the i^th particle, respectively, w is the inertia weight coefficient, c₁ and c₂ are learning factors, rand() is a positive random number within [0,1] , pbest_i denotes the best position of the i^th particle from the beginning to current iteration step, and pbest shows the global optimal position of the particle group. 

Finally, optimal solution is achieved by minimizing the fitness function according to Eqs. (9) and (10) such that all  \(x_{i}^{M S M}\) in our MSM extremely approach the corresponding \(x_{i}^{F E M}\) in the reference model FEM with an appropriate set of parameters k^s, k^f and d .

2.3 Deformation Calculation

According to Newton’s second law, the dynamics equation of the MSM can be described as:

m_ia_i = F_i       (11)

\(F_{i}=f_{\text {int }}^{i}+f_{\text {ext }}^{i}\)       (12) 

where m_i , a_i , and F_i are the mass, acceleration and resultant force of point x_i ,respectively. \(f_{\text {int }}^{i}\)  represents an internal force acting on point x_i , and  \(f_{\text {ext }}^{i}\) denotes an external force applied on point x_i by user. 

Since point x_i is connected with n different points through spring units, the structure spring force \(f_{i}^{s}\), the damping force\(f_{i}^{d}\) and the flexion spring force \(f_{i}^{f}\) acting on point x_i are computed as follows, respectively:

\(f_{i}^{s}=\sum_{j=1}^{n} f_{i j}^{s}=\sum_{j=1}^{n} k_{i j}^{s}\left(\left\|l_{i j}\right\|-\left\|l_{i j}\right\|^{0}\right) \frac{l_{i j}}{\left\|l_{i j}\right\|}\)       (13)

\(f_{i}^{d}=\sum_{j=1}^{n} f_{i j}^{d}=\sum_{j=1}^{n} c_{i j}^{d} \Delta l_{i j}\)       (14)

\(f_{i}^{f}=\sum_{j=1}^{n} f_{i j}^{f}=\sum_{j=1}^{n} k_{i j}^{f} \cdot \frac{\theta_{i j}}{\left|x_{j}(0)-x_{i}(0)\right|} \cdot \frac{d_{i j}}{\left|d_{i j}\right|}\)       (15)

where n is the number of points linked with point x_i . Based on the Eqs. (13), (14) and (15), the internal force \(f_{\text {int }}^{i}\)  is computed as follows:

\(f_{\text {int }}^{i}=f_{i}^{s}+f_{i}^{d}+f_{i}^{f}\)       (16)

Therefore, Eq. (11) is transformed as Eq.(17)

\(m_{i} \frac{d^{2} x_{i}}{d t^{2}}+f_{i}^{d}+f_{i}^{s}+f_{i}^{f}=f_{e x t}^{i}\)       (17)

Since the displacement of point x_i is unknown, we introduce the velocity variable v_i and transform the second-order ordinary differential equation of Eq. (17) into a first-order ordinary differential equation, which is given as Eq. (18).

\(\left\{\begin{array}{l} v=\frac{d x}{d t} \\ \frac{d v}{d t}=\frac{d^{2} x}{d t^{2}} \\ \frac{d v}{d t}=\frac{f_{e x t}^{i}-f_{i}^{d}-f_{i}^{s}-f_{i}^{f}}{m} \end{array}\right.\)       (18)

A popular method to numerically solve the second-order ordinary differential equation is the Euler method [3]. We adopt the implicit Euler method to solve Eq. (17), as is shown in Eq. (19).

v^t = dt·f(v,t) + v^t-dt       (19) 

where v^t and v^t-dt are current and previous velocity, respectively. dt denotes the time step, and f(v,t) denotes the derivative of velocity.

In summary, the steps of the deformation calculation for our MSM are shown as follows. First, initialize the mass m_i , position x_i(t₀) and velocity v_i(t₀) of points x_i ,i =1,2,...,N . Then, obtain the internal force \(f_{\text {int }}^{i}(t)\)and external force \(f_{\text {ext }}^{i}(t)\) for each point at each time step using Eqs. (13), (14), (15) and (16). Finally, the corresponding velocity v_i and position x_i can be iteratively updated in the next time step according to Eqs. (18) and (19).

2.4 Volume Conservation Constraint

The optimized MSM incorporating a flexion spring is only used to simulate the deformation on the surface; it cannot characterize the volumetric characteristic of the soft tissue. Therefore, we employ the PBD approach to implement volume conservation constraint based on tetrahedron elements inside the soft tissue so as to simulate the volumetric characteristic and incompressibility, i.e., maintain the volume of the deformable object.

The PBD approach can be depicted in the following steps [19]:

(i) Initialize the position p_i and velocity v_i of each vertex in the tetrahedron.

(ii) Update the velocity and position of each vertex according to the external force and internal damping of the system at each time step ∆t , which is shown in Eqs. (20) and (21).

\(v_{i}^{n e w}=v_{i}+\Delta t \cdot f_{\text {ext }} \cdot w_{i}+\operatorname{Damp}\left(v_{i}\right)\)       (20)

\(p_{i}^{\text {new }}=p_{i}+\Delta t \cdot v_{i}^{\text {new }}\)       (21)

where \(w_{i}=\frac{1}{m_{i}}\)

(iii) Applying volume conservation constraint to the estimated position \(p_{i}^{\text {new }}\) , the updated position \(p_{i}^{s o l}\) can be obtained by Gauss-Seidel iteration. 

\(p_{i}^{s o l}=p_{i}^{n e w}+\Delta p\)       (22)

(iv) The final position \(p_{i}^{f i n a l}\) and velocity \(v_{i}^{\text {final }}\)of each vertex can be obtained. 

\(p_{i}^{\text {final }}=p_{i}^{\text {sol }}\)      (23)

\(v_{i}^{\text {final }}=\frac{p_{i}^{\text {final }}-p_{i}}{\Delta t}\)       (24)

The interior of the soft tissue is represent by a set of P vertices and Q constraints, and the tetrahedron meshes are composed of these vertices. Let p be the vertices vector [\(p_{1}^{T}, p_{2}^{T}, p_{3}^{T}, p_{4}^{T}\)]^T of a tetrahedron. For step (iii), given p we want to find a correction factor ∆p such that C(p +∆p) = 0, i.e., move p to a valid position which can satisfy the defined volume conservation constraint. 

For the correction factor of an individual point p_{i ,}i = 1,2,3,4 , we have

\(\Delta p_{i}=-\frac{C(p)}{\sum_{j} w_{j}\left|\nabla_{p_{j}} C(p)\right|^{2}} w_{i} \nabla_{p_{i}} C(p)\)       (25) 

where C(p) and ∇_pC(p) denote the constraint function and its gradient, respectively. 

The volume conservation constraint of a single tetrahedron element can be described as:

\(C\left(p_{1}, p_{2}, p_{3}, p_{4}\right)=\frac{1}{6}\left(p_{2}-p_{1}\right) \times\left(p_{3}-p_{1}\right) \cdot\left(p_{4}-p_{1}\right)-V_{0}\)       (26) 

where V₀ is the initial volume of a single tetrahedron element. 

As shown in Fig. 6, the tetrahedron element inside the model will be deformed accordingly when the deformation is generated on the surface, resulting in the change of volume of the soft tissue. However, after we define the volume conservation constraint, the four vertices of the tetrahedron will perform corresponding adjustments to ensure that its volume remains unchanged, so that we can maintain the volume of the soft tissue.

Fig. 6. The structure of a tetrahedron element 

In order to obtain the correction factor ∆p_i, the unknown parameter ∇_pi C(p) in Eq. (25) need be calculated. Therefore, the corresponding gradients of Eq. (26) are deduced as follows:

\(\nabla_{p_{1}} C\left(p_{1}, p_{2}, p_{3}, p_{4}\right)=\frac{1}{6}\left(p_{4}-p_{2}\right) \times\left(p_{3}-p_{2}\right)\)       (27)

\(\nabla_{p_{2}} C\left(p_{1}, p_{2}, p_{3}, p_{4}\right)=\frac{1}{6}\left(p_{3}-p_{1}\right) \times\left(p_{4}-p_{1}\right)\)       (28)

\(\nabla_{p_{3}} C\left(p_{1}, p_{2}, p_{3}, p_{4}\right)=\frac{1}{6}\left(p_{4}-p_{1}\right) \times\left(p_{2}-p_{1}\right)\)       (29)

\(\nabla_{p_{4}} C\left(p_{1}, p_{2}, p_{3}, p_{4}\right)=\frac{1}{6}\left(p_{2}-p_{1}\right) \times\left(p_{3}-p_{1}\right)\)       (30) 

Now all the necessary parameters in Eq. (25) are obtained, we can calculate the correction factor by substituting Eqs. (27), (28), (29) and (30) into Eq. (25), thus the final position of the vertices satisfying the volume conservation constraint can be achieved.

3. Experiment

3.1 Experiment Environment

We adopt the PHANTOM OMNI hand controller produced by Senable Technologies Company for force tactile interaction, and conduct experiments on the Windows 2010 operating system with the VC++ 2017 and OpenGL. We realize the deformation simulation of the virtual liver on the following hardware platform: Intel Core i9-7920X CPU, 2.90 GHz,32GB RAM and NVIDA GeForce RTX 2080Ti which is shown in Fig. 7. 

Fig. 7. Simulation environment 

To simulate our proposed soft tissue deformation model, we firstly obtain the medical liver image by CT scanning, then employ OpenGL to reconstruct the 3D geometric model of the image, and use illumination, texture mapping, etc. to render virtual surgery scene. Then the operators interact with the soft tissue by the PHANTOM OMNI device, and the force feedback information will feedback to them at the same time.

3.2 Simulation Results

To verify the validity of the proposed model and the realism of the soft tissue deformation,we regard the area around the action point of the external force as the deformation area, in which 40 points are selected and marked. The parameters in this experiment are set as follows: the number of surface points N₁ =40, the mass of each point m_i =0.66, the elastic coefficient of each structure spring k^s =0.035, the elastic coefficient of each flexion spring k^f =0.32, and the damping coefficient d =0.004. We simulate the deformation of the virtual liver by using the virtual probe to apply the stress of 0.8N, 1.6N and 2.2N to the action point,as shown in Fig. 8. We also can observe that the deformation effect is continuous, realistic and meet the authenticity of the human-computer interaction.

Fig. 8. The deformation of virtual liver under different stress 

3.3 Experimental Analysis

3.3.1 Analysis of Shape Restoration

In this paper, the deformation average error \(\bar{x}\) , deformation standard deviation _s and shape restoration time t_e are regarded as evaluation criterions to evaluate the shape restoration ability of the traditional MSM [8], volume MSM [3] and our MSM under the same external conditions.

(1) Deformation Average Error 

The deformation average error x is used to evaluate the accuracy of the shape restoration ability after the external force is withdrawn, as defined:

\(\bar{x}=\frac{\sum_{i=1}^{N_{1}}\left|x_{i}^{0}-x_{i}^{\infty}\right|}{N_{1}}\)       (31) 

where \(x_{i}^{0}\) is the original position of the i^th point, \(x_{i}^{\infty}\) is the position of the i^th point when the soft tissue deformation ends. 

(2) Deformation Standard Deviation

The deformation standard deviation s is employed to evaluate the stability of the shape restoration ability, as defined:

\(s=\sqrt{\frac{\sum_{i=1}^{N_{1}}\left(\left|x_{i}^{0}-x_{i}^{\infty}\right|-\bar{x}\right)^{2}}{N_{1}-1}}\)       (32)

(3) Shape Restoration Time

The shape restoration time t_e is applied to describe the time when the virtual probe removes from the model surface until the shape restoration stops. 

We select and mark points on the traditional MSM and volume MSM surface as the action point of the external force and deformation area, respectively, and these points correspond with the ones on our MSM. Two different operation modes, i.e., pull mode and push mode, are adopted to the virtual soft tissue surface for simulation. The specific experimental process of evaluating the shape restoration ability for the traditional MSM, volume MSM and our MSM is shown as follows: Firstly, we use the virtual probe to the action point by same external force and record the original position of all marked points. Secondly, we remove the virtual probe and record the final position of all marked points when the deformation of virtual soft tissue ends. Finally, we calculate the deformation average error \(\bar{x}\) ,the deformation standard deviation _s and the shape restoration time t_e for the three different models. The corresponding results are given in Table 1. 

Table 1. Comparisons of shape restoration ability

Table 1 indicates that the deformation average error \(\bar{x}\) , deformation standard deviation _s, and shape restoration time t_e of our proposed model are smaller than these of the traditional MSM and volume MSM in the pull mode and push mode, which means that our proposed MSM outperforms the traditional MSM and volume MSM.

3.3.2 Accuracy of the Model

We choose the FEM [16] as the reference model, and then evaluate the accuracy and realism of our proposed model by comparing the deformation position between the proposed model and the FEM model.

We select and mark the points on the FEM surface corresponding to the proposed model as the action point of the external force and deformation area, and use the virtual probe to the action point by the same external force. The speed of the probe is 1 cm/s. Besides, we record the deformation position of all marked points at 0.2s, 0.4s, 0.6s, 0.8s, 1.0s, 1.2s, 1.4s, 1.6s,1.8s, 2.0s, and calculate the average positional error, the maximal positional error and the minimal positional error of the corresponding points, as shown in Fig. 9.

Fig. 9. Position errors of corresponding points under the same external force

Fig. 9 shows that position errors between the corresponding points of the FEM and the proposed model under the same external force. Fig. 9 indicates the average error in corresponding positions between the FEM and the proposed model is less than 1mm, which is imperceptible by human eyes. Therefore, the deformation based on the proposed model approximates to the reference model, which illustrates that our proposed model can effectively improve the accuracy and realism of the soft tissue deformation. 

3.3 Incompressibility of the Model

In this paper, we compare the MSM with and without volume conservation constraint to evaluate the incompressibility of the soft tissue. Fig. 10 exhibits the curves of volume change based on these two different methods. The comparison shows that our proposed method has good volume conservation and more stable and closer to the initial volume value, which means our proposed model has good incompressibility.

Fig. 10. The curves of volume change based on these two different methods

3.3.4 Comprehensive Evaluation of the Simulation System

To verify the interaction performance of the simulation system based on our proposed model,four evaluation indicators are selected, which are interactive naturalness, visual fluency,presence authenticity, and system stability. We make a comparison between our model and the traditional MSM [8], the PBD method [14] and the FEM [16], and grade the four models based on the above four indicators.

Since the above four evaluation indicators cannot quantitatively provide the actual standard to measure the interaction performance of simulation system, and there is a fuzzy relationship between these indicators, we employ the fuzzy comprehensive evaluation method to grade the interaction performance of four models, which is described as the following three steps:

(1) Determine fuzzy relationship matrix

Determine the evaluation indicator set U = {U₁,U₂,U₃,U₄} and the evaluation grade set V = {V₁,V₂,V₃,V₄,V₅} of the evaluation method, in which U₁ , U₂ , U₃ and U₄ represent the interactive naturalness, visual fluency, presence authenticity and system stability,respectively; V₁ , V₂ , V₃ , V₄ and V₅ indicate five grades of interaction performance of simulation system, i.e., excellent, good, general, bad and worse, respectively. 

We randomly select 27 doctors from First Affiliated Hospital of Nanjing Medical University, including 9 interns, 8 residents, 6 associate chief physicians, and 4 chief physicians. The doctors interact with the virtual liver through PHANTOM OMNI device,based on our proposed model, the traditional MSM, the PBD method and the FEM for 10 minutes respectively, and then evaluate the four indicators after the interaction. Construct the membership of the evaluation indicator to evaluation grade based on evaluation results and determine the fuzzy relationship matrix R. The membership is calculated as follows:

\(r_{i j}=\frac{n_{i j}}{n}(i=1, \ldots, 4, j=1, \ldots, 5)\)       (33)

R = (r_ij)_4x5       (34)

where r_ij represents the membership of the i^th evaluation indicator to the j^th evaluation grade, n_ij represents the number of doctors who make j^th evaluation grade to the i^th evaluation indicator, and n represents the number of doctors. 

(2) Establish Comprehensive Evaluation Matrix 

Determine the weight of each evaluation indicator by the expert physician according to subjective weight. The subjective weights are 0.2, 0.25, 0.25 and 0.3, respectively, and constitute a weight matrix W. Then a comprehensive evaluation matrix B is established based on the weight matrix W and fuzzy relationship matrix R. The comprehensive evaluation matrix B is calculated as follows:

W = (w₁,w₂,w₃,w₄) = (0.2,0.25,0.25,0.3)       (35)

\(B=W \circ R=\left(w_{1}, w_{2}, w_{3}, w_{4}\right) \circ \left(\begin{array}{ccccc} r_{11} & r_{12} & r_{13} & r_{14} & r_{15} \\ r_{21} & r_{22} & r_{23} & r_{24} & r_{25} \\ r_{31} & r_{32} & r_{33} & r_{34} & r_{35} \\ r_{41} & r_{42} & r_{43} & r_{44} & r_{45} \end{array}\right)=\left(b_{1}, b_{2}, b_{3}, b_{4}, b_{5}\right)\)       (36)

where o denotes fuzzy operator, the weighted average operator M (g⊕ ) is selected as fuzzy operator in this paper, as follows:

\(b_{j}=\min \left\{1, \sum_{i=1}^{4} w_{i} r_{i j}\right\}, j=1, \ldots, 5\)       (37)

(3) Comprehensive Evaluation Result

We adopt the weighted average principle and normalizes the comprehensive evaluation matrix. Furthermore, we regard the numerical value {5,4,3,2,1} as the rank of evaluation grade, and then the rank is weighted to the corresponding evaluation indicator so that the final comprehensive evaluation score S can be achieved. The weighted average principle is computed as follows:

\(S=\frac{\sum_{j=1}^{5} b_{j}^{\beta} \cdot j}{\sum_{j=1}^{5} b_{j}^{\beta}}\)       (38) 

where β is the weighting coefficient, and β =1 is selected in this paper.

According to the above steps, we can obtain the comprehensive evaluation score of the interaction performance of simulation system based on the four models, as shown in Table 2.

Table 2. Comparison of the comprehensive interaction performance

Table 2 shows that compared to the other three models, our proposed model achieved the best final comprehensive evaluation score. It indicates that our proposed model has a better interaction performance in terms of naturalness, fluency, authenticity, and stability, which provides a virtual environment with a sense of realism.

4. Conclusion

In this paper, we proposed an optimized mass-spring model with shape restoration ability based on volume conservation to simulate the soft tissue of the virtual liver. The virtual liver deformation simulation system is built on PHANTOM OMNI force tactile interaction device with VC++ 2017 and OpenGL. Our system incorporates flexion spring into the MSM to restore its original shape in constructing the soft tissue surface model and adopts the particle swarm optimization algorithm to optimize the model parameters. Besides, the volume conservation constraint is considered into the PBD approach to maintain the volume of the virtual liver for constructing soft tissue volumetric model. Experimental results show that our proposed model provides a good shape restoration ability as well as incompressibility, enhances the accuracy and realism of deformation. 

Our proposed model simulates the linear elasticity of soft tissue. In future work, we will focus on other complex biological characteristics of soft tissue, such as non-linearity, viscoelasticity, etc. We believes that the deformation can be modeled much closer to the real liver under the same external force.