Integrating viscoelastic mass spring dampers into position-based dynamics to simulate soft tissue deformation in real time

We propose a novel method to simulate soft tissue deformation for virtual surgery applications. The method considers the mechanical properties of soft tissue, such as its viscoelasticity, nonlinearity and incompressibility; its speed, stability and accuracy also meet the requirements for a surgery simulator. Modifying the traditional equation for mass spring dampers (MSD) introduces nonlinearity and viscoelasticity into the calculation of elastic force. Then, the elastic force is used in the constraint projection step for naturally reducing constraint potential. The node position is enforced by the combined spring force and constraint conservative force through Newton's second law. We conduct a comparison study of conventional MSD and position-based dynamics for our new integrating method. Our approach enables stable, fast and large step simulation by freely controlling visual effects based on nonlinearity, viscoelasticity and incompressibility. We implement a laparoscopic cholecystectomy simulator to demonstrate the practicality of our method, in which liver and gallbladder deformation can be simulated in real time. Our method is an appropriate choice for the development of real-time virtual surgery applications.


LX, 0000-0003-4326-181X
We propose a novel method to simulate soft tissue deformation for virtual surgery applications. The method considers the mechanical properties of soft tissue, such as its viscoelasticity, nonlinearity and incompressibility; its speed, stability and accuracy also meet the requirements for a surgery simulator. Modifying the traditional equation for mass spring dampers (MSD) introduces nonlinearity and viscoelasticity into the calculation of elastic force. Then, the elastic force is used in the constraint projection step for naturally reducing constraint potential. The node position is enforced by the combined spring force and constraint conservative force through Newton's second law. We conduct a comparison study of conventional MSD and position-based dynamics for our new integrating method. Our approach enables stable, fast and large step simulation by freely controlling visual effects based on nonlinearity, viscoelasticity and incompressibility. We implement a laparoscopic cholecystectomy simulator to demonstrate the practicality of our method, in which liver and gallbladder deformation can be simulated in real time. Our method is an appropriate choice for the development of real-time virtual surgery applications.

