1. A method for soft tissue modeling based on position constraints and non-linear springs, the method comprising the steps of:

step 1: constructing a soft tissue model for the acquired medical image information by utilizing a three-dimensional modeling technology, wherein the model is composed of tetrahedral units formed by nonlinear springs and virtual body springs, and the tetrahedral model structure is used for acquiring surface node information and model internal volume information;

step 2: simulating the deformation behavior of a soft tissue model, wherein in the deformation process of the model, the nonlinear spring enables the stress and the strain of the model to present a linear and nonlinear relation, the virtual body spring can inhibit the change of the volume of the model, and the virtual body spring can act together to simulate the biomechanical characteristics of the soft tissue;

and step 3: estimating the approximate position of each node through a kinetic equation, and then judging whether the state of a node constraint equation changes or not, wherein the node constraint comprises spring length constraint and spring curvature constraint; if the initial state of the equation changes, correcting the node to a correct position within a reasonable range; if the initial state of the equation is not changed, the nonlinear spring and the virtual body spring simulate the deformation effect of soft tissues;

and 4, step 4: and ending the single-step iteration of calculating the deformation, and entering the next round of circulation.

2. The soft tissue modeling method based on position constraint and nonlinear spring as claimed in claim 1 is characterized in that said step 3 specifically comprises:

3.1 for any node i in the model, calculating the resultant force received by the node and preliminarily estimating the position of the node;

the calculation formula of the resultant force received by the node i is as follows:in the formula M_iRepresenting the quality of the node i in terms of,represents the acceleration of node i, F (u)_i) Representing the sum of internal forces, F, experienced by node i_{ext_i}An external force indicative thereof;

the position expression of the node i is as follows: u. of_t+Δt＝u_t+v_i·Δt；

3.2 judging whether the position of the node i meets the constraint equation or not, and if the position of the node meets the constraint equation, reserving the estimated node position in 3.1; if not, correcting the node position according to a constraint equation; the constraint equations comprise spring length constraint and spring bending constraint equations;

the spring length constraint is suitable for the distance between adjacent nodes, and if the deformation rate of the spring length between the nodes is greater than the preset deformation rate, the positions of the nodes are directly corrected to meet the critical deformation rate;

if the length of the spring exceeds the threshold value and the state of the constraint equation changes or the bending degree of the deformation of the spring is larger than the critical value, the nodes p at the two ends of the spring_iAnd p_jMaking the corresponding corrections.

3. The soft tissue modeling method based on position constraint and nonlinear spring as claimed in claim 2 is characterized in that the approximate position estimation derivation process of each node in the step 3.1 is as follows:

F(u_i) The spring consists of a spring force and a spring damping force:

the expression shows that the internal force applied to the node i is equal to the sum of the spring force, the damping force and the volume force of all the adjacent nodes; wherein k is_ijRepresenting the elastic coefficient between node i and node j, c_ijRepresenting the spring damping coefficient between node i and node j,representing the first derivative of the displacement, i.e. the velocity, l, of node i_ij＝u_i-u_jRepresenting the distance vector of the spring between node i and node j,representing the initial distance of the spring between node i and node j;

the nonlinear spring force and its corresponding coefficient are expressed as follows:

in the formula, k_ijThe spring elastic coefficient connecting the node i and the node j; k is a radical of₁And k₂Is a constant; l_ijShowing the change in length of the spring, u_iIs the position of node i, | u_i-u_j| | is the distance between the node i and the node j; when the displacement is small, the relationship between the spring force and the length change is a cubic polynomial, and when the displacement is large, the relationship is a linear relationship;

the formula for the volume force is as follows:

in the formula, V_jRepresenting the volume, k, of tetrahedral units containing a node i₃Is the virtual body spring coefficient;andrespectively representing the current volume and the initial volume of the tetrahedral unit;

the damping force of the spring is used for simulating the viscoelasticity of soft tissues, and the formula of the damping force of the spring is as follows:

deriving a system of differential equations for node i:

and obtaining the position of the node i according to the differential equation system.

4. A method of soft tissue modelling based on position constraints and non-linear springs according to claim 2 wherein the rate of deformation of the spring length in step 3.2 is as follows:

wherein τ represents the deformation rate of the spring length, l_curAnd l_initRespectively representing the current length and the initial length of the spring; judging the degree of compression or extension of the spring through the deformation rate, and further using a constraint equation to generate position correction for the node;

for each node, the spring length constraint equation is:

in the formula, τ₀A threshold value representing the rate of spring length deformation;

if the spring length exceeds the threshold value and the state of the constraint equation changes, the node p at the two ends of the spring_iAnd p_jMaking corresponding corrections:

5. the soft tissue modeling method based on position constraint and nonlinear spring as claimed in claim 2 characterized in that the spring bending degree δ expression in step 3.2 is as follows:

wherein, p'_iAnd p'_jIs the initial position of the node at both ends of the spring, p_iAnd p_jRespectively are the current positions of the nodes at the two ends;

δ₀representing the critical value of the spring bending, the constraint equation of the spring bending is as follows:

if the bending degree of the spring deformation is larger than the critical value, the nodes p at the two ends of the spring_iAnd p_jThe position is revised to satisfy the constraint equation:

Background

Virtual surgery is an important surgical training aid system that can closely simulate the surgical environment in a low-cost mode. Virtual surgery not only can provide visualization requirements for users, but also enables doctors to become familiar with the use of surgical instruments and accumulates surgical skills. Currently, a perfect virtual surgery system should have three requirements at the same time: real-time, accuracy, stability. Soft tissue modeling is one of the core modules in virtual surgery, and the fineness of the model plays a decisive role in the reliability of the system.

The soft tissue modeling method mainly comprises a geometric modeling method and a physical modeling method. The geometric modeling method redraws a geometric model of soft tissue by three-dimensional reconstruction of medical images, and the physical modeling method builds a physical model according to biomechanics and kinematics characteristics of the soft tissue. The geometric modeling method is simple in calculation, but the constitutive relation of object materials is not considered, so that the accuracy of the model cannot be guaranteed (soft tissue deformation and tearing model research in brain surgery virtual surgery, Nanchang university, 2016). Currently, common physical modeling Methods include Finite Element Method (FEM) and Mass-Spring Method (MSM) (Soft Tissue Deformation and Optimized Data Structures for Mass-Spring Methods, IEEE International Conference on Bioinformatics & Bioinformatics. IEEE, 2009: 22-24).

FEM is a modeling Method based on continuous medium theory for discretizing a continuous domain problem into a plurality of Finite elements, The deformation of The model is solved based on biomechanics, and The soft tissue model based on FEM has higher accuracy (An Introduction to The Fine Element Method, Topics in Engineering materials, 1992 (81): 37-60). Therefore, the finite element method is widely applied to soft tissue modeling. Although the FEM can accurately simulate mechanical characteristics, the biological characteristics of human Soft tissues are generally complex and have fine internal structures, which causes the solving of the FEM to consume a large amount of resources and cannot meet the requirements of virtual surgery on Real-Time Interaction (A Surface Mass-Spring Model With New flexibility Springs and collagen Detection on Volume Structure Interaction, IEEE Access, 2018, 6: 75572-75597).

Unlike the continuous physical model, MSM discretizes soft tissue Models into multiple particles, and the interaction between the particles is achieved by springs connecting the particles (Calibration of Mass Spring Models for organic Simulations, IEEE/RSJ International Conference on Intelligent Robots & systems. IEEE, 2007: 370-375). Compared with FEM, MSM has high calculation efficiency and is very suitable for a model with a frequently changed topological structure. In the field of computer graphics, MSM-based soft tissue modeling has been widely used. Currently, continuous collision detection algorithms have been used by researchers in virtual surgery.

However, the conventional MSM still has the defects (a support ball spring model for soft tissue deformation simulation, computer application and software, 2013 (01): 109 and 112): 1. the accuracy of the model is closely related to the spring coefficient, and the coefficient size is usually determined only by experience; 2. the mechanical principle of the traditional MSM is based on Hooke's law, so that the biological characteristics of biological soft tissues, such as nonlinearity, volume retentivity and the like, cannot be accurately simulated. 3. The calculation result of the differential equation brings inevitable errors, and the accuracy and the stability of the model are influenced. Researchers have made many contributions to improving traditional MSM.

