Soft Tissue Deformation Model Based on Marquardt Algorithm and
Enrichment Function

In order to solve the problem of high computing cost and low simulation accuracy caused by discontinuity of incision in traditional meshless model, this paper proposes a soft tissue deformation model based on the Marquardt algorithm and enrichment function. The model is based on the element-free Galerkin method, in which Kelvin viscoelastic model and adjustment function are integrated. Marquardt algorithm is applied to fit the relation between force and displacement caused by surface deformation, and the enrichment function is applied to deal with the discontinuity in the meshless method. To verify the validity of the model, the Sensable Phantom Omni force tactile interactive device is used to simulate the deformations of stomach and heart. Experimental results show that the proposed model improves the real-time performance and accuracy of soft tissue deformation simulation, which provides a new perspective for the application of the meshless method in virtual surgery.


Introduction
In medical education, interns learn from experienced experts to improve their expertise. Traditional surgical training is often carried out on the remains of animals or humans, however, for it there still exist defects such as long training period, non-reusability, and high cost. With the development of computer science and technology, emerging virtual surgery simulation system has gradually been able to solve the above problems [1,2].
In the virtual surgery simulation system, interns operate on the virtual soft tissue model repeatedly through the force and tactile feedback device to improve their surgical skills. Therefore, the soft tissue model is required to be accurate and with good real-time performance, and be able to handle topological changes like cutting and suturing [3][4][5][6]. Common soft tissue modeling methods include mass-spring model [7][8][9][10], finite element model [11][12][13] and meshless model [14][15][16]. The structure and calculation of the mass-spring model is simple, but its parameters are difficult to select and the iterative calculation is unstable [17][18][19]. Kot et al. [20] proposed the mass-spring model with adjustable Poisson's ratio, the model does not need to introduce additional structures and can freely represent homogeneous isotropic materials, but the stability is poor. Li et al. [21] proposed an improved surface mass-spring model, which develops new flexion springs and surface triangle topological element to improve the deformation accuracy and calculation efficiency, but its parameters are difficult to select and determine. The finite element model has high precision and adaptability, but it is difficult to realize real-time simulation due to its high computational complexity [22][23][24]. Wang et al. [25] simulated subcutaneous adipose tissue based on a linear elastic and hyper-plastic finite element model, this model has high simulation accuracy, but the computational complexity is high and the real-time performance is poor. Marinkovic et al. [26] proposed a geometrically nonlinear corotational finite element model to simulate liver and stomach, which has high accuracy but poor real-time performance. Compared with the previous two models, the meshless model only needs a set of discrete nodes and does not need to process the mesh data in advance, so some problems such as mesh distortion and mesh entanglement will not appear, and it has great advantages in processing large deformation of soft tissue [27][28]. Belinha et al. [29] combined the natural neighborhood radial point interpolation method with fracture growth algorithm to simulate fracture propagation in brittle materials, which is faster and more efficient than the mesh-based method. Dehghan et al. [30] used element-free Galerkin method to simulate the invasion of cancer cells into surrounding tissues, and the effect of the simulation is realistic. Cheng et al. [31] proposed a novel interactive meshless cutting model, which divided the model cutting into three stages with high simulation accuracy. Zhou et al. [32] proposed a meshless local radial point interpolation method for soft tissue modeling, which solved the difficult problem of solving partial differential equations in traditional soft tissue models and improved computational efficiency. Although the meshless model has the above advantages, it still has some deficiencies, for example, it cannot accurately describe the biomechanical characteristics of soft tissue and has poor real-time performance during soft tissue simulation [33]. To solve above challenges, this paper proposes a new soft tissue deformation model based on Marquardt algorithm and enrichment function. The model simulates the deformation of soft tissue by introducing the element-free Galerkin method with fast convergence. Its shape function is constructed based on the moving least square method [34,35], and the Kelvin viscoelastic model and adjustment function are integrated to improve the simulative realism. Meanwhile, Marquardt algorithm is applied to fit the mathematical relation in advance between the force on the soft tissue and the displacement of each node. When the soft tissue during virtual surgery is subjected to the action of force, the deformation can be immediately shown by calling the fitting relation. The enrichment function is applied to deal with the discontinuity that occurs when soft tissue is used to perform interaction operation.
The organization of this paper is as follows: Section 2 introduces the nonlinear viscoelastic meshless method. Section 3 describes the fitting relation between the force and the displacement of each node based on Marquardt algorithm. Section 4 introduces the processing method for soft tissue discontinuity. Section 5 gives the experimental process and results. In the last section, we make a summary of our work.