Introduction
Surgery simulators can help novices become familiar with real surgical procedures without human loss. They enable scientific and repeatable training through visual rendering via software and interactive feedback via a hardware-based manipulator. Compared with actual surgical training, simulators reduce the cost of learning and guarantee the training's effectiveness [1]. Simulating the deformation of soft tissue is necessary for rendering three key aspects of a surgery simulator: the surgical environment, organs' interactions with the manipulator and force feedback [2]. A well-balanced deformation method can more closely simulate clinical manifestations. Further, accurate simulation of deformation, tearing, bleeding and force feedback make training more efficient. Modelling soft tissue as a viscoelastic body represents more complex nonlinear behaviours than a simple elastomer and viscous body, including geometric features such as creep, stress relaxation and incompressibility [3]. Therefore, it is crucial to develop a stable and fast method with which to simulate the complex deformation behaviours of biological soft tissue in real time.
The development of a real-time deformation algorithm has been accompanied by improvements in computing hardware. In early work, poor computational performance required the use of geometric techniques to simulate deformation, including spline and patch, freeform deformation and skin-based animation. In these techniques, physical accuracy is sacrificed for computational efficiency, and the system has no knowledge about the material being deformed [4]. With advances in computing power, two types of physics-based approaches have flourished. One solves partial differential equations based on a continuum model. This approach is typified by a finite-element method based on meshing object space, and uses the same strategy as the finite-volume and finite-differential methods [5]. The second approach, a meshless method based on shape functions and support domains has been used increasingly widely, such as with smoothed particle hydrodynamics and radial basis function methods [6]. The advantage of continuum methods is that they are sufficiently accurate to reflect an object's physical characteristics. They are widely used in structural analysis and electromagnetic radiation calculation [7]. Owing to deformation introducing the mesh distortion which reduced accuracy, finite-element method is naturally difficult to deal with large deformation. But the meshless method does not require connection information between simulation node, it works well under large deformation. Also, the material point method uses particles called material point to keep the physical properties and calculate the deformation gradient under background grid, it has been widely used in large deformation, multiphase simulation and fluid dynamics [8]. However, these methods are prohibitively time consuming, and are usually employed in offline applications. Our simulation aims to reflect the significant characteristics of soft tissue (such as viscoelasticity and incompressibility) in surgical applications; for real-time applications, it is appropriate to choose a physics-like method, using mass spring dampers (MSD) [9] and position-based dynamics (PBD) [10]. By abstracting typical features into this model, we achieve high computational efficiency.
Unlike the finite-element method, which uses the discrete solution domain of a partial difference equation, the mass spring method models discrete objects directly. An object consists of node masses with no size and elastic springs with no mass. Deformation is generated by applying elastic force to a node. The state of the system can then be determined by solving algebraic equations derived from Hooke's Law and Newton's second law. The mass spring method is simple to calculate and has a clearly structured topology; however, because nodes can only affect adjacent nodes, force propagation is limited and somewhat delayed, which causes local superelasticity [11]. A volume mesh topology and nonlinear spring can mitigate this effect [12]. When updating the node position, numerical integration is required for iteration. An explicit Euler step would cause the result to fail to converge, owing to second-order error. Using implicit Verlet or Runge-Kutta integration can improve solution accuracy [13]. Deformation in the model is directly related to spring stiffness, but the physical properties only indicate blurry correlation with spring stiffness. Therefore, it is difficult to adjust the parameters of MSD, which is a pressing concern for their application [14].
As with MSD, PBD discretizes an object as a set of node masses (figure 1), although the relationship between nodes is described as a constraint function. The constraint function is a scalar function that takes node positions as variables. When the constraint equilibrium state changes, node position is modified via constraint projection iteration to satisfy the constraint function [15]. Therefore, deformation calculation is a constraint-function optimization problem. Different constraint functions enable the realization of multiple deformation features. For example, distance and dihedral angle constraints are used to simulate cloth [16]. Strain energy constraint is used to simulate soft bodies [17]. A Darboux frame is used to simulate an elastic rod [18]. Depending on the Gauss-Seidel method used to directly modify the node position, the PBD method is unconditionally stable. Despite its many advantages, the conventional PBD method has two notable problems. One is that constraint stiffness and deformation effects depend on iteration count and timestep, while the iteration count determines the constraint convergence and computation costs. Hence, it is necessary to decouple these parameters to more easily control the model.  The second problem is that no force or time parameters participate in the constraint projection step. Thus, it is difficult to simulate the viscoelasticity and nonlinearity of soft tissue, as these are time-related force characteristics. The distance-dependent variable constraint parameter can be used to fit the nonlinearity by controlling the node position directly [19]. However, due to the lack of time and force parameters in constraint projection step, the viscoelasticity cannot be integrated into.
This study proposes a new method for integrating viscoelastic and nonlinear springs into PBD. It realizes fast and stable simulation of rich and complex deformation behaviour, such as the viscoelasticity, nonlinearity and incompressibility of biological soft tissue. In this method, deformation effects do not depend on the iteration count and timestep. The contributions of this study are as follows: -a new method of integrating viscoelastic and nonlinear springs into PBD; -examples that demonstrate the validity of the proposed method compared with analytical solutions of viscoelasticity; -a simulation of liver and gallbladder deformation, demonstrating the different structures of the organs; -a laparoscopic cholecystectomy simulator that demonstrates the practical applicability of the proposed method.