Qin et al propose a multi-layered MSM model Framework that simulates the biomechanics of human skin and skeletal muscle by introducing a dual-spring system (A Novel Modeling frame for multi-layered Soft Tissue Deformation in visual organic Surgery, Journal of medical Systems, 2010, 34 (3): 261-271). Duan et al propose an MSM model for Liver Gallbladder that uses implicit integration and imposes constraints on the length of the spring to maintain the volume invariance of the model (Synchronous Simulation for Deformation of Liver and Galllayder with Stretch and Compression, Conference recommendations:. annular International Conference of the IEEE Engineering in Medicine and Biology society. Conference, 2013, 2013 (2013): 4941) 4944. Li et al propose that the bending Spring improves the accuracy of the traditional MSM under large deformations (A Surface Mass-Spring Model With New Deformation Springs and precision Detection Algorithms Based on Volume Structure for read-Time Soft-Tissue Deformation Interaction, IEEE Access, 2018, 6: 75572-. However, none of their solutions completely addresses the drawbacks of the conventional models described above.

In summary, due to the complexity of the soft tissue structure of the human body, the modeling method of the soft tissue model is always a key and difficult part of the virtual surgery system, and there are few soft tissue modeling methods that can meet the requirements of the virtual surgery reality.

Disclosure of Invention

Aiming at the defects and problems in the prior art, the invention aims to provide a soft tissue modeling method based on position constraint and a nonlinear spring, which is suitable for tissue organs with various shapes and sizes and can enable a model to better meet the practical requirements of virtual surgery.

The invention is realized by the following technical scheme:

a method of soft tissue modeling based on position constraints and non-linear springs, the method comprising the steps of:

step 1, constructing a soft tissue model for acquired medical image information by using a three-dimensional modeling technology, wherein the model is composed of tetrahedral units formed by nonlinear springs and virtual body springs, and the tetrahedral model structure is used for acquiring surface node information and model internal volume information;

step 2, simulating deformation behavior of the soft tissue model, wherein in the deformation process of the model, the nonlinear spring enables the stress and the strain of the model to present linear and nonlinear relations, the virtual body spring can inhibit the change of the volume of the model, and the virtual body spring can act together to simulate the biomechanical characteristics of the soft tissue;

step 3, estimating the approximate position of each node through a kinetic equation, and then judging whether the state of a node constraint equation changes or not, wherein the node constraint comprises spring length constraint and spring curvature constraint, the two new constraints are introduced into the positions of the nodes to increase the stability of the model, and the new constraints not only compensate the limit distance between the nodes, but also correct the relative positions of the nodes according to the bending angles of the springs;

if the initial state of the equation changes, correcting the node to a correct position within a reasonable range; if the initial state of the equation is not changed, the nonlinear spring and the virtual body spring simulate the deformation effect of soft tissues;

and 4, finishing the single-step iteration of deformation calculation and entering the next round of circulation.

Further, in the deformation calculation process in the step 3, the resultant force received by the node is calculated according to Newton kinematics, the node position is preliminarily estimated, whether the node position meets the constraint equation of the length and the bending angle of the spring or not is judged according to the preliminarily calculated node position, and if the node position meets the constraint equation, the nonlinear spring and the virtual body spring simulate the deformation effect of the soft tissue; otherwise the position and velocity of the node will be modified to satisfy the constraint equations, which gives the model better controllability.

The specific steps of step 3 are:

3.1 for any node i in the model, estimating the approximate position of p according to a mechanical differential equation:

here, M_iRepresenting the quality of the node i in terms of,represents the acceleration of node i, F (u)_i) Representing the sum of internal forces, F, experienced by node i_{ext_i}Indicating its external force. F (u)_i) The spring consists of a spring force and a spring damping force:

the expression above indicates that the node i is subjected to an internal force equal to the sum of the spring force and the damping force, and the volume force of all the nodes adjacent to the node i. Wherein k is_ijRepresenting the elastic coefficient between node i and node j, c_ijRepresenting the spring damping coefficient between node i and node j,representing the first derivative to the displacement, i.e. the velocity of node i. l_ij＝u_i-u_jRepresenting the distance vector of the spring between node i and node j,representing the initial distance of the spring between node i and node j.

In order to achieve a non-linear deformation behavior of biological soft tissue, the expression of the non-linear spring force and its corresponding coefficient should be as follows:

in the formula, k_ijThe spring elastic coefficient connecting the node i and the node j; k is a radical of₁And k₂Is a constant. Δ l_ijIndicating the change in length of the spring. It can be seen that the relationship between spring force and length change is a cubic polynomial at smaller displacements and linear at larger displacements.

The formula for the volume force is as follows:

here V_jRepresenting the volume, k, of tetrahedral units containing a node i₃Is the virtual body spring rate.Andrepresenting the current volume and the initial volume of the tetrahedral unit, respectively.

The damping force of the spring can be used to simulate the viscoelasticity of soft tissue. To achieve this viscoelastic characteristic, the formula for the damping force is as follows:

in order to realize the volume retentivity of the soft tissue, not only the stress between the surface nodes of the soft tissue model but also the influence of the internal structure on the surface nodes are considered.

In the above, the system of differential equations for node i can be derived:

thus, the positions of the available nodes i are:

u_t+Δt＝u_t+v_i·Δt

3.2 judging whether the position of the node i meets the constraint equation or not, and if the position of the node meets the constraint equation, reserving the estimated node position in 3.1; and if not, correcting the node position according to the constraint equation. In particular, the invention proposes a length constraint and a curvature constraint of the spring.

