Soft tissue cutting simulation method based on self-adaptive hexahedral mesh subdivision

Abstract

The invention provides a soft tissue cutting simulation method based on self-adaptive hexahedral mesh subdivision, which comprises the following steps of: S1, constructing a soft tissue model and a surgical instrument model by using three-dimensional modeling software; s2, performing collision detection on the soft tissue model and the surgical instrument model; s3, marking the initial collision position of the surgical instrument and the cutting model, and determining a cutting surface; s4, copying all hexahedrons and grid units intersected with the cutting surface, and redefining the connection relationship between the hexahedrons and the curved surface grids; according to the method, details of a complex cutting path can be caught more finely through self-adaptive hexahedron mesh subdivision, the sense of reality of cutting simulation is improved, in the cutting process, by copying or segmenting hexahedron units and updating quality and stiffness matrixes, it is ensured that the physical performance of the model is kept unchanged, unnecessary calculation overhead is reduced, and the method is suitable for large-scale popularization and application. And the simulation authenticity is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of virtual surgery simulation, and particularly to a soft tissue cutting simulation method based on adaptive hexahedron mesh subdivision. Background Art

[0002] Cutting techniques in soft tissue interaction simulation are widely used in fields such as medical simulation, virtual surgery, and virtual reality. In virtual surgery systems, common cutting algorithms include surface mesh-based cutting algorithms and polyhedron-based cutting methods.

[0003] In the prior art, common surface mesh-based cutting algorithms include the element removal method and the vertex splitting method, which are two algorithms commonly used to simulate the cutting effect of soft tissues. The element removal method, as a cutting algorithm proposed earlier by Bro-Nielsen et al., generates a cut by removing the primitives in contact with the surgical instrument in the soft tissue. This method is relatively simple to implement, but may appear less refined when dealing with complex cutting paths or when a smooth cutting surface needs to be maintained. The vertex splitting method is more refined. According to the cutting path of the interaction between the surgical instrument and the soft tissue, the cutting points of the model are split into two new vertices symmetric about the cutting path, and then these new vertices are connected with new connecting lines to simulate a smoother and more realistic cutting effect. To prevent the generation of low-quality reconstructed triangular patches during cutting, the vertex splitting algorithm often replaces the intersection points with new vertices to maintain a reasonable topological structure of the model after cutting. These two algorithms have their own advantages and disadvantages and are suitable for different application scenarios. The element removal method is more suitable for situations where the cutting accuracy requirement is not high or the computing resources are limited, while the vertex splitting method is more suitable for application scenarios that require high-precision cutting simulation and good visual effects.

[0004] Currently, hexahedron meshes are usually used for soft tissue deformation calculation, and surface meshes are embedded in the hexahedron meshes for visualization processing. Although this method has a certain degree of flexibility, the following problems also exist:

[0005] Insufficient capture of details of complex cutting paths: In the prior art, the generation of nodes in hexahedron meshes depends on established mesh division methods. When the cutting path is a curve or a non-planar path, due to the limited mesh resolution, the accuracy of the cutting path decreases, and the realism of the simulation results is poor. Especially in virtual surgery, it is difficult to accurately reproduce curved or fine anatomical cutting paths.

[0006] Efficiency and consistency issues in local mesh adjustment: To improve the cutting accuracy, increasing the overall mesh density will significantly increase the computational cost, and during the local mesh adjustment process, duplicate elements may be generated, resulting in non-conservation of the volume, mass, and stiffness of the model. Summary of the Invention

[0007] In view of this, in order to solve the problems of non-conservation of quality, volume, and stiffness in mesh reconstruction and non-smooth cutting paths in the technical background, the present invention proposes a soft tissue cutting simulation method based on adaptive hexahedron mesh subdivision. The specific technical solution is as follows:

[0008] A soft tissue cutting simulation method based on adaptive hexahedron mesh subdivision, comprising

[0009] Step S1: Use 3D modeling software to construct a soft tissue model and a surgical instrument model;

