CN113409443B - Soft tissue modeling method based on position constraint and nonlinear spring - Google Patents

Abstract

The invention discloses a soft tissue modeling method based on position constraints and nonlinear springs. (3) Estimate the approximate position of each node, and then judge whether the state of the node constraint equation has changed, including the spring length and curvature constraints; if the initial state of the equation changes, the node will be corrected to the correct position within a reasonable range ; If the initial state of the equation does not change, the nonlinear spring and virtual body spring will simulate the deformation effect of the soft tissue; (4) The single-step iteration of calculating the deformation ends, and the next cycle is entered. The present invention adopts a nonlinear spring system and a virtual body spring to simulate the biomechanical properties of soft tissues, and also adopts a new constraint equation to limit the node displacement to enhance the stability of the model, and improves the traditional mass spring model with poor accuracy and instability. defect.

Description

A Soft Tissue Modeling Method Based on Position Constraints and Nonlinear Springs

technical field

The invention belongs to the technical field of soft tissue deformation model modeling, in particular to a soft tissue modeling method based on position constraints and nonlinear springs, which is used to solve the soft tissue modeling problem in virtual surgery experiments.

Background technique

Virtual surgery is an important surgical assistant training system that can realistically simulate the surgical environment in a low-cost mode. Virtual surgery can not only provide users with visualization needs, but also make doctors familiar with the use of surgical instruments and accumulate surgical skills. At present, a perfect virtual surgery system should have three requirements at the same time: real-time, accuracy and 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.

Soft tissue modeling methods mainly include geometric modeling and physical modeling. The geometric modeling method redraws the geometric model of the soft tissue by 3D reconstruction of the medical image, while the physical modeling method establishes the physical model according to the biomechanical and kinematic properties of the soft tissue. The geometric modeling method is simple in calculation, but does not consider the constitutive relationship of the object material, so the accuracy of the model cannot be guaranteed (Research on Soft Tissue Deformation and Tear Model in Virtual Brain Surgery, Nanchang University, 2016). At present, 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 & Bioengineering. IEEE , 2009: 22-24).

FEM is a modeling method based on continuum theory that discretizes continuous domain problems into several finite elements. The deformation of the model is obtained based on biomechanics. The soft tissue model based on FEM has high accuracy (An Introduction to The Finite Element Method, Topics in EngineeringMathematics, 1992(81):37-60). Therefore, the finite element method has been widely used in soft tissue modeling. Although FEM can accurately simulate the mechanical properties, the biological properties of human soft tissues are usually complex and the internal structure is fine, which will cause the FEM solution to consume a lot of resources and cannot meet the real-time interaction requirements of virtual surgery (A SurfaceMass-Spring Model With New Flexion Springs and Collision Detection AlgorithmsBased on Volume Structure for Real-Time Soft-Tissue Deformation Interaction, IEEE Access, 2018, 6:75572-75597).

Different from the continuous physical model, the MSM discretizes the soft tissue model into multiple particles, and the interaction between the particles is realized by the springs connecting the particles (Calibration of Mass Spring Models for Organ Simulations, IEEE/RSJ International Conference on Intelligent Robots&Systems. IEEE, 2007 : 370-375). Compared with FEM, MSM has high computational efficiency and is very suitable for models with frequently changing topology. In the field of computer graphics, MSM-based soft tissue modeling has been widely used. At present, some researchers have used the continuous collision detection algorithm in virtual surgery.

However, the traditional MSM still has defects (a supported ball spring model for soft tissue deformation simulation, Computer Application and Software, 2013(01): 109-112): 1. The accuracy of the model is closely related to the spring coefficient, and the coefficient The size is usually determined only by experience; 2. The mechanical principle of traditional MSM is based on Hooke's law, so it cannot accurately simulate the biological properties of biological soft tissues, such as nonlinearity, volume retention, etc. 3. The calculation results of differential equations bring inevitable errors, which affect the accuracy and stability of the model. Researchers have made many contributions to improve traditional MSM.