The spring length constraint applies to the distance between adjacent nodes. If the length deformation rate of the spring between the nodes is larger than the preset deformation rate, the positions of the nodes are directly corrected to meet the critical deformation rate. The rate of deformation of the spring length is as follows:

τ＝|l_cur-l_init|/l_init

here, τ represents the deformation rate l of the spring_initDenotes the initial length of the spring,/_curIndicating the current length of the spring. The over-deformation rate can judge the compression or the extension of the springAnd then using a constraint equation to generate a position correction for the node. For each node, the spring length constraint equation is defined as:

if the spring length exceeds the threshold value and the state of the constraint equation changes, the node p at the two ends of the spring_iAnd p_jMaking corresponding corrections:

similar to the spring length constraint, define δ as the spring bending, as follows:

wherein, p'_iAnd p'_jIs the initial position of the node at both ends of the spring, p_iAnd p_jThe current positions of the nodes at the two ends are respectively. If delta₀Representing the critical value of the spring bending, the spring bending constraint equation is:

if the bending degree of the spring deformation is larger than the critical value, the nodes p at the two ends of the spring_iAnd p_jThe position is revised to satisfy the constraint equation:

compared with the prior art, the invention adopts a new nonlinear spring force system and a node position constraint mechanism, so that the model can effectively feed back the biomechanical characteristics of soft tissues. The new spring constraint can limit the threshold value of the spring deformation strength, such as the ultimate extension length or bending angle. The system after the position constraint is unconditionally stable, and the defect that the traditional MSM model is unstable is overcome.

Drawings

FIG. 1 is a flow chart of the MSM model deformation calculation of the present invention.

FIG. 2 is a schematic diagram of the basic tetrahedral unit structure forming the soft tissue model of the present invention.

Fig. 3 is a schematic view of the virtual body spring in the tetrahedral unit of the present invention.

FIG. 4 is a schematic length constraint diagram of the spring of the present invention.

FIG. 5 is a schematic view of the curvature constraint of the spring of the present invention.

Detailed Description

The invention will be further described with reference to the accompanying drawings.

A soft tissue modeling method based on position constraint and nonlinear springs is shown in a flow chart of FIG. 1, and comprises the following steps:

s1, constructing a soft tissue model for acquired medical image information by using a three-dimensional modeling technology, wherein the model is composed of tetrahedral units formed by nonlinear springs and virtual body springs as shown in a figure 2 and a figure 3;

s2, simulating deformation behavior of the soft tissue model, wherein in the deformation process of the model, the nonlinear spring enables the stress and the strain of the model to present linear and nonlinear relations, the virtual body spring can inhibit the change of the volume of the model, and the virtual body spring can act together to simulate the biomechanical characteristics of the soft tissue;

s3, estimating the approximate position of each node through a kinetic equation, and then judging whether the state of a node constraint equation changes or not, wherein the node constraint comprises spring length constraint shown in figure 4 and spring bending constraint shown in figure 5; if the initial state of the equation changes, correcting the node to a correct position within a reasonable range; if the initial state of the equation is not changed, the nonlinear spring and the virtual body spring simulate the deformation effect of soft tissues;

3.1 for any node i in the model, estimating the approximate position of p according to a mechanical differential equation:

here, M_iRepresenting the quality of the node i in terms of,represents the acceleration of node i, F (u)_i) Representing the sum of internal forces, F, experienced by node i_{ext_i}An external force indicative thereof;

F(u_i) The spring consists of a spring force and a spring damping force:

the expression above indicates that the node i is subjected to an internal force equal to the sum of the spring force and the damping force, and the volume force of all the nodes adjacent to the node i. Wherein k is_ijRepresenting the elastic coefficient between node i and node j, c_ijRepresenting the spring damping coefficient between node i and node j,representing the first derivative to the displacement, i.e. the velocity of node i. l_ij＝u_i-u_jRepresenting the distance vector of the spring between node i and node j,representing the initial distance of the spring between node i and node j.

In order to achieve a non-linear deformation behavior of biological soft tissue, the expression of the non-linear spring force and its corresponding coefficient should be as follows:

in the formula, k_ijThe spring elastic coefficient connecting the node i and the node j; k is a radical of₁And k₂Is a constant. Δ l_ijIndicating a change in length of the spring; it can be seen that the relationship between spring force and length change is a cubic polynomial at smaller displacements and linear at larger displacements;

the formula for the volume force is as follows:

here V_jRepresenting the volume, k, of tetrahedral units containing a node i₃Is the virtual body spring rate.Andrespectively representing the current volume and the initial volume of the tetrahedral unit;