[0010] Step S2: Perform collision detection on the soft tissue model and the surgical instrument model;

[0011] Step S3: Mark the position where the surgical instrument and the cutting model have an initial collision and determine the cutting surface;

[0012] Step S4: Copy all hexahedrons and mesh cells that intersect the cutting surface, and redefine the connection relationship between the hexahedrons and the surface meshes.

[0013] The specific operation of constructing the soft tissue model in step S1 is as follows:

[0014] The soft tissue model is composed of regular hexahedron elements and surface mesh cells. The surface mesh structure is embedded on the surface of the regular hexahedron. The vertices of the surface mesh are approximately obtained from the positions of the vertices of the regular hexahedron, and the model boundary is represented by the surface mesh.

[0015] The specific operation of the collision detection in step S2 is as follows:

[0016] S21: The contact between the soft tissue model and the surgical instrument is simulated in two steps. First, project the surface edges intersecting the tool onto a plane, then assign the constraint conditions to the tool end, and use the tool's axis-aligned bounding box (AABB) for rough collision testing to find potential collision pairs; the margins of the AABB are the average lengths of the surface along the X, Y, and Z axes in the initial state, and then perform fine collision testing on each collision pair to obtain the intersecting surface edges.

[0017] S22: Use the least squares method for surface representation. The calculation formula for the i th th surface vertex is:

[0018]

[0019] where is the initial position, is a group of adjacent hexahedron nodes, is the midpoint of a group of adjacent hexahedron nodes, ξ ij and ξ i′jis the moving least squares approximation function, u j is the position of the j-th node.

[0020] In the projection step of step S21, the initial positions of each intermediate surface vertex are updated as follows:

[0021]

[0022] where is the updated position, x0 s is the initial position, p 1 is the position of a vertex of the cutting tool, and n is the normalized normal of the cutting tool.

[0023] The projected surface edges will be parallel to the cutting tool, so continuous collisions may not be detected correctly. This problem can be solved by using center edges. Connect the center points of the mesh cells to the midpoints of each edge. The connecting line is called the center edge. A one-to-one matching can be achieved between the center edge and the surface edge. If a certain center edge collides with the cutting tool, then the corresponding surface edge is regarded as having collided.

[0024] A soft tissue cutting simulation method based on adaptive hexahedral mesh subdivision, characterized in that the specific operation of the cutting surface in step S3 is as follows:

[0025] The cutting tool is simplified to a quadrilateral. During the cutting process, the surface mesh of the deformable model is presented in wireframe form, and the intersecting edges and intermediate surface vertices are presented as thicker lines and points respectively.

[0026] The cutting surface is represented by generating additional surface meshes inside the model. This area is defined by the initial position x0 of the mid-vertex s and the current position First, additional surface vertices A are generated along the tip of the tool. The negative cutting surface is represented by the vertices in S and A, while the positive cutting surface is represented by the vertices in S' and A. These two surfaces have common edges and vertices at the tip of the tool. The determination of the positive and negative cutting surfaces is as follows:

[0027]

[0028] where N f is a set of triangular faces around a vertex or an edge, m is the number of elements in N f is the center position of face i th p 1 is the position of a vertex of the cutting tool, and n is the normalized normal direction of the tool.

[0029] ​A method for simulating soft tissue cutting based on adaptive hexahedral mesh subdivision, characterized in that there are three special cases to be considered during the cutting process of the cutting surface. The first case is when a tip of the cutting tool is inside the deformable model, add this tip to A to handle this case. The second case is when the model is completely cut, the two sides of the cutting surface must be separated from each other. In this case, each cutting surface is only represented by the vertices in S and S′, and there is no need to generate an additional group A. The third case is that multiple Ss may appear when cutting a concave geometric body. i For the third case, it is handled by considering each group separately.

[0030] A method for simulating soft tissue cutting based on adaptive hexahedral mesh subdivision, characterized in that the specific operation of the replication unit in step S4 is as follows:

[0031] S41. Define the volume ratio of all hexahedral elements, denoted as α ∈ R mx1 , where m is the number of hexahedrons, the value of a complete hexahedron is 1, and the value of a partial hexahedron is between 0 and 1. When the hexahedron intersects the cutting path completely, execute the element replication process. All the regular hexahedrons intersected by the cutting path are replicated into two divided parts, and the replicated intersecting edges and vertices are recorded as S and S′. S and S′ represent the initial group and the replicated group respectively, and S is the positive side of the surgical instrument.

[0032] S42. A single hexahedral unit x k is divided into two parts by the cutting plane. The cutting plane is represented by the normal vector n and a vertex of the plane. Generate n 3 distribution points inside all the cut regular hexahedrons. The distribution points are located on both sides of the cutting plane respectively, and both sides of the cutting plane can be determined by the formula

[0033] (P - x i ) · n (4)

[0034] where x i is the position of the i-th distribution point, and the volume ratios of the two separated hexahedrons are n neg / n 3 and n pos / n 3 respectively, where n neg and n pos are the number of points on the positive and negative sides of the plane respectively. When , that is, the number of negative side faces is less than the number of points of the original hexahedron then only keep the replicated hexahedron on the positive side and delete the hexahedron on the negative side face The quantity is where x k is the center point of the hexahedron where the vertex is located, M fFor the hexahedron near the vertex, when the condition is not satisfied, that is, the number of points on both sides of the plane is close, then both hexahedrons are retained and recorded as . For the global volume of the model, keeping the quantities of m and m' equal can maintain the overall conservation of the model;

[0035] The calculation of the hexahedron mass matrix is expressed as

[0036]

[0037] where ρ is the mass density, N is the shape function, and Ω e is the element domain, and M e represents the mass moment.

[0038] The calculation of the hexahedron stiffness matrix is calculated by using the shape function and the constitutive relationship of the elastic material, and is expressed as

[0039]

[0040] where Ω e is the hexahedron element volume region, B is the derivative of the shape function N, E is the elastic matrix, and Ω e is the element domain.

[0041] S43. Calculate the stiffness matrix of the nodal element:

[0042]

[0043] where N i is a group of hexahedrons around the i-th node, represents the number of hexahedrons around the i-th node, and K e is the hexahedron stiffness matrix, T1 is the nodal domain. Considering the influence of mass on the calculation of the nodal element stiffness matrix, there is:

[0044]

[0045] where N i is a group of hexahedrons around the i-th node, where x k is the center point of the hexahedron where the vertex is located, represents the number of hexahedrons around the i-th node, and K e is the hexahedron stiffness matrix, T1 is the nodal domain, where x k is the center point of the hexahedron where the vertex is located, and M f is the hexahedron near the vertex.

[0046] The nodal unit stiffness term of the model is represented by f i . The calculation of the overall stiffness term of the model:

[0047]

[0048] is the stability matrix, f′ i is the stiffness matrix of the nodal element, x j and are respectively the current position and the stationary position of the j-th th hexahedral element. α is the volume ratio of the hexahedral element. The mass and stiffness matrices of the entire model must be updated according to α. The mass and stiffness matrices of the entire model must be updated according to the volume ratio α of the hexahedral element, and the overall conservation is maintained.

[0049] Adopting the above technical solution, the following beneficial effects are achieved:

[0050] Through adaptive hexahedral mesh subdivision, the present invention can capture the details of complex cutting paths more precisely, improve the realism of cutting simulation. During the cutting process, by replicating or splitting hexahedral elements and updating the mass and stiffness matrices according to the volume ratio, it ensures that the overall performance of the model remains unchanged. By optimizing the collision detection and cutting surface determination algorithms, unnecessary computational overhead is reduced, and the real-time performance of the simulation is improved. Description of the Drawings

[0051] Figure 1 is a flowchart of a soft tissue cutting method based on adaptive hexahedral mesh subdivision according to the present invention. Detailed Embodiments