Nonlinear Viscoelastic Meshless Method
The meshless method is different from the mesh-based method in that it only relies on a group of nodes, which the nodes are uniformly or randomly distributed on the internal and boundary, so the meshless method is suitable for large deformation and progressive cutting of soft tissue. The flow chart of the proposed method in this paper is shown in Fig. 1.
The deformation of soft tissue in the meshless method is approximately expressed by field variable u. Assuming that u h X ð Þ is the moving least squares approximation of the field function u X ð Þ at position X in the domain of the problem , and u X ð Þ can be expressed as: where P T X ð Þ represents the polynomial basis function vector and it is usually composed of single algebraic terms in three-dimensional space, let P T X ð Þ ¼ 1; x; y; z ½ ; and a X ð Þ represents the nonconstant coefficient vector.
To solve the nonconstant coefficient vector a X ð Þ, we consider the weight function w i X ð Þ ¼ w X À X i ð Þ of each node and minimize the Euclidean norm of approximation error we can obtain a X ð Þ as follows: where A À1 X ð Þ is the inverse of the weighted instantaneous matrix A X ð Þ, A X ð Þ is defined as: where W X ð Þ and P X ð Þ represent the weighted function matrix and polynomial basis function matrix, respectively. W X ð Þ and P X ð Þ are defined as follows: The moving least square approximation shape function È X ð Þ is defined as: Finally, the field function can be expressed as: The smoothness of weight function w s i ð Þ plays an important role in the continuity of shape function [36]. In general, cubic spline weight function is adopted, which is expressed as follows: where is the radius of influence of the dimensionless weight function, and q i indicates the radius of influence sphere of the point i.
In order to improve the authenticity of the simulation and present the biomechanical characteristics of the real nonlinear viscoelasticity of soft tissue, the Kelvin viscoelasticity model is applied to construct the mechanical model of soft tissue because its two springs and a damper can be used to represent linear viscoelastic properties [37], at the same time, an adjustment function is integrated to the viscoelastic relaxation constitutive relation to show the nonlinear relation between stress and strain. At last, we establish the nonlinear viscoelastic model and integrate it into the element-free Galerkin method to solve the equation. The Kelvin viscoelastic model is shown in Fig. 2.
In Fig. 2, r, r d , g, K 1 and K 2 represent the stress, stress time derivative, damper, the stiffness of the first spring and the stiffness of the second spring in the model, respectively.
The relaxation constitutive relation of soft tissue under stress loading in the Kelvin viscoelastic model can be expressed as: where r t ð Þ, E t ð Þ and e represent relaxation response, relaxation modulus and strain of the material, respectively. In order to make the soft tissue model has nonlinear characteristics, an adjustment function N e ð Þ ¼ e þ ke 2 is integrated to the viscoelastic relaxation constitutive relation, so the nonlinear viscoelastic relaxation constitutive relation is defined as: We divide the deformation simulation time T into n time slices t 1 ; t 2 ; . . . ; t n , Dt ¼ T=n represent a time increment. Du n , Dr n and De n represent the increment of displacement, stress and strain during the time period from moment t n to moment t nþ1 , respectively. We can obtain the stress of the Kelvin viscoelastic model at moment t n and t nþ1 through the follow equations: The incremental form of the nonlinear viscoelastic model is expressed as follows: where E k and r 0;n . represent the nonlinear relaxation coefficient and the initial stress, respectively.
At moment t nþ1 , displacement, stress and strain are expressed as follows: The relation between strain and displacement can be expressed as follows: where L represents partial differential operator, which can be defined as: Since the function obtained by the least square approximation is a smooth curve and does not pass the node value, the moving least square approximation shape function cannot characterize the Kronecker delta function. In order to enforce the essential boundary conditions, we use the weak form of the constrained Galerkin method as follows: where , d, u, D, u T , b, À t , t, À u , u and a represent analysis domain, Kronecker function value, displacement, elastic constant matrix, displacement vector, physical vector, natural boundary condition, corresponding surface force, essential boundary condition, displacement corresponding to essential boundary condition, and penalty factor, respectively.
By substituting Eqs. (13)- (17) into Eq. (19), the incremental form of nonlinear viscoelastic meshless solution equation can be obtained: where K n represents the viscoelastic stiffness matrix, it is defined as: where B i represents the matrix of the derivative of the shape function of node i, and it can be expressed as: where Φ i;x , Φ i;y and Φ i;z represent the derivative of the shape function with respect to x, y and z at node i, respectively. Similarly, we can get B j .
K a n is the penalty stiffness matrix determined by the shape function and the derivative of the shape function in the element-free Galerkin method, and DR n is the vector of unbalanced force. They are expressed as follows: Given the material parameters of the soft tissue model and time step, we can calculate Du n according to Eq. (20). Finally, the new displacement, stress and strain at each point can be figured out.

