Method for calculating image force feedback based on discrete elastic rod physical modeling guide wire

Abstract

The invention discloses a method for calculating image force feedback by physically modeling a guide wire based on a discrete elastic rod, which comprises the following steps of: modeling the guide wire into a non-telescopic deformable linear object, representing the center line of the guide wire by adopting an arc length parameterized three-dimensional curve, and describing bending and torsion characteristics by combining an orthogonal material frame; a curvature vector and a torsion angle are calculated through a discretized particle model to simulate the internal state of the guide wire, and the motion of the guide wire is predicted based on a semi-implicit integration method through a Lagrange motion equation; then, coarse detection is carried out in combination with the preprocessed mask image of the blood vessel wall, whether the guide wire makes contact with the blood vessel wall or not is judged, and due to the existence of a preprocessing structure, the detection is faster than a common bounding box which is divided through a tree structure; then, whether the tip of the guide wire penetrates through the blood vessel wall or not is judged through the blood vessel center line and kd-tree, a force feedback calculation model is constructed by combining a spring damping system according to deformation caused to the blood vessel wall, and corresponding feedback force after collision is calculated.

Description

Technical Field

[0001] The present invention relates to the field of guidewire control technology, and in particular to a method for calculating image force feedback of a guidewire based on discrete elastic rod physical modeling. Background Art

[0002] In cardiovascular interventional surgery, guidewires serve as the core tool for guiding instruments to accurately traverse complex vascular pathways. The real-time modeling of their spatial morphology and mechanical interaction directly affects surgical safety and operational efficiency. In recent years, with the advancement of medical imaging technology, virtual surgical simulation systems based on physical modeling have gradually become a research hotspot. Among them, the discrete elastic rod model has the ability to efficiently represent slender deformable objects. Based on the Lagrange equation of motion, the system can solve the force balance equations of each particle. Traditional methods such as explicit integration have strict restrictions on the time step, are prone to unstable phenomena such as divergence or oscillation, and have poor convergence. Implicit integration has a high computational cost because each calculation requires solving a large matrix.

[0003] In virtual surgery simulation systems, other common collision detection methods, such as the hierarchical bounding box algorithm, can achieve logarithmic time complexity, but due to its use of triangle facet comparisons, it is still relatively slow. Spatial hashing methods require the selection of a suitable hash function to avoid collisions, but hash collisions can still occur. Regarding collision response, the current common method uses a penalty method to return a feedback force after a collision, but this model is overly simplistic, resulting in different response characteristics between different materials and failing to fully represent the properties of the blood vessel wall.

[0004] Therefore, it is necessary to develop a method for calculating image force feedback based on discrete elastic rod physical modeling of guidewires to solve the above problems. Summary of the Invention

[0005] The purpose of the present invention is to design a method for calculating image force feedback of a guidewire based on discrete elastic rod physical modeling in order to solve the above problems.

[0006] The present invention achieves the above-mentioned purpose through the following technical solutions:

[0007] A method for calculating image force feedback of a guidewire based on discrete elastic rod physical modeling, comprising:

[0008] S1, electromagnetic sensor obtains guidewire information;

[0009] S2, guidewire geometry modeling;

[0010] S3, guidewire physical modeling and calculation;

[0011] S4, guidewire motion prediction;

[0012] S5, coarse detection: Use three-dimensional voxels to create a mask image to determine whether the guidewire tip has collided with the vessel wall; based on the collision information, use the method of S4 to perform motion prediction;

[0013] S6, fine detection: collision detection is performed based on the vascular wall model, vascular centerline and kd-tree;

[0014] S7, judging whether the guidewire tip will pass through the blood vessel wall based on the relationship between the guidewire tip and the blood vessel centerline, if yes, proceeding to the next step, otherwise proceeding to step S9;

[0015] S8. Calculate the collision feedback force;

[0016] S9: Check whether the simulation is finished. If not, go to step S1 and continue the loop. If yes, then the simulation is finished.