[0052] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0053] Embodiment 1: Refer to Figure 1 As shown, a soft tissue cutting method based on adaptive hexahedral mesh subdivision includes: Step S1, using 3D modeling software to construct a soft tissue model and a surgical instrument model; Step S2, performing collision detection on the soft tissue model and the surgical instrument model; Step S3, marking the position where the surgical instrument and the cutting model initially collide and determining the cutting surface; Step S4, replicating all hexahedrons and mesh cells intersecting with the cutting surface and redefining the connection relationship between the hexahedrons and the surface meshes;

[0054] The specific operation of constructing the soft tissue model according to Step S1 is as follows:

[0055] The soft tissue model is composed of hexahedron elements and curved surface mesh elements. The curved surface mesh structure is embedded on the surface of the hexahedron. The vertices of the curved surface mesh are approximately obtained from the positions of the vertices of the hexahedron. The model boundary is represented by the curved surface mesh;

[0056] The specific operation of collision detection described in step S2 is as follows:

[0057] S21. The contact between the soft tissue model and the surgical instrument is simulated in two steps. First, the curved surface edges intersecting with the tool are projected onto a plane, and then the constraint conditions are assigned to the tool end. The axial bounding box (AABB) of the tool is used for rough collision testing to find potential collision pairs; the margins of the AABB are the average lengths of the surface along the X, Y, and Z axes in the initial state. Then, fine collision testing is performed on each collision pair to obtain the intersecting surface edges.

[0058] S22. The least squares method is used for surface representation. In the projection step, the initial positions of each intermediate surface vertex are updated. After projection, the surface edges will be parallel to the cutting tool, so continuous collisions may not be correctly detected. This problem can be solved by using center edges. The center point of the grid cell is connected to the midpoint of each edge, and this connecting line is called the center edge. A one-to-one match can be achieved between the center edge and the surface edge. If a certain center edge collides with the cutting tool, then the corresponding surface edge is considered to have collided.

[0059] The specific operation for determining the cutting surface is as follows:

[0060] The cutting tool is simplified to a quadrilateral. During the cutting process, the surface mesh of the deformable model is presented in wireframe form, and the intersecting edges and intermediate surface vertices are presented as thicker lines and points respectively.

[0061] The cutting surface is represented by generating additional surface meshes inside the model. This area is defined by the initial position x0 of the midpoint s The current position x s First, additional surface vertices A are generated along the tip of the tool. The negative cutting surface is represented by the vertices in S and A, while the positive cutting surface is represented by the vertices in S′ and A. These two surfaces have common edges and vertices at the tip of the tool.

[0062] There are three special cases to consider during cutting. The first case is when one tip of the cutting tool is inside the deformable model. In this case, this tip is added to A to handle this situation. The second case is when the model is completely cut. The two sides of the cutting surface must be separated from each other. In this case, each cutting surface is only represented by the vertices in S and S′, and no additional set A needs to be generated. The third case is that multiple Ss may appear when cutting a concave geometry. iGroup. For the third case, it is handled by considering each group separately.

[0063] The specific operations of the replication unit described in step S4 are as follows:

[0064] S41. Define the volume ratio of all hexahedron elements, denoted as α ∈ R mx1 , where m is the number of hexahedrons. The value of a complete hexahedron is 1, while the value of a partial hexahedron is between 0 and 1. When the hexahedron intersects the cutting path completely, the element replication process is executed. All the regular hexahedrons intersected by the cutting path are replicated into two divided parts, and the replicated intersecting edges and vertices are recorded as S and S′. S and S′ represent the initial group and the replicated group respectively, and S is the positive side of the surgical instrument.

[0065] S42. A single hexahedron unit x k is divided into two parts by the cutting plane. The cutting plane is represented by the normal vector n and a vertex of the plane. n 3 distribution points are generated inside all the cut regular hexahedrons, and the distribution points are located on both sides of the cutting plane.