Relation Fitting of Force and Displacement Based on Marquardt Algorithm
Compared with the method based on mesh, meshless method does not need to maintain the topological information between data points, and can avoid complex topological problems of the model, but its computation is time-consuming. Therefore, we calculate the displacement of all nodes according to the nonlinear viscoelastic meshless solution equation. Then, we use the Marquardt algorithm to fit the relation in advance between the force applied to the soft tissue and the displacement of each node [38,39], and we can obtain the fitting surface to show the overall deformation of the model, which improves the real-time simulation.
In order to establish the relation between the applied force and the deformed surface, as well as the fitting relation between the applied force r and the stressed node d, it is assumed that n spatial forces r are applied to the node which can be decomposed into component forces r x , r y , r z in three directions. The corresponding displacement components at the stressed node d are Du x , Du y , Du z and the corresponding displacement function of each component is: . . . ; l z þ1 ½ are the system parameter corresponding to three components. x , y and z are the highest order of independent variables of the system which are determined by cross validation.
When a node is subject to the force, neighbor points will also move accordingly. Therefore, the relation between the displacement of the stressed node and that of other nodes is established to represent the deformation surface. For simplification, we assume that the surface nodes have the same coordinates z. Under the force r z , the induced surface function is: where Dx, Dy, and Dz represent the induced displacement of other points in three directions under the force r z Á c i1 ¼ c i11 ; c i12 ; . . . ; c i1a 1 þb 1 þ1 ½ , c i2 ¼ c i21 ; c i22 ; . . . ; c i2a 2 þb 2 þ1 ½ and c i3 ¼ c i31 ; c i32 ; . . . ; c i3a 3 þb 3 þ1 ½ ði ¼ 1; 2; …; nÞ represent the parameters of three surface deformation functions. a 1 , a 2 , a 3 , b 1 , b 2 and b 3 are the order surface deformation which are determined by cross validation. Similarly, the induced surface functions of component force r x and r y can be obtained.

Discontinuity Processing in Meshless Method
In virtual surgery, interactive operation such as cutting and suturing on soft tissues is necessary. Thus, the weight enrichment function is introduced to deal with discontinuous cracks on the tissue surface caused by cutting and other operations [40]. To this end, we take the cutting treatment as a piecewise linear segment, and calculate the absolute distance between meshless nodes and the segment. According to the distance field, we obtain the enrichment function. Let the enrichment function multiplies the weight function of the corresponding node to replace the original weight function. In order to reduce the calculation cost, only the shape functions of nodes whose distance is less than their expansion parameters q i are recalculated. The shape functions of other nodes remain unchanged. Fig. 3 is the schematic diagram of the soft tissue with discontinuous cracks.
In Fig. 3, s, f, À u , u, À t and r represent normal coordinates, local coordinates, essential boundary conditions, displacement boundary conditions, natural boundary conditions and force boundary conditions, respectively.
For a given point x; y ð Þ, its two-dimensional distance function d 2 x; y ð Þ can be calculated as follows: where d s f ð Þ is a one-dimensional signed distance function with endpoints f 1 ; f 2 ð Þ.
We solve the partial derivative of the distance function d 2 x; y ð Þ with respect to the normal coordinate s, and we can obtain the discontinuous function ' x; y ð Þ as follows: Finally, enrichment function h x; y ð Þ can be expressed as:

Data Acquisition
Meshless method is used for soft tissue modeling, which relies on point cloud data. Therefore, CT images of soft tissue are imported into the software of Mimics to generate 3D images, and then the STL format files of the 3D model are exported. The software MeshLab is used to convert STL format files into OBJ format files to obtain vertex information and point cloud data. And the CT data in this paper is provided by the First Affiliated Hospital of Nanjing Medical University and the consent of patients is obtained.