Integration method
The objects that we aim to simulate are each composed of a set of mass points, constraints and springs. The goal of deformation calculation is to find the actual position of every mass point at each step. The position of a mass point is controlled by its total force load. For our method, based on Newton's second law of motion To determine f constraint , constraint energy potential is defined as potential, the force generated by the energy potential is in the direction of the maximum energy change to the mass point position. Thus, the constraint force is defined as the gradient of the mass point position of energy potential: To solve the mass point position state at each timestep, we discretize equation (2.1) with timestep t, using n as an index of the current timestep: Similar to extended position-based dynamics (XPBD) [20], we introduce a Lagrange multiplier to decompose the force into its directional and scalar components. Further, we use the time interval and the inverse stiffness of constraint together to obtainα = α/ t 2 for simplifying the description The right-hand side of the equation, which represents the predicted position, can be simplified as and is calculated prior to the constraint projection step. To determine the mass point position at timestep n + 1, we must solve the following constraint-functionbased minimization problem. and Analysing the above equations, our goal is to determine λ and x for timestep n + 1. The number of functions is equal to the number of variables, so there is a unique solution. Because the constraint function is abstracted on the physical deformation behaviour, it is often not the linear form of mass point position. There is no analytical solution for mass point position directly. To solve these nonlinear equations, we linearize them with the Newton-Raphson method and solve them iteratively. We set the iteration count to i and eliminate the n + 1 timestep: Our goal is to obtain x and λ in every iteration, ultimately using these values to update λ i+1 = λ i + λ and x i+1 = x i + x. This system requires two approximations (as does XPBD) that do not significantly influence the final solution. First, consider that when iteration count i = 0, g(x i , λ i ) = 0, and for large iteration counts, g(x i , λ i ) remains small, making g(x i , λ i ) negligible. Second, it is expensive to determine the second-order partial derivative of the constraint function: whereas rejecting the first part of ∂g/∂x has little influence on the convergence speed of the result. Therefore, it is reasonable to let ∂g/∂x = −M. For more detail, see [20]. Next, we can solve the linearized equations (2.9) and (2.10): and Then, we get x = M −1 ∇C(x i ) T λ from equation (2.12), and substituting it into equation (2.13), get Once λ is solved, x can be determined easily.  (2.15) and f damper as This direct velocity damp force to the point mass mimics viscoelasticity [21]. The spring elastic force is defined as k ij corresponds to l ij = ||x i − x j || − l 0 ij . Therefore, spring stiffness is a third-degree polynomial at low displacements, and linear at higher displacements.
The concept of our new method is to separate the behaviour model into two components: nonlinear viscoelastic MSD reflect the general physical characteristics of soft tissue; and XPBD handles complex effects, such as incompressibility and maintains the models' stability. Control parameters are listed in table 1 for clarification.

The algorithm
Using the derivations above, we can summarize the flow of our algorithm as follows. First, nonlinearity and viscoelasticity are reflected in the spring force calculation step and the extended spring force is introduced as external force to estimate the mass position. Then, the constraint projection loop step begins to modify the position to satisfy the constraint-function equation. We use XPBD to keep the deformation effect from depending on the iterator count. We benefit from using the Gauss-Seidel method, which maintains the stability of the final mass point's position.
From the above steps, compared to conventional PBD, the only part that is added is the spring force calculated in the predicting position step. Thus, the time cost per timestep of this method is slightly higher than PBD. However, it is still acceptable and can meet the demand of realtime deformation applications, while achieving unconditional stability and modelling more complex deformation behaviours.

Collision detection and response
A precise and fast collision detection method is crucial for avoiding penetration and triggering the deformation calculation, lest the trainee's immersion and concentration be disrupted. For this purpose, we use the modified space hash algorithm [22]. As shown in figure 2, in our model, the simulated objects are constructed as triangles and spheres.

1:
for all particles i do 2: initialize for all springs j do 4: initialize l restlength , b 0 , b 1 , k 1 , k 2 5: for all constraints k do 6: initialize α k and other constraint specific parameters 7: loop (timestep = t, timestep count = n + 1) 8: for all particles i do 9: predict position for all particles i do 11 generate collision constraint C coll (α = 0) 12: for all constraints k do 13: reset λ k = 0 14: loop max iteration count times 15: solve constraint function C to get λ k , x 16: update for all particles i do 18: update    to the axial plane in the x-, y-and z-directions. We maintain a two-dimensional array for each primitive detection state to avoid repeat detection. After collision detection is completed, collision response moves the primitives to a collision free state. In conventional PBD, two types of distance constraints are generated for collision response.
Sphere-sphere collisions, as in figure 3, can be handled as a vertex distance constraint. For sphere centre points q and p, and radii r 1 and r 2 , the constraint function is (2.20) Sphere-triangle collisions and triangle-triangle collisions, demonstrated in figure 4, can be handled as a vertex and triangle distance constraint. The triangle-triangle collisions are resolved using three separate vertex and triangle distance constraints. For sphere centre point q; radius r and triangle points q 1 , q 2 and q 3 , the constraint function should be After collision constraints are generated, they are added to later constraint projection steps for collision response calculation. For all collision constraint, α = 0 to guarantee that the constraint equation is fulfilled as PBD.