[0066] The volume ratios of the two separated hexahedrons are n neg / n 3 and n pos / n 3 , where n neg and n pos are the number of points on the positive and negative sides of the plane respectively. When , that is, the number of negative side faces is less than the number of points of the original hexahedron , then only the replicated hexahedron on the positive side is retained, and the hexahedron on the negative side face is deleted The quantity is where x k is the center point of the hexahedron where the vertex is located, and M f is the hexahedron near the vertex. When the condition is not met, that is, the number of points on the positive and negative sides of the plane is close, then both hexahedrons are retained and recorded as For the global volume of the model, keeping the quantities of m and m′ equal can maintain the overall conservation of the model.

[0067] Calculate the hexahedron mass matrix of the model by using the shape function and the constitutive relation of the elastic material.

[0068] S43. After executing the replication of the hexahedron unit, the mass and stiffness matrices of the whole model must be updated to maintain the overall physical mass and stiffness conservation of the model.

[0069] The basic principles and main features of the present invention have been described above. Those skilled in the art should understand that the present invention is not limited by the above embodiments. What is described in the above embodiments and the specification only illustrates the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention will have various changes and improvements, and these changes and improvements all fall within the scope of the present invention claimed. The scope of the invention claimed is defined by the appended claims and their equivalents.

Claims

1. A soft tissue cutting simulation method based on adaptive hexahedron mesh subdivision, characterized in that, including Step S1: Use 3D modeling software to construct a soft tissue model and a surgical instrument model; Step S2: Perform collision detection on the soft tissue model and the surgical instrument model; Step S3: Mark the position where the surgical instrument and the cutting model initially collide and determine the cutting surface; Step S4: Copy all hexahedrons and mesh elements that intersect the cutting surface and redefine the connection relationship between the hexahedrons and the surface mesh. The specific operation of constructing the soft tissue model in Step S1 is as follows: The soft tissue model is composed of hexahedron elements and surface mesh elements. The surface mesh structure is embedded on the surface of the hexahedron. The vertices of the surface mesh are approximately obtained from the vertex positions of the hexahedron. The model boundary is represented by the surface mesh; The specific operation of the collision detection in Step S2 is as follows: S21: The contact between the soft tissue model and the surgical instrument is simulated in two steps. First, project the surface edges intersecting the cutting tool onto a plane, and then assign constraint conditions to the tool end. Use the tool's axis-aligned bounding box (AABB) for a rough collision test to find potential collision pairs; the margins of the AABB are the average lengths of the surface along the X, Y, and Z axes in the initial state. Then, perform a fine collision test on each collision pair to obtain the intersecting surface edges; S22. Use the least squares method for surface representation. The calculation formula for the th ith surface vertex is as follows: Among them, is the initial position, is a group of adjacent hexahedral nodes, is the middle node of a group of adjacent hexahedral nodes, ξ ij and ξ i′j are moving least squares approximation functions, u j is the position of the j-th node; In the projection step of Step S21, the initial position of each intermediate surface vertex is updated as follows: Among them, is the updated position, x0 s is the initial position, p 1 is the vertex position of a cutting tool, and n is the normalized normal of the cutting tool; The projected surface edges will be parallel to the cutting tool, so continuous collisions may not be correctly detected. This problem can be solved by using center edges. Connect the center points of the mesh elements to the midpoints of each edge. This connecting line is called the center edge. A one-to-one match can be achieved between the center edge and the surface edge. If a certain center edge collides with the cutting tool, then the corresponding surface edge is considered to have collided.

2. The soft tissue cutting simulation method based on adaptive hexahedral mesh subdivision according to claim 1, wherein The specific operation of the cutting surface in Step S3 is as follows: Simplify the cutting tool into a quadrilateral. During the cutting process, the surface mesh of the deformable model is presented in wireframe form, and the intersecting edges and intermediate surface vertices are presented as thicker lines and points respectively; The cutting surface is represented by generating additional surface meshes inside the model, and this area is defined by the initial position x0 of the mid-vertices s and the current position First, additional surface vertices A are generated along the tip of the tool. The negative cutting surface is represented by the vertices in S and A, while the positive cutting surface is represented by the vertices in S' and A. These two surfaces have common edges and vertices at the tip of the tool. The determination of the positive and negative cutting surfaces is as follows: where N f is a set of triangular faces surrounding a vertex or an edge, m is the number of elements in N f , is the center position of face i th , p 1 is a vertex position of the cutting tool, and n is the normalized normal direction of the tool.