Experimental Environment
In this paper, Sensable Phantom Omni hand controller is used as force and tactile interactive device, other hardware includes a computer with Intel Core i9-7960X CPU, NVIDIA GeForce 1080Ti graphics  Fig. 4 shows the simulation environment.

Deformation Simulation of Soft Tissue
During soft tissue deformation simulation, the tissue surface is considered as the problem domain with scattered nodes, and its displacement can be calculated by meshless method. The proposed model in this paper is used to establish the fitting relation between forces and soft tissue surface deformation, and the original displacement and induced displacement can be obtained through this relation when different forces are applied to the nodes. To verify the validity of the proposed model, simulating the cutting operation based on the model on stomach and heart is performed. As shown in Figs. 5 and 6, surgical instruments are used to perform incisions on stomach and heart with forces of 0.6 N, 1.2 N, 1.8 N and 2.4 N. By observing the deforming effect of stomach and heart under different cutting forces, we can see that the incision of the soft tissue model is smooth. At the same time, according to the surgeon's feedback, we can know that the proposed model can effectively simulate the deformation responses of human soft tissue to some extent. And in the future research, we will compare the deformation data of the proposed model with the deformation data of real human soft tissue and improve it so as to make the proposed model more real for human soft tissue deformation simulation and apply it to the interactive training of virtual surgery.

Accuracy Verification
In virtual surgery simulation, the accuracy of the model is crucial. Since the accuracy of the finite element model is high, it is used as a reference model. Under the same force, we calculate the deformation displacement of the reference model, as well as the deformation displacement of the proposed model and the interactive meshless cutting model [31]. The specific process is as follows: Firstly, the finite element model, the proposed model and the interactive meshless cutting model are used to establish three lung models, and the same stress points are selected. And the intention of modeling the lung model in accuracy verification is to show the generality of the proposed model. Secondly, choose 10 points near the stress points of the models as sample points, applying cutting force to them and make the models have small deformation (no cracks) and large deformation (to generate cracks). Finally, the coordinates of the sample points after deformation of each model are recorded, and the deformation displacement of the three models are calculated by using Eq. (33), respectively.
where x s i ð Þ, y s i ð Þ and z s i ð Þ represent the coordinates of the sample points i after the deformation of the three simulation models in x, y and z directions, respectively; x 0 i ð Þ, y 0 i ð Þ and z 0 i ð Þ represent the initial coordinates of the sample points i in x, y and z directions, respectively.
We can see that Figs. 7 and 8 are deformation displacement radar diagrams of the three models under small deformation and large deformation, respectively.
At the same time, we calculate the RMSE (Root Mean Square Error) values of the deformation displacement of the proposed model and the interactive meshless cutting model. Note that the RMSE values are obtained based on the reference model.
where n represents the number of sample points, x s i ð Þ, y s i ð Þ and z s i ð Þ represent the coordinates of the sample points i after the deformation of the proposed model and the interactive meshless cutting model in x, y and z directions, respectively; x c i ð Þ, y c i ð Þ and z c i ð Þ represent the coordinates of the sample points i after the deformation of the reference model in x, y and z directions, respectively. Tab. 1 shows the RMSE of It can be seen from Figs. 7, 8 and Tab. 1 that the deformation displacement of the proposed model is more similar to the reference model, and the RMSE values of the proposed model is smaller regardless of small deformation or large deformation, which verifies the accuracy of the proposed model in this paper.

Real-Time Verification
In order to verify the real-time performance of the proposed model, we use the local radial point interpolation model [32] with better real-time performance in virtual surgical training and the proposed The finite element model The proposed model The interactive meshless cutting model  model to establish the stomach model, and we select 8 groups of spatial forces within the range of 0-4.0 N and apply them to the same stress point. After that, we compare the time consumed by the two models in calculating the deformation results. As shown in Fig. 9, the time consumed by the two models under different forces is recorded as histogram. It can be seen that the time consumed by the proposed model in calculating the deformation results is shorter and it has better real-time performance. Therefore, the proposed model meets the real-time requirement of virtual surgical training.