the damping force of the spring can be used to simulate the viscoelasticity of soft tissue, and to achieve this viscoelasticity characteristic, the formula for the damping force is as follows:

in order to realize the volume retentivity of the soft tissue, not only the stress between the surface nodes of the soft tissue model but also the influence of the internal structure on the surface nodes are considered.

In the above, the system of differential equations for node i can be derived:

thus, the positions of the available nodes i are: u. of_t+Δt＝u_t+v_i·Δt

3.2 judging whether the position of the node i meets the constraint equation or not, and if the position of the node meets the constraint equation, reserving the estimated node position in 3.1; if not, the node position is corrected according to the constraint equation.

The spring length constraint applies to the distance between adjacent nodes. If the length deformation rate of the spring between the nodes is larger than the preset deformation rate, the positions of the nodes are directly corrected to meet the critical deformation rate. The rate of deformation of the spring length is as follows: τ ═ l_cur-l_init|/l_init

Here, τ represents the deformation rate of the spring, l_initDenotes the initial length of the spring,/_curIndicating the current length of the spring. The over-deformation rate can determine the degree of compression or extension of the spring, which in turn uses a constraint equation to generate a position correction for the node. For each node, the spring length constraint equation is defined as:

if the spring length exceeds the threshold value and the state of the constraint equation changes, the node p at the two ends of the spring_iAnd p_jMaking corresponding corrections:

similar to the spring length constraint, define δ as the spring bending, as follows:

wherein, p'_iAnd p'_jIs the initial position of the node at both ends of the spring, p_iAnd p_jThe current positions of the nodes at the two ends are respectively. If delta₀Representing the critical value of the spring bending, the spring bending constraint equation is:

if the bending degree of the spring deformation is larger than the critical value, the nodes p at the two ends of the spring_iAnd p_jThe position is revised to satisfy the constraint equation:

and S4, finishing single-step iteration of deformation calculation and entering the next round of circulation.

The MSM soft tissue model provided by the embodiment is composed of 2032 tetrahedral units, 722 triangular surfaces and 515 nodes. The mass of each node in the model is 1g, and the parameter k of the nonlinear spring₁Is 4, k₂Is 2, Δ l_c1.5mm, A3.5 and B4. Damping coefficient b₀Is 0.5, b₁0.1, parameter k of the virtual body spring₃Is 1. Spring length constrained critical deformation rate tau₀20% of the spring bending constraint bending threshold delta₀Is 30 degrees; for the adjacent node p₁、p₂If p is₁The volume change of the tetrahedral unit is 5mm³The coordinates of the center of the tetrahedron are (1, 0, 0); p is a radical of₁At a speed of 4mm/s and a direction ofThe initial coordinate is (2, 0, 0), and the current coordinate is (3, 0, 0); p is a radical of₂The speed is 0mm/s, the initial coordinate and the current coordinate are both (0, 0, 0), the time step length is 0.1s, and p is calculated₁Position at time t. In practical applications, the calculation process for all nodes in the model is the same as in this example.

First, p is estimated according to a mechanical differential equation₁Approximate position of (2):

(1) the deformation of the spring is 1mm and less than delta l_cThen node p₁The non-linear spring force of (2) is 4 x 1+2 x 1³6N in the direction of

(2)p₁The damping force is (0.5+4 x 0.1) x 1 ═ 0.9N, and the direction is

(3)p₁The spring force of the virtual body is 1 × 5 ═ 5N, and the direction is

(4)p₁The speed change of (d) is:

(5)p₁the new coordinate change of (2) is:

second, judge the node p₁Whether the position of (d) satisfies the constraint equation:

(1) judging node p₁Whether the position of (a) satisfies the spring length constraint. The length variation of the spring is 21%, if the length constraint of the spring is not satisfied, the node p is₁Correction of (Δ p)₁Is 0.01mm and has a direction ofNode p₂Correction of (Δ p)₂Is 0.01mm and has a direction of

(2) Judging node p₁Whether the position of (a) satisfies the spring bending constraint. The bending degree of the spring is 0 degrees and does not reach a critical value, so that a bending constraint equation is satisfied.

Finally, the node p is calculated₁And node p₂The final positions of (2.20, 0, 0), (0.01, 0, 0).

The foregoing merely represents preferred embodiments of the invention, which are described in some detail and detail, and therefore should not be construed as limiting the scope of the invention. It should be noted that, for those skilled in the art, various changes, modifications and substitutions can be made without departing from the spirit of the present invention, and these are all within the scope of the present invention. Therefore, the protection scope of the present patent shall be subject to the appended claims.