Model example
There is a vital link between the structure of soft tissue and organ deformation behaviour. We can naturalistically simulate effects involving an organ's topological structure by introducing corresponding geometric constraints on mass point position. This technique enables our method to quickly and stably simulate complex behaviours. As an example, we demonstrate two typical organ models: the liver and gallbladder. The liver is an entity organ, meaning that it is composed of closely arranged cells. Conversely, the gallbladder is a cavity organ, composed of stratiform distribution surface cells. We use a closed triangular surface mesh to construct the topology of the gallbladder, and a tetrahedral volume mesh to construct the topology of the liver. (d-f ) geometry data, surface mesh and physical data of the gallbladder. Figure 6. Overstretching constraint enforces p 1 and p 2 during compression (a) and stretching (b), respectively.
In more detail, model expression is separating into three steps (figure 5). First, the organ surface mesh models are constructed from the 'Virtual Chinese Human' colour slice dataset [23]. Second, for entity organs such as the liver, we generate a tetrahedral volume mesh. Finally, the basic elements for deformation calculating are generated from topology information such as nodes, springs and constraints. Different volume constraints are enforced for volume preservation when modelling the structural differences in the liver and gallbladder. Given the local volume preservation of liver, tetrahedral volume constraint is applied: Similarly, to model the volume preservation of gallbladder, a closed surface mesh volume constraint is applied, as below: (2.23) In the viscoelasticity calculation step of the mass spring model (MSM), there would be superelasticity effects on the spring that load a large force near the contact point ( figure 6). To mitigate these effects, we introduce the below overstretching and compression constraints:

Viscoelasticity and nonlinearity
Soft tissue, as a polymer, has nonlinear, creeping and stress-relaxing characteristics. To verify that the proposed integration model can represent these features, we built a simple cube model with 60 mm long edges for testing ( figure 7). Each direction is composed of six cube primitives. The centre point of the upper surface was selected as a test point. The corresponding force load and displacement constraints are applied to it. First, we tested the nonlinearity of the model. We applied a constant displacement velocity constraint to the test point. The direction of velocity was perpendicular to the upper surface. Then, through the force balance state of the test point, we obtain the force state. Figure 8 indicates that while the displacement is small, the relationship between displacement and force is nonlinear, and then it gradually transitions to a relatively linear relationship. When the displacement triggers the overstretching constraint, the slope of the figure is flat. Meanwhile, as the displacement velocity increases, the force increases. This phenomenon is consistent with the nonlinear characteristics of soft tissue.
Second, we tested creeping (figure 9). We still applied a constant force to the test point perpendicular to the upper surface. For the same values of b 1 , the displacement of the test point exhibits nonlinear growth over time. As b 1 increases, the displacement of the test point is smaller, and finally, stabilizing to equal the displacement sample. This phenomenon is consistent with the creeping characteristics of soft tissue.
Third, we tested stress relaxation ( figure 10). Before force sampling, the test point is moved to a same displacement of 20 mm over 0.5 s. Then, keeping the test point fixed, its force is sampled for 2 s. In the beginning of the test, due to the momentary displacement applied on test point, the instant elastic force loaded on the connected mass points around test point are much larger than the damping force, it causes the slight oscillatory response of the total loading force. But it quickly stabilized under the constraint of damping force by viscoelasticity. With the increase of the viscoelasticity coefficient b 1 , the trend of force declining becomes slower and the oscillatory response becomes weaker. The oscillatory response also can be optimized by variable timestep simulation. On the general trend, the model maintains the tendency of forcing decline over time, and the decay rate becomes slower as the viscoelasticity coefficient b 1 increases. This phenomenon is consistent with the stress-relaxing characteristics of soft tissue.
It is not easy to find viscoelastic model parameter to directly correspond with the actual physical parameters. However, the effects of the model have the characteristic of real organs. Further, by adjusting the control parameters, the extent of the modelled organs' viscoelasticity changes consistently. Therefore, the proposing model meets the requirement of visual plausibility for visual surgery applications.    11g). In order to measure the global effect of each model, the whole surface node distance difference parameter Q is calculated, which is defined by the adding each node distance to the FEM standard.