Multi-Performance Indicators Analysis
In order to verify the superiority of the multi-performance indicators of the proposed model, we invite 28 doctors from the First Affiliated Hospital of Nanjing Medical University to evaluate the heart models which are established by the proposed model, mass-spring model and finite element model. There are 4 chief doctors marked D a , 7 associate chief doctors marked D b , 10 attending doctors marked D c and 7 interns marked D d . According to the professional level of doctors, we set their score weights as 0.35, 0.3, 0.25 and 0.1. The evaluation indicators include force feedback ability I 1 , visual fluency I 2 , deformation efficiency I 3 , immersion I 4 , texture characteristics I 5 , interaction naturalness I 6 and system stability I 7 . The score of each indicator is within the range of [0,10]. The score Sc ij ði ¼ 1; 2; 3; j ¼ 1; 2; 3; …; 7Þ of 7 indicators of the three models can be calculated as follows: The scoring results of Sc ij ði ¼ 1; 2; 3; j ¼ 1; 2; 3; …; 7Þ are shown in Tab. 2.
We use the comprehensive evaluation method of RSR (Rank Sum Ratio) to analyze each indicator score of the three models [41], and then we can determine the overall performance of the three models. The details of the calculation procedure are as follows.  Step 1: rank the indicators to be evaluated. Since the 7 indicators evaluated by the three models are all high-quality indicators, we rank the indicators to be evaluated according to the size of the indicator scores, and the result is shown in Tab. 3.
Step 2: calculate the rank sum ratio RSR i of each indicator. And RSR value is the average rank of seven indicators of each model. The larger the RSR value is, the better the indicators' performance of the model will be. The equation is as follows: where R ij , m and n are the rank, number of indicators and number of groups of each indicator, respectively.
The RSR values and ranking of model evaluation indicators of three models are shown in Tab. 4.
Step 3: determine the RSR distribution. We rank the RSR values from small to large and list the frequency f of each group. Meanwhile, we determine the rank range R and average rank of each group R by calculating the downward cumulative frequency P f , and calculate the cumulative frequency R=n À Á Â 100%. Finally, the Probit value of the corresponding probability unit is obtained by referring to the comparison table of percentages and probability units, as shown in Tab. 5.
Step 4: calculate the linear regression equation RSR ¼ a þ b Â Probit. We set the probability unit value Probit i corresponding to the cumulative frequency as the independent variable and the rank sum ratio RSR i as the dependent variable, the linear regression equation is calculated as:   Step 5: grade the three models according to the RSR estimate corresponding to the regression equation. We divide the three models into three grades according to the grading table: excellent, good and poor. And the grading table is based on the grading standards provided by statistician professor Tian. The grading is shown in Tab. 6 below.
Tab. 6 shows that the proposed model in this paper has better overall performance.

Conclusion
In virtual surgery, the real-time performance and accuracy of the soft tissue model are crucial. Therefore, this paper proposes a soft tissue deformation model based on Marquardt algorithm and enrichment function. The model uses the nonlinear viscoelastic meshless method based on element-free Galerkin to simulate the soft tissue deformation. At the same time, Marquardt algorithm is applied to fit the relation between the force and displacement of surface deformation on soft tissue, which improves the authenticity and real-time performance of simulation. Furthermore, the enrichment function is applied to deal with the discontinuity of soft tissue in the interactive operation, which makes the interactive simulation between soft tissue and surgical instruments more stable and effective.
To verify the validity of the model, gastric and cardiac dissection operations are simulated. Finite element model is regarded as a reference model, the RMSE of deformation displacement of the proposed model and the interactive meshless cutting model based on the reference model are compared, so as to verify the accuracy of the model. In terms of the real-time performance of the model, the time consumed by the local radial point interpolation model and the proposed model in calculating the deformation is compared. In addition, the comprehensive evaluation method of RSR is used to evaluate the proposed model, the mass-spring model and the finite element model. The experimental results show that the deformation of the proposed model is valid and the accuracy is equivalent to the finite element model. Moreover, the model has good real-time performance, and nice overall performance in force feedback ability, visual fluency, deformation efficiency, immersion, texture characteristics, interactivity and system stability.
Acknowledgement: The authors would like to thank the editors and two anonymous reviewers for their constructive comments on our manuscript, as well as the 28 doctors from the First Affiliated Hospital of Nanjing Medical University for their contributions to the experimental verification of our manuscript, which makes our work much better.