3. A soft tissue cutting simulation method based on adaptive hexahedron mesh subdivision according to claim 1, characterized in that, There are three special cases to consider during the cutting process of the cutting surface: One case is when a tip of the cutting tool is inside the deformable model, and this tip is added to A to handle this case; the second case is when the model is completely cut, and the two sides of the cutting surface must be separated from each other. In this case, each cutting surface is represented only by the vertices in S and S′, and no additional group A needs to be generated; the third case is that multiple S i groups may occur when cutting a concave geometric body. For the third case, it is handled by considering each group separately.

4. A soft tissue cutting simulation method based on adaptive hexahedron mesh subdivision according to claim 1, characterized in that, The specific operation of copying the elements in Step S4 is as follows: S41. Define the volume ratio of all hexahedron elements, denoted as α ∈ R mx1 , where m is the number of hexahedrons. The value of a complete hexahedron is 1, while the value of a partial hexahedron is between 0 and 1. When the hexahedron intersects the cutting path completely, the element replication process is executed. The regular hexahedron intersected by the cutting path is replicated into two divided parts, and the replicated intersecting edges and vertices are recorded as S and S′. S and S′ represent the initial group and the replicated group respectively, and S is the positive side of the surgical instrument; S42. Single hexahedron unit x k is divided into two parts by the cutting plane, and the cutting plane is represented by the normal vector n and a vertex of the plane. Generate n 3 distribution points within all the cut regular hexahedrons. The distribution points are located on both sides of the cutting plane, and both sides of the cutting plane can be determined by the formula as follows: (P - x i )·n (4) where x i is the position of the i-th distribution point, and the volume ratios of the two separated hexahedrons are n neg / n 3 and n pos / n 3 , where n neg and n pos are the number of points on the positive and negative sides of the plane respectively. When , that is, the number of negative side faces is less than the number of points of the original hexahedron then only the replicated hexahedron on the positive side is retained, and the hexahedron on the negative side is deleted The quantity is where x k is the center point of the hexahedron where the vertex is located, and M f is the hexahedron near the vertex. When the condition is not met, that is, the number of points on the positive and negative sides of the plane is close, then both hexahedrons are retained and recorded as For the global volume of the model, keeping the quantities of m and m' equal can maintain the overall conservation of the model; The calculation of the hexahedron mass matrix of the model is expressed as: where ρ is the mass density, N is the shape function, and Ω e is the element domain, and M e represents the mass moment; The calculation of the hexahedron stiffness matrix is calculated by using the shape function and the constitutive relationship of the elastic material, and is expressed as: where Ω e is the volume region of a hexahedral element, B is the derivative of the shape function N, E is the elastic matrix, and Ω e is the element domain; S43: Calculate the nodal element stiffness matrix where N i is a group of hexahedrons around the i-th node, represents the number of hexahedrons around the i-th node, K e is the hexahedron stiffness matrix, T1 is the node domain. Considering the influence of the quantity on the calculation of the node element stiffness matrix, there are: Among which N i is a group of hexahedrons around the i-th node, where x k is the center point of the hexahedron where the vertex is located, represents the number of hexahedrons around the i-th node, K e is the hexahedron stiffness matrix, T1 is the node domain, where x k is the center point of the hexahedron where the vertex is located, M f is the hexahedron near the vertex; The nodal unit stiffness term of the model is represented by f i The calculation of the overall stiffness term of the model is as follows: is the stability matrix, f i ′ is the stiffness matrix of the nodal element, x j and are respectively the current position and the static position of the j th -th hexahedral element, α is the volume ratio of the hexahedral element. The mass and stiffness matrices of the whole model must be updated according to α. The mass and stiffness matrices of the whole model must be updated according to the volume ratio α of the hexahedral element, and the overall conservation is maintained.

Images