Qin et al. proposed a multi-layer MSM model framework to simulate the biomechanical properties of human skin and skeletal muscle by introducing a double-spring system (A Novel Modeling Framework for Multilayered Soft Tissue Deformation in Virtual Orthopedic Surgery, Journal of Medical Systems, 2010, 34 (3): 261-271). Duan et al. proposed an MSM model suitable for the liver gallbladder. They employed implicit integration and imposed constraints on the length of the spring to keep the volume invariant of the model (Synchronous Simulation for Deformation of Liver and Gallbladder with Stretch and Compression Compensation, Conference procedures: . Annual International Conference of the IEEE Engineering in Medicine and Biology Society. IEEE Engineering in Medicine and Biology Society. Conference, 2013, 2013(2013): 4941-4944). Li et al. proposed that the bending spring improves the accuracy of traditional MSM under large deformation (A Surface Mass-Spring Model With New Flexion Springs and Collision Detection Algorithms Based on Volume Structure for Real-TimeSoft-Tissue Deformation Interaction, IEEE Access, 2018, 6:75572-75597). However, none of their schemes fully address the above-mentioned shortcomings of traditional models.

To sum up, due to the complexity of human soft tissue structure, the modeling method of soft tissue model has always been a key and difficult part of virtual surgery system, and there are few soft tissue modeling methods that can meet the needs of virtual surgery reality.

SUMMARY OF THE INVENTION

In view of the deficiencies and difficulties in the prior art, the present invention aims to provide a soft tissue modeling method based on position constraints and nonlinear springs. Meet the real needs of virtual surgery.

The present invention is achieved through the following technical solutions:

A soft tissue modeling method based on position constraints and nonlinear springs, the method includes the following steps:

Step 1. Use 3D modeling technology to build a soft tissue model for the collected medical image information, the model is composed of tetrahedral elements formed by nonlinear springs and virtual body springs, and the tetrahedral model structure is used to obtain surface node information and the internal volume of the model information;

Step 2. Simulate the deformation behavior of the soft tissue model. During the deformation process of the model, the nonlinear spring makes the stress and strain of the model show a linear and nonlinear relationship. The virtual body spring can suppress the change of the model volume, and they will work together to simulate Biomechanical properties of soft tissue;

Step 3. Estimate the approximate position of each node through the dynamic equation, and then judge whether the state of the node constraint equation has changed. The node constraint includes the spring length constraint and the spring bending degree constraint. The above two new constraints are introduced at the node position. , to increase the stability of the model, the new constraint not only compensates the limit distance between nodes, but also corrects the relative position of the nodes according to the bending angle of the spring;

If the initial state of the equation changes, the node will be corrected to the correct position within a reasonable range; if the initial state of the equation does not change, the nonlinear spring and virtual body spring will simulate the deformation effect of soft tissue;

Step 4. The single-step iteration of calculating the deformation ends, and the next cycle is entered.

Further, in the deformation calculation process of step 3, according to Newtonian kinematics, calculate the resultant force received by the node and preliminarily estimate the node position, and then judge whether it satisfies the constraint equation of the spring length and bending angle according to the preliminarily calculated node position, if If satisfied, the nonlinear spring and virtual body spring will simulate the deformation effect of soft tissue; otherwise, the positions and velocities of the nodes will be modified to satisfy the constraint equations, which makes the model more controllable.

The specific steps of step 3 are:

3.1 For any node i in the model, estimate the approximate position of p according to the mechanical differential equation:

表示节点i的加速度，F(u_i)表示节点i所受内力和，F_{ext_i}表示它的外力。F(u_i)由弹簧力和弹簧阻尼力构成：Here, M _i represents the mass of node i,

represents the acceleration of node i, F(u _i ) represents the sum of internal forces on node i, and F _{ext_i} represents its external force. F(u _i ) consists of spring force and spring damping force:

表示对位移的一阶导数即节点i的速度。l_ij＝u_i-u_j表示节点i和节点j之间弹簧的距离矢量，

表示节点i和节点j之间弹簧的初始距离。The above formula indicates that the internal force of node i is equal to the sum of the spring force, damping force and body force of all its adjacent nodes. where k _ij represents the elastic coefficient between node i and node j, c _ij represents the spring damping coefficient between node i and node j,

Represents the first derivative of the displacement, the velocity of node i. l _ij =u _i -u _j represents the distance vector of the spring between node i and node j,

represents the initial distance of the spring between node i and node j.

In order to realize the nonlinear deformation behavior of biological soft tissue, the expressions of nonlinear spring force and its corresponding coefficient should be as follows:

In the formula, k _ij is the spring elastic coefficient connecting node i and node j; k ₁ and k ₂ are constants. Δl _ij represents the change in spring length. It can be seen that 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, it is a linear relationship.

The formula for body force is as follows:

和

分别表示四面体单元的当前体积和初始体积。Here _Vj represents the volume of the tetrahedral element containing node i, and _k3 is the virtual volume spring coefficient.

and

represent the current volume and initial volume of the tetrahedral element, respectively.

The damping force of the spring can be used to simulate the viscoelasticity of soft tissue. To achieve this viscoelastic feature, the damping force is formulated as:

In order to achieve the volume retention of soft tissue, not only the force between the surface nodes of the soft tissue model, but also the influence of the internal structure on the surface nodes should be considered.

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

In this way, the position of node i can be obtained as:

u _t+Δt = _{u t} +vi ·Δt

3.2 Determine whether the position of node i satisfies the constraint equation, if the node position satisfies the constraint equation, keep the node position estimated in 3.1; if not, correct the node position according to the constraint equation. Specifically, the present invention proposes a length constraint and a curvature constraint for the spring.

Spring length constraints apply to the distance between adjacent nodes. If the deformation rate of the spring length between the nodes is greater than the preset deformation rate, the node position is directly modified 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 _init , the initial length of the spring, and l _cur , the current length of the spring. The rate of deformation can determine how much the spring is compressed or stretched, and then use the constraint equations to generate position corrections for the nodes. For each node, define the spring length constraint equation as:

If the length of the spring exceeds the threshold and the state of the constraint equation changes, the nodes p _i and p _j at both ends of the spring are corrected accordingly:

Similar to the spring length constraint, δ is defined as the spring curvature, as follows:

Among them, p′ _i and p′ _j are the initial positions of the nodes at both ends of the spring, and p _i and p _j are the current positions of the nodes at both ends respectively. If δ ₀ represents the critical value of spring bending, the spring bending constraint equation is:

If the curvature of the spring deformation is greater than the critical value, the positions of the nodes p _i and p _j at both ends of the spring are re-corrected to satisfy the constraint equation:

Compared with the prior art, the present invention adopts a new nonlinear spring force system and a mechanism of node position constraint, so that the model can effectively feed back the biomechanical properties of the soft tissue. Instead, it prejudges whether the node positions satisfy the proposed constraint equations—including the new spring length constraint and spring bending constraint. New spring constraints can limit the critical value of spring deformation strength, such as the ultimate stretch length or bending angle. The system after the position constraint is unconditionally stable, which makes up for the instability of the traditional MSM model.

Description of drawings

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

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

FIG. 3 is a schematic diagram of a virtual body spring in a tetrahedral unit of the present invention.

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

FIG. 5 is a schematic diagram of the bending degree constraint of the spring of the present invention.

Detailed ways

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

The soft tissue modeling method based on position constraints and nonlinear springs, as shown in the flowchart of Figure 1, the method steps include:

S1. Use three-dimensional modeling technology to construct a soft tissue model for the collected medical image information, and the model is composed of tetrahedral elements formed by nonlinear springs and virtual body springs as shown in Figure 2 and Figure 3;