[0017] Furthermore, step S1 specifically includes: obtaining the precise spatial position information of the guidewire through the electromagnetic sensor, which provides the position of the guidewire ) and the rotation quaternion ( , converted to the coordinates of the world coordinate system through the transformation matrix: ( ), and then the guidewire is modeled according to the obtained three-dimensional position coordinates.

[0018] Furthermore, step S2 specifically includes: describing the guide wire as a coordinate three-dimensional curve: ;in It is a parameterized curve describing the arc length of the guidewire centerline in three-dimensional space. is an orthogonal material frame assigned to each point on the centerline, which contains all the information needed to calculate torsion and bending; the material frame satisfies Represents the tangent direction of the center line, and the two are perpendicular to each other. , Represents the The derivative of the centerline curvature vector .

[0019] Furthermore, step S3 specifically includes: performing physical modeling on the guide wire tip, based on Kirchhoff's elastic rod theory, converting the energy Divided into three categories: , and , the total energy is:

[0020] Since the guidewire is inextensible, the total energy is obtained by excluding the guidewire stretching energy and applying a PBD (Position-Based Dynamics) auxiliary constraint to enforce the guidewire's inextensibility:

[0021] The PBD constraint uses the Guass-Seidel method, which processes each edge separately and converges after more than 5 iterations, which can achieve non-scalability;

[0022] The bending energy and torsional energy can be calculated by the following formula:

[0023]

[0024] in represents the guidewire length, represents the curvature vector of the curve, represents the first part of the curvature vector, represents the second part of the curvature vector, represents the inherent curvature of the guidewire tip, represents the bending Young's modulus, represents the guidewire radius, represents the shear elastic modulus, represents the guidewire torsion, and D is the differential operator symbol;

[0025] Perform numerical calculations and discretize the guidewire into particles for calculation. For bending force, first calculate the curvature of each discretized particle. , The calculation depends on two adjacent edges. In the case of continuity, the size is ,in Represents the angle between two sides, and the discrete numerical calculation uses the following formula:

[0026]

[0027] in, Represents the number corresponding to each particle after discretization, ; The discretization of bending energy is expressed as:

[0028]

[0029] in ;

[0030] For the calculation of discrete torsional force, calculate the torsion angle of each point , According to the principle of minimum energy, we only need to know the torsion angle of the tip , use Newton iteration method to solve the following formula:

[0031]

[0032] in for matrix, is the gradient; is calculated as follows:

[0033]

[0034] in For rotation The rotation matrix of is a symmetric positive definite 2D matrix, matrix is calculated as follows:

[0035]

[0036] in represents the torsional modulus, The calculation method is:

[0037]

[0038] Finally updated all , The twist angle is obtained:

[0039]

[0040] The scalar , represents a discrete torsion.

[0041] Furthermore, step S4 specifically includes: using the Lagrange equation of motion to calculate the bending force of each particle of the guide wire and torsional force :

[0042]

[0043] Since the time step of the system is very short, there will be no numerical instability. The semi-implicit integration method is directly used to calculate the position change at the next moment. and speed change :

[0044]

[0045] in Represents quality, Represents the current speed, Represents the current position, Represents the current force.

[0046] Furthermore, step S5 specifically includes: in the coarse detection stage, directly obtaining whether the front end of the guide wire has collided with the blood vessel wall through voxels, and if the position Value , indicating that a collision has occurred, the guidewire motion is predicted based on the established physical model, and the speed, displacement, and force status of the guidewire at the next position are updated; if no collision occurs, the next step will continue to maintain the current state and move forward.

[0047] Furthermore, step S6 specifically includes: in the fine detection stage, the next step position is already known, and it is necessary to combine the blood vessel centerline and kd-tree to detect whether it penetrates the blood vessel wall. The kd-tree construction process will select the dimension with the largest variance each time. As a splitting axis:

[0048]

[0049] in Represents the mean of the i-th dimension data. There are k dimensions and n points in total. The randomness of the data distribution is guaranteed by maximizing the variance. Then, according to the median greed, the corresponding median of the selected dimension is selected. To divide:

[0050]

[0051] It represents taking the median. The greedy median ensures the balance of data volume in the left and right subtrees, keeps the tree height as low as possible, and achieves the log-level query distance to the nearest center line position.

[0052] Furthermore, steps S7 and S8 specifically include:

[0053] Through the above integral formula and guidewire model, the next displacement information of the guidewire has been obtained. According to the relationship between the guidewire tip and the centerline of the blood vessel, it can be judged whether the guidewire tip will penetrate the blood vessel wall. If penetration occurs, the distance between the tip and the blood vessel wall can be found. ; Because different tissues have different elastic coefficients and damping coefficient , the simple spring-damper model is:

[0054]

[0055] Combining the spring damping system with the structure of the blood vessel wall, the calculation model of the feedback force after collision can be obtained as follows:

[0056]

[0057] in represents the Young's modulus of the blood vessel wall, Represents the contact area between the blood vessel wall and the guidewire in the image. is the thickness of the blood vessel wall.

[0058] The beneficial effects of the present invention are:

[0059] In the coarse detection stage, the present invention uses image masks for collision detection. In the precise detection stage, the centerline is used to determine whether the guidewire tip will penetrate the blood vessel wall. Furthermore, a new model for calculating feedback force is proposed. Compared with the traditional penalty method that only uses simple elastic force, we combine the relevant parameters of the blood vessel wall and the linear spring damping system related to the tissue elastic coefficient to achieve the simulation of the differentiated mechanical responses of different soft and hard tissues, making the speed faster and the calculation more accurate.

[0060] The present invention utilizes the semi-implicit time integration method and the Verlet method to solve the particle motion state, which has high speed, numerical stability and stronger scalability. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 This is a flowchart of the application. DETAILED DESCRIPTION

[0062] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more apparent, the technical solutions of the embodiments of the present invention will be described clearly and completely below in conjunction with the accompanying drawings of the embodiments of the present invention. It should be understood that the described embodiments are only a portion of the embodiments of the present invention, not all of them. Generally, the components of the embodiments of the present invention described and illustrated in the drawings herein may be arranged and designed in a variety of different configurations.

[0063] Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention as claimed, but merely represents selected embodiments of the present invention. All other embodiments derived by persons of ordinary skill in the art based on the embodiments of the present invention without creative effort shall fall within the scope of protection of the present invention.

[0064] It should be noted that similar reference numerals and letters denote similar items in the following drawings, and therefore, once an item is defined in one drawing, it does not need to be further defined or explained in subsequent drawings.

[0065] In the description of the present invention, it should be understood that the terms "upper", "lower", "inside", "outside", "left", "right", etc. indicate orientations or positional relationships based on the orientations or positional relationships shown in the accompanying drawings, or are the orientations or positional relationships in which the inventive product is conventionally placed when in use, or are the orientations or positional relationships conventionally understood by those skilled in the art. These are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or component referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be understood as a limitation on the present invention.

[0066] Furthermore, the terms “first”, “second”, etc. are merely used for distinguishing descriptions and should not be understood as indicating or implying relative importance.

[0067] In the description of the present invention, it should also be noted that, unless otherwise expressly specified or limited, terms such as "disposed" and "connected" should be understood in a broad sense. For example, "connected" can mean a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection or an indirect connection through an intermediate medium; it can also mean internal communication between two components. Those skilled in the art will be able to understand the specific meanings of the above terms in the present invention based on specific circumstances.

[0068] The specific embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0069] like Figure 1 As shown, a method for calculating image force feedback of a guidewire based on discrete elastic rod physical modeling includes:

[0070] S1. Obtaining guidewire information; Step S1 specifically includes: obtaining the precise spatial position information of the guidewire through an electromagnetic sensor, which provides the position of the guidewire. ) and the rotation quaternion ( , converted to the coordinates of the world coordinate system through the transformation matrix: ( ), and then the guidewire is modeled according to the obtained three-dimensional position coordinates.

[0071] S2. Guidewire geometry modeling. Step S2 specifically includes: for a deformable linear object such as a guidewire, which has a small diameter and volume, the guidewire is described as a coordinate three-dimensional curve: ;in It is a parameterized curve describing the arc length of the guidewire centerline in three-dimensional space. is an orthogonal material frame assigned to each point on the centerline, which contains all the information needed to calculate torsion and bending; the material frame satisfies Represents the tangent direction of the center line, and the two are perpendicular to each other. , Represents the The derivative of the centerline curvature vector .

[0072] S3, guidewire physical modeling and calculation; Step S3 specifically includes: in order to calculate the elastic force on the guidewire tip when it moves in the blood vessel, the guidewire tip is physically modeled, and the energy is converted to Divided into three categories: , and , the total energy is:

[0073]

[0074] Since the guidewire is inextensible, the total energy is obtained by excluding the guidewire stretching energy and applying a PBD (Position-Based Dynamics) auxiliary constraint to enforce the guidewire's inextensibility:

[0075]

[0076] The PBD constraint uses the Guass-Seidel method, which processes each edge separately and converges after more than 5 iterations, which can achieve non-scalability;

[0077] The bending energy and torsional energy can be calculated by the following formula:

[0078]

[0079] in represents the guidewire length, represents the curvature vector of the curve, represents the first part of the curvature vector, represents the second part of the curvature vector, represents the inherent curvature of the guidewire tip, represents the bending Young's modulus, represents the guidewire radius, represents the shear elastic modulus, represents the guidewire torsion, and D is the differential operator symbol;

[0080] In order to perform numerical calculations, the guide wire needs to be discretized into particles for calculation. For the bending force, the curvature of each discretized particle is first calculated. , The calculation depends on two adjacent edges. In the case of continuity, the size is ,in Represents the angle between two sides, and the discrete numerical calculation uses the following formula:

[0081]

[0082] Among them, i represents the number corresponding to each particle after discretization, ; The discretization of bending energy is expressed as:

[0083]

[0084] in ;

[0085] For the calculation of discrete torsional force, it is necessary to calculate the torsion angle of each point , According to the principle of minimum energy, we only need to know the torsion angle of the tip , we can use Newton's iteration method to solve the following formula:

[0086]

[0087] in for matrix, g is the gradient; g is calculated as follows:

[0088]

[0089] in For rotation The rotation matrix of is a symmetric positive definite 2D matrix, matrix is calculated as follows:

[0090]

[0091] in represents the torsional modulus, The calculation method is:

[0092]

[0093] Finally updated all , The twist angle is obtained:

[0094]

[0095] The scalar , represents a discrete torsion.

[0096] S4, guidewire motion prediction; Step S4 specifically includes: using the Lagrange equation of motion to calculate the bending force of each particle of the guidewire and torsional force :

[0097]

[0098] Since the time step of the system is very short, there will be no numerical instability. The semi-implicit integration method is directly used to calculate the position change at the next moment. and speed change :

[0099]

[0100] in Represents quality, Represents the current speed, Represents the current position, and f represents the current force.

[0101] S5, coarse detection: according to the vascular wall model, a mask image is created using three-dimensional space voxels to obtain whether the front end of the guide wire has collided with the vascular wall; based on the information of whether the collision occurs, the method of S4 is used to perform motion prediction; Step S5 specifically includes: in the coarse detection stage, according to the vascular wall model, directly using voxels to obtain whether the front end of the guide wire has collided with the vascular wall, if the position Value , indicating that a collision has occurred, the guidewire motion is predicted based on the established physical model, and the speed, displacement, and force status of the guidewire at the next position are updated; if no collision occurs, the next step will continue to maintain the current state and move forward.

[0102] S6, fine detection: collision detection is performed based on the vascular wall model, the vascular centerline, and the kd-tree. Step S6 specifically includes: in the fine detection stage, the next step position is already known, and it is necessary to combine the vascular centerline and the kd-tree to detect whether it penetrates the vascular wall. The kd-tree construction process will select the dimension with the largest variance each time. As a splitting axis:

[0103]

[0104] in Represents the mean of the i-th dimension data. There are k dimensions and n points in total. The randomness of the data distribution is guaranteed by maximizing the variance. Then, according to the median greed, the corresponding median of the selected dimension is selected. To divide:

[0105]

[0106] It represents taking the median. The greedy median ensures the balance of data volume in the left and right subtrees, keeps the tree height as low as possible, and achieves the log-level query distance to the nearest center line position.

[0107] S7, judging whether the guidewire tip will pass through the blood vessel wall based on the relationship between the guidewire tip and the blood vessel centerline, if yes, proceeding to the next step, otherwise proceeding to step S9;

[0108] S8, calculating the collision feedback force; Steps S7 and S8 specifically include:

[0109] Through the above integral formula and guidewire model, the next displacement information of the guidewire has been obtained. According to the relationship between the guidewire tip and the centerline of the blood vessel, it can be judged whether the guidewire tip will penetrate the blood vessel wall. If penetration occurs, the distance between the tip and the blood vessel wall can be found. ; Because different tissues have different elastic coefficients and damping coefficient , the simple spring-damper model is:

[0110]

[0111] Combining the spring damping system with the structure of the blood vessel wall, the calculation model of the feedback force after collision can be obtained as follows:

[0112]

[0113] in represents the Young's modulus of the blood vessel wall, Represents the contact area between the blood vessel wall and the guidewire in the image. is the thickness of the blood vessel wall.

[0114] S9: Check whether the simulation is finished. If not, go to step S1 and continue the loop. If yes, then the simulation is finished.

[0115] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the technical principles of the present invention. These improvements and modifications should also be regarded as within the scope of protection of the present invention.

Claims

1. A method for calculating image force feedback of a guidewire based on discrete elastic rod physical modeling, characterized in that: include: S1, electromagnetic sensor obtains guidewire information; S2, guidewire geometry modeling; S3, guidewire physical modeling and calculation; S4, guidewire motion prediction; S5, coarse detection: Use three-dimensional voxels to create a mask image to determine whether the guidewire tip has collided with the vessel wall; based on the collision information, use the method of S4 to perform motion prediction; S6, fine detection: collision detection is performed based on the vascular wall model, vascular centerline and kd-tree; S7, judging whether the guidewire tip will pass through the blood vessel wall based on the relationship between the guidewire tip and the blood vessel centerline, if yes, proceeding to the next step, otherwise proceeding to step S9; S8. Calculate the collision feedback force; S9: Check whether the simulation is finished. If not, go to step S1 and continue the loop. If yes, then the simulation is finished.

2. The method for calculating image force feedback of a guidewire based on discrete elastic rod physical modeling according to claim 1, characterized in that: Step S1 specifically includes: obtaining the precise spatial position information of the guidewire through an electromagnetic sensor, which provides the position of the guidewire. ) and the rotation quaternion ( , converted to the coordinates of the world coordinate system through the transformation matrix: ( ), and then the guidewire is modeled according to the obtained three-dimensional position coordinates.

3. The method for calculating image force feedback of a guidewire based on discrete elastic rod physical modeling according to claim 2, characterized in that: Step S2 specifically includes: describing the guide wire as a coordinate three-dimensional curve: ;in It is a parameterized curve describing the arc length of the guidewire centerline in three-dimensional space. is an orthogonal material frame assigned to each point on the centerline, which contains all the information needed to calculate torsion and bending; the material frame satisfies Represents the tangent direction of the center line, and the two are perpendicular to each other. , Represents the The derivative of ; the curvature vector of the center line .

4. The method for calculating image force feedback of a guidewire based on discrete elastic rod physical modeling according to claim 3, characterized in that: Step S3 specifically includes: physically modeling the guidewire tip, and converting the energy Divided into three categories: , and , the total energy is: ； Since the guidewire is inextensible, the total energy is obtained by excluding the guidewire stretching energy and applying a PBD auxiliary constraint to enforce the guidewire inextensibility: ； The PBD constraint uses the Guass-Seidel method, which processes each edge separately and converges after more than 5 iterations, which can achieve non-scalability; The bending energy and torsional energy can be calculated by the following formula: ； Where L represents the guide wire length, represents the curvature vector of the curve, represents the first part of the curvature vector, represents the second part of the curvature vector, represents the inherent curvature of the guidewire tip, represents the flexural Young's modulus, represents the guidewire radius, represents the shear elastic modulus, represents the guidewire torsion, and D is the differential operator symbol; Perform numerical calculations and discretize the guidewire into particles for calculation. For bending force, first calculate the curvature of each discretized particle. , The calculation depends on two adjacent edges. In the case of continuity, the size is ,in Represents the angle between two sides, and the discrete numerical calculation uses the following formula: ； in, Represents the number corresponding to each particle after discretization, ; The discretization of bending energy is expressed as: ； in ; For the calculation of discrete torsional force, calculate the torsion angle of each point , According to the principle of minimum energy, we only need to know the torsion angle of the tip , use Newton iteration method to solve the following formula: ； in for matrix, is the gradient; is calculated as follows: ； in For rotation The rotation matrix of is a symmetric positive definite 2D matrix, The matrix H is calculated as follows: ； in represents the torsional modulus, The calculation method is: ； Finally updated all , The twist angle is obtained: ； The scalar , represents a discrete torsion.

5. The method for calculating image force feedback of a guidewire based on discrete elastic rod physical modeling according to claim 4, characterized in that: Step S4 specifically includes: using the Lagrange equation of motion to calculate the bending force of each particle of the guide wire and torsional force : ； Since the time step of the system is very short, there will be no numerical instability. The semi-implicit integration method is directly used to calculate the position change at the next moment. and speed change : ； Where M represents the quality, Represents the current speed, Represents the current position, Represents the current force.

6. The method for calculating image force feedback of a guidewire based on discrete elastic rod physical modeling according to claim 5, characterized in that: Step S5 specifically includes: in the coarse detection stage, directly obtaining whether the front end of the guide wire has collided with the blood vessel wall through voxels, and if the position Value , indicating that a collision has occurred, the guidewire motion is predicted based on the established physical model, and the speed, displacement, and force status of the guidewire at the next position are updated; if no collision occurs, the next step will continue to maintain the current state and move forward.

7. The method for calculating image force feedback of a guidewire based on discrete elastic rod physical modeling according to claim 6, characterized in that: Step S6 specifically includes: in the fine detection stage, the next step position is already known, and it is necessary to combine the blood vessel centerline and kd-tree to detect whether it penetrates the blood vessel wall. The kd-tree construction process will select the dimension with the largest variance each time. As a splitting axis: ； in Represents the mean of the i-th dimension data. There are k dimensions and n points in total. The randomness of the data distribution is guaranteed by maximizing the variance. Then, according to the median greed, the corresponding median of the selected dimension is selected. To divide: ； It represents taking the median. The greedy median ensures the balance of data volume in the left and right subtrees, keeps the tree height as low as possible, and achieves the log-level query distance to the nearest center line position.

8. The method for calculating image force feedback of a guidewire based on discrete elastic rod physical modeling according to claim 7, characterized in that: Steps S7 and S8 specifically include: Through the above integral formula and guidewire model, the next displacement information of the guidewire has been obtained. According to the relationship between the guidewire tip and the centerline of the blood vessel, it can be judged whether the guidewire tip will penetrate the blood vessel wall. If penetration occurs, the distance between the tip and the blood vessel wall can be found. ; Since different tissues have different elastic coefficients K and damping coefficients D, the simple spring damping model is: ； Combining the spring damping system with the structure of the blood vessel wall, the calculation model of the feedback force after collision can be obtained as follows: ； in represents the Young's modulus of the vascular wall, S represents the contact area between the vascular wall and the guidewire in the image, is the thickness of the blood vessel wall.

Images