Comparative analysis of liver model
Smaller Q demonstrates a smaller overall variance to the FEM model. It is found that the value of Q of our model is 2.57, that of PBD is 4.75 and that of MSM is 6.20. The above results show that our model has higher accuracy, while ensuring visual effects.

Volume preservation
The structural features of the entity organs and cavity organs can be expressed by two different volume constraints. For the volume preservation effects of entity organs, we apply a local volume constraint. The tetrahedral volume constraint only influences the node position of related tetrahedra. Here, we use the cantilever beam model for testing by detecting the volume change of a cantilevered beam as it sags under gravity ( figure 12). From the initial state (t = 0 s) to the natural sag (t = 1 s), the volume change ratio remained relatively low ( figure 13). For volume preservation in cystic organs, we apply a global volume constraint. The closed surface volume constraint influences the position of all surface nodes at the same time. In this study, we use a falling gallbladder model for testing, by detecting the volume change of a gallbladder falling and pressed by a plane ( figure 14). Owing to the enforcement of the global volume constraint, although the volume change ratio increases, it maintains an acceptable level ( figure 15).
The two volume constraints above can constrain the volume of the model within acceptable levels, thus reflecting the structural features of the corresponding organ.

Stability
It is critical that soft tissue deformation methods for real-time virtual surgery applications can meet stability requirements. If unrealistic organ explosion occurs during interaction between the organs and instruments, the trainee's experience and immersion will be poor. Numerical integration errors make it easy to trigger such phenomena in traditional mass spring methods. By integrating PBD constraints into the spring force solution, we improve the stability of the model. In figure 16, we use an abnormal

Virtual laparoscopic cholecystectomy application
To verify the practicality of our method, we constructed a laparoscopic virtual surgery application. We built physical models for two important organs, the liver and gallbladder, using the proposed method.     figure 17. Surgical scene information was recorded by the Unity3D frame rate profiler tool (table 3), indicating that a high frame rate (greater than 30 fps) was maintained during the operation. In figure 18, the average time cost per frame of three steps in figure 17 was recorded. From step (a) to step (b), and step (b) to step (c), the contacts between instruments and organs are increasing. Owing to collision response, the increase in constraint numbers leads to a slight increase in time cost. Nevertheless, the total time cost per frame is always below the real-time standard requirements (30 fps). Meanwhile, the fluctuation of time cost during the operation is small.

Conclusion
To simulate the viscoelasticity, nonlinearity and incompressibility of soft tissue deformation in real time, we have derived a new PBD method that integrates viscoelastic MSD. In this method, viscoelasticity and other time-related characteristics are simulated based on MSD. Then, the spring and external forces are combined with a constraint force generated by a PBD constraint function to correct the movement of particles. This method succeeds in controlling complex deformation behaviour through multiple constraint types, and is independent of iteration count and unconditionally stable. We generated a series of simple examples, which indicated our ability to control the extent of viscoelasticity through parameter adjustment. Finally, we applied this method to simulate soft tissue deformation in a laparoscopic cholecystectomy. This simulation showed the realistic deformation effects of the liver and gallbladder, and verified the practicality of the proposed method. The new method can meet the demand of modelling soft tissue deformation in real-time applications such as virtual surgery.
In future work, we plan to build a viscoelastic spring based on PBD through more controllable distance constraints. The XPBD framework has shown that the periodicity of a spring oscillator can be simulated by adjusting the parameters corresponding to the inverse stiffness of a realistic spring. Under extreme deformations, such as momentary displacement or force, the system has slight oscillations, which can be improved by introducing variable timestep simulation. Meanwhile, the variable timestep can be used for reducing time cost. Also, force feedback plays a significant role in measuring the effectiveness of the operation and enhancing the immersion of the trainers. Using the method proposed above, the feedback force in the deformation process can be easily solved by analysing the motion state of mass points. In the future application, we will add force feedback devices such as servo motors to produce reaction forces on the medical instrument and study the numerical accuracy of feedback force. Further, it would be worthwhile to build an overall relationship between model and physical parameters, as in the case of Young's modulus or Poisson's ratio.