S2. Simulate the deformation behavior of the soft tissue model. During the deformation process of the model, the nonlinear spring makes the stress and strain of the model show a linear and nonlinear relationship. The virtual body spring can suppress the change of the model volume, and they will work together to simulate the biological behavior of the soft tissue. mechanical properties;

S3. Estimate the approximate position of each node through the dynamic equation, and then judge whether the state of the node constraint equation has changed. The node constraint includes the spring length constraint as shown in Figure 4 and the spring bending degree constraint as shown in Figure 5; if If the initial state of the equation changes, the node will be corrected to the correct position within a reasonable range; if the initial state of the equation does not change, the nonlinear spring and virtual body spring will simulate the deformation effect of soft tissue;

3.1 For any node i in the model, estimate the approximate position of p according to the mechanical differential equation:

表示节点i的加速度，F(u_i)表示节点i所受内力和，F_{ext_i}表示它的外力；Here, M _i represents the mass of node i,

represents the acceleration of node i, F(u _i ) represents the sum of the internal forces on node i, and F _{ext_i} represents its external force;

F(u _i ) consists of spring force and spring damping force:

表示对位移的一阶导数即节点i的速度。l_ij＝u_i-u_j表示节点i和节点j之间弹簧的距离矢量，

表示节点i和节点j之间弹簧的初始距离。The above formula indicates that the internal force of node i is equal to the sum of the spring force, damping force and body force of all its adjacent nodes. where k _ij represents the elastic coefficient between node i and node j, c _ij represents the spring damping coefficient between node i and node j,

Represents the first derivative of the displacement, the velocity of node i. l _ij =u _i -u _j represents the distance vector of the spring between node i and node j,

represents the initial distance of the spring between node i and node j.

In order to realize the nonlinear deformation behavior of biological soft tissue, the expressions of nonlinear spring force and its corresponding coefficient should be as follows:

In the formula, k _ij is the spring elastic coefficient connecting node i and node j; k ₁ and k ₂ are constants. Δl _ij represents the change of spring length; it can be seen that 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, it is a linear relationship;

The formula for body force is as follows:

和

分别表示四面体单元的当前体积和初始体积；Here _Vj represents the volume of the tetrahedral element containing node i, and _k3 is the virtual volume spring coefficient.

and

represent the current volume and initial volume of the tetrahedral element, respectively;

The damping force of the spring can be used to simulate the viscoelasticity of soft tissue. To achieve this viscoelasticity, the damping force is formulated as:

In order to achieve the volume retention of soft tissue, not only the force between the surface nodes of the soft tissue model, but also the influence of the internal structure on the surface nodes should be considered.

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

In this way, the position of node i can be obtained as: u _t+Δt =u _t +v _i ·Δt

3.2 Determine whether the position of node i satisfies the constraint equation. If the node position satisfies the constraint equation, keep the node position estimated in 3.1; if not, correct the node position according to the constraint equation.

Spring length constraints apply to the distance between adjacent nodes. If the deformation rate of the spring length between the nodes is greater than the preset deformation rate, the node position is directly modified to meet the critical deformation rate. The deformation rate of the spring length is as follows: τ=|l _cur -l _init |/l _init

Here, τ is the deformation rate of the spring, l _init is the initial length of the spring, and l _cur is the current length of the spring. The rate of deformation can determine how much the spring is compressed or stretched, and then use the constraint equations to generate position corrections for the nodes. For each node, define the spring length constraint equation as:

If the length of the spring exceeds the threshold and the state of the constraint equation changes, the nodes p _i and p _j at both ends of the spring are corrected accordingly:

Similar to the spring length constraint, δ is defined as the spring curvature, as follows:

Among them, p′ _i and p′ _j are the initial positions of the nodes at both ends of the spring, and p _i and p _j are the current positions of the nodes at both ends respectively. If δ ₀ represents the critical value of spring bending, the spring bending constraint equation is:

If the curvature of the spring deformation is greater than the critical value, the positions of the nodes p _i and p _j at both ends of the spring are re-corrected to satisfy the constraint equation:

S4. The single-step iteration of calculating the deformation ends, and the next cycle is entered.

初始坐标为(2，0，0)，当前坐标为(3，0，0)；p₂速度为0mm/s，初始坐标和当前坐标均为(0，0，0)，时间步长为0.1s，计算p₁在t时刻的位置。在实际应用中，模型内所有节点的计算过程均与本例中相同。The MSM soft tissue model proposed in this embodiment consists of 2032 tetrahedral elements, 722 triangular faces, and 515 nodes. The masses of the nodes in the model are all 1g, the parameters _k1 of the nonlinear spring are 4, k2 is ₂ , _Δlc is 1.5mm, A is 3.5, and B is 4. The damping force coefficient b ₀ is 0.5, b ₁ is 0.1, and the parameter k ₃ of the virtual body spring is 1. The critical deformation rate τ ₀ of the spring length constraint is 20%, and the bending critical value δ ₀ of the spring bending constraint is 30°; for the adjacent nodes p ₁ and p ₂ , if the volume of the tetrahedral element where p ₁ is located is 5mm ³ , The coordinates of its tetrahedron center are (1, 0, 0); the speed of p ₁ is 4mm/s, and the direction is

The initial coordinates are (2, 0, 0), the current coordinates are (3, 0, 0); the speed of p ₂ is 0mm/s, the initial and current coordinates are (0, 0, 0), and the time step is 0.1 s, calculate the position of p ₁ at time t. In practical applications, the calculation process of all nodes in the model is the same as in this example.

The first step is to estimate the approximate position of p ₁ according to the mechanical differential equation:

(1) The deformation of the spring is 1mm, which is less than Δl _c , then the nonlinear spring force of the node p ₁ is 4*1+2*1 ³ =6N, and the direction is

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

(3) The virtual body spring force on p ₁ is 1*5=5N, and the direction is

(4) The speed change of p ₁ is:

(5) The new coordinate change of p ₁ is:

The second step is to determine whether the position of node p ₁ satisfies the constraint equation:

节点p₂的修正Δp₂为0.01mm，方向为

(1) Determine whether the position of the node p ₁ satisfies the spring length constraint. The length change of the spring is 21%, and the length constraint of the spring is not satisfied, then the correction Δp ₁ of the node p ₁ is 0.01mm, and the direction is

The correction Δp ₂ for node p ₂ is 0.01mm and the direction is

(2) Determine whether the position of node p ₁ satisfies the spring bending constraint. The bending of the spring is 0°, which does not reach the critical value, so the bending constraint equation is satisfied.

Finally, the final positions of node p ₁ and node p ₂ are calculated as (2.20, 0, 0), (0.01, 0, 0).

The above description only expresses the preferred embodiments of the present invention, and the description thereof is relatively specific and detailed, but should not be construed as a limitation on the scope of the patent of the present invention. It should be pointed out that for those of ordinary skill in the art, without departing from the concept of the present invention, several modifications, improvements and substitutions can be made, which all belong to the protection scope of the present invention. Therefore, the protection scope of the patent of the present invention should be subject to the appended claims.

Claims (3)

1. a soft tissue modeling method based on position constraint and nonlinear spring, is characterized in that, described method comprises the following steps: Step 1: Use 3D modeling technology to build a soft tissue model for the collected medical image information, the model is composed of tetrahedral elements formed by nonlinear springs and virtual body springs, and the tetrahedral model structure is used to obtain surface node information and the internal volume of the model information; Step 2: Simulate the deformation behavior of the soft tissue model. During the deformation process of the model, the nonlinear spring makes the stress and strain of the model show a linear and nonlinear relationship. The virtual body spring can suppress the change of the model volume, and they will work together to simulate Biomechanical properties of soft tissue; Step 3: Estimate the approximate position of each node through the dynamic equation, and then judge whether the state of the node constraint equation has changed. The node constraint includes the spring length constraint and the spring bending degree constraint; if the initial state of the equation changes, correct the node to a reasonable value. correct position within the range; nonlinear springs and virtual body springs will simulate the deformation effects of soft tissue if the initial state of the equation is unchanged; specifically: 3.1 For any node i in the model, calculate the resultant force on the node and preliminarily estimate the node position; The formula for calculating the resultant force on node i is:

where M _i represents the quality of node i,

represents the acceleration of node i, F(u _i ) represents the sum of the internal forces on node i, and F _{ext_i} represents its external force; The position expression of node i is: u _t+Δt = u _t +v _i ·Δt, u _t+Δt and u _t represent the displacement vectors at time t+Δt and time t respectively, and v _i represents the speed of node i; The approximate location estimation derivation process of each node is as follows:

The above formula indicates that the internal force of node i is equal to the sum of spring force, damping force and body force of all its adjacent nodes; among them, k _ij represents the elastic coefficient between node i and node j, and c _ij represents node i and node the spring damping coefficient between j,

represents the first derivative of the displacement, i.e. the velocity of node i, l _ij =u _i -u _j represents the distance vector of the spring between node i and node j,

represents 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 _ij is the spring elastic coefficient connecting node i and node j; k ₁ and k ₂ are constants; l _ij represents the change of spring length, ui is the position of node _{i, ||u i -u j} || is the distance between node i and 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, it is a linear relationship; The formula for body force is as follows:

In the formula, V _j represents the volume of the tetrahedral element including node i, k ₃ is the spring coefficient of the virtual body; V _j ^cur and V _j ^init represent the current volume and initial volume of the tetrahedral element, respectively, X _i represents the coordinate of node i , X _j represents the coordinates of node j; The damping force of the spring is used to simulate the viscoelasticity of soft tissue. The formula of the spring damping force is as follows:

Derive the system of differential equations for node i:

In the formula, F _i represents the sum of the spring force and the body force of the node _i , vi represents the speed of the node i,

is the acceleration of node i; According to the above differential equation system, the position of node i is obtained; 3.2 Determine whether the position of node i satisfies the constraint equation, if the node position satisfies the constraint equation, keep the node position estimated in 3.1; if not, correct the node position according to the constraint equation; the constraint equation includes spring length constraint and spring bending degree constraint equation; The spring length constraint is applicable to the distance between adjacent nodes. If the spring length deformation rate between the nodes is greater than the preset deformation rate, the node position is directly modified to meet the critical deformation rate; If the length of the spring exceeds the threshold and the state of the constraint equation changes or the curvature of the spring deformation is greater than the threshold, the nodes p _i and p _j at both ends of the spring are corrected accordingly; Step 4: The single-step iteration for calculating the deformation ends, and the next cycle is entered. 2. A soft tissue modeling method based on position constraints and nonlinear springs according to claim 1, wherein the deformation rate of the spring length in the step 3.2 is as follows:

In the formula, τ represents the deformation rate of the spring length, and l _cur and l _init represent the current length and initial length of the spring, respectively; the degree of compression or stretching of the spring is judged by the deformation rate, and then the constraint equation is used to generate position corrections for the nodes; For each node, the spring length constraint equation is:

In the formula, τ ₀ represents the critical value of the spring length deformation rate; If the length of the spring exceeds the threshold and the state of the constraint equation changes, the nodes p _i and p _j at both ends of the spring are corrected accordingly:

3. a kind of soft tissue modeling method based on position constraint and nonlinear spring according to claim 1, is characterized in that, in described step 3.2, spring bending degree δ expression is as follows:

Among them, p _i ' and p _j ' are the initial positions of the nodes at both ends of the spring, and p _i and p _j are the current positions of the nodes at both ends respectively; δ ₀ represents the critical value of the spring tortuosity, then the spring tortuosity constraint equation is:

If the curvature of the spring deformation is greater than the critical value, the positions of the nodes p _i and p _j at both ends of the spring are re-corrected to satisfy the constraint equation:

Images