A virtual surgery simulation method based on grid adaptive cutting and subdivision

Abstract

The present invention discloses a virtual surgery simulation method based on grid adaptive cutting and subdivision. First, the input object to be cut is tetrahedralized, and the position of the scalpel is recorded to form a cutting surface. Secondly, a bounding box and a bsp tree A are constructed, and the splicing surface and the simulation surface are marked. The common cutting line is generated by triangulation to obtain a simulated cutting surface. Then the movable vertices and the immovable vertices are marked, and the positions of the new vertices are solved according to Hooke's law. The vertices are moved to simulate the shape after cutting, and the local area is re-tetrahedralized and bsp subdivided to generate a bsp tree B and A merged. Finally, all bsp nodes are traversed, and the non-tetrahedral nodes are tetrahedralized using the center of gravity of the non-tetrahedron, the simulated cutting is completed, and optimization is performed. The present invention accurately cuts the object to be cut, improves the accuracy of the cutting plane, and ensures the stability and accuracy of the cutting simulation.

Description

Technical Field

[0001] The invention belongs to the technical field of virtual surgery cutting simulation, and in particular relates to a virtual surgery simulation method based on grid adaptive cutting and subdivision. Background Art

[0002] Virtual surgery technology provides doctors with an interactive surgical environment and a comprehensive operating platform, which is widely used in many fields such as surgical planning, training and simulation teaching. Compared with traditional surgical training and teaching methods, virtual surgery is not only low-cost but also repeatable, which makes it more and more used in medical education and practice. With the continuous advancement of virtual reality technology and the continuous improvement of computer hardware performance, virtual surgery has gradually become a hot research field. Among them, the cutting simulation of deformable objects has been regarded as a basic and key technology, which requires high accuracy and real-time performance to ensure the correctness of surgical simulation results.

[0003] At present, most cutting simulation methods for deformable objects are based on mesh construction, among which the cutting simulation method based on tetrahedral mesh is the most widely used one. However, when fitting the cutting plane through mesh subdivision, tetrahedral mesh may have errors and ill-conditioned units. On the one hand, the inaccuracy of the cutting plane will significantly affect the accuracy of the physical simulation; on the other hand, the generation of ill-conditioned units will reduce the computational efficiency and stability of the deformation simulation. Therefore, how to effectively process the mesh and control its quality has always been a research hotspot in the field of cutting simulation. Summary of the invention

[0004] The present invention aims to solve the accuracy problem of deformable object cutting simulation in virtual surgery and the generation of pathological units, thereby improving the mesh quality of the model and enhancing the efficiency and stability of cutting simulation. To this end, the present invention proposes a virtual surgery simulation method based on mesh adaptive cutting and subdivision. The method adopts a new cutting surface generation technology, solves the accuracy problem of cutting surface generation, and effectively reduces the generation of pathological units.

[0005] In order to solve the above technical problems, the technical solution of the present invention is:

[0006] A virtual surgery simulation method based on grid adaptive cutting and subdivision comprises the following steps:

[0007] Step 1: When the input object to be cut is a non-tetrahedral mesh, it is re-tetrahedralized to obtain a tetrahedron corresponding to the object to be cut;

[0008] Step 2: Use line segments to simulate a scalpel. Each time the scalpel moves for a time step, the position of the scalpel is recorded. The starting position and the position after moving for a time step together constitute the cutting surface on the tetrahedron.

[0009] Step 3: construct a suitable bounding box according to the cutting surface and the maximum side length of the tetrahedron, and construct a fork space partition bsp tree A for the input tetrahedron model of the cut object;

[0010] Step 4: Use the bounding box and winding number calculation to mark the BSP nodes involved in the bounding box, and mark the splicing surface and simulation surface;

[0011] Step 5: Taking the simulated surface and the cutting surface as input, a constrained Delaunay triangulation method is used to generate a common cutting line, and then the common cutting line and the cutting surface are preprocessed by physical simulation to form a surface that can be simulated, thereby obtaining a simulated cutting surface;

[0012] Step 6. Construct a mass-spring model for the simulated cutting surface: mark the movable vertices and immovable vertices respectively, then use Newton's iteration method according to Hooke's law to solve the position of the new vertex, and move the vertex to simulate the shape after cutting.

[0013] Step 7: Re-tetrahedralize and BSP subdivide the local area based on the simulated shape, generate a new BSP tree B, and merge BSP trees A and B;

[0014] Step 8: After merging, traverse all BSP nodes, and for non-tetrahedral nodes, use the center of gravity of the non-tetrahedron to re-tetrahedralize to complete the simulated cutting;

[0015] Step 9, optimizing the tetrahedron model after simulated cutting;

[0016] In the step 1, when the input data is a triangular mesh surface in obj format, the tetgen algorithm is automatically executed to generate a tetrahedral mesh.

[0017] In step 2, the mouse click point and click duration are obtained through the qt interactive component, and a corresponding starting line segment is generated according to the mouse click point and its direction, and its length depends on the click duration. The mouse is clicked again to obtain the end line segment, and its length is equal to that of the starting line segment.

[0018] In the step 2, four endpoints A, B, C, and D constituting the starting line segment AB and the ending line segment CD are obtained, and a cutting plane consisting of two triangles is generated using the four endpoints, where the triangles are ABC and ACD.

[0019] In step 3, since a bounding box of a suitable size needs to be constructed to extract a local area that meets the requirements, the side length of the tetrahedron and the length of the cutting plane need to be considered.

[0020] The specific method of step 3 is:

[0021] Step 3.1: By traversing the edge lengths of all tetrahedrons, get the maximum edge length len_ A ; Traverse the edge lengths on the cutting plane and get the maximum edge length len_ of the cutting plane B . Take the maximum side length of the tetrahedron len_ A and the maximum side length of the cutting plane len_ B The maximum value in len_ C , and then give len_ C Multiply it by an appropriate factor and use it as the offset value of the diagonal of the bounding box;

[0022] Step 3.2: According to the positions of the four vertices of the cutting plane ABCD, calculate the maximum value of these four points in the three-dimensional coordinate axis (x max ,y max , z max ) and minimum value (x min ,y min , z min ).

[0023] where x max ,y max , z max The maximum values ​​of the four vertices ABCD in the x, y, and z axes, x min ,y min , z min is the minimum value;

[0024] Step 3.3, the bounding box is composed of 8 vertices, and its diagonal is (x max + offset,y max +offset,z max + offset) and (x min - offset, y min -offset,z min -offset), so the 8 vertices are (x max + offset,y max + offset, z max + offset), (x max + offset,y min - offset, z max +offset)、(x max + offset,y max + offset, z min - offset), (x max + offset,y min - offset, zmin -offset)、(x min - offset, y max + offset, z max + offset)、(x min - offset, y max + offset, z min -offset)、(x min - offset, y min - offset, z max + offset)、(x min - offset, y min - offset, z min -offset);

[0025] Step 3.4, traverse the cut tetrahedron, where each triangular surface is regarded as a plane, and each tetrahedron is regarded as a leaf node of a bsp unit. In this way, the bsp tree A can be created;

[0026] The specific method of step 4 is:

[0027] Step 4.1, use Jacobson algorithm to construct the inner-outer function of the constructed bounding box;

[0028] Step 4.2, traverse all BSP nodes and calculate the center of gravity of each node;

[0029] Step 4.3, use the inner-outer function to determine whether the centroid is contained in the bounding box. If it is in the bounding box, mark the bsp node. If it is not in the bounding box, do not process it.

[0030] Step 4.4, traverse all marked bsp nodes, if the face in the node is adjacent to an unmarked bsp node, and the face is not the surface of the cut object, then mark the face as a splicing face;

[0031] Step 4.5, traverse all the faces in the marked bsp nodes, if the face is the surface of the cutting object, mark the face as a simulation surface;

[0032] The specific method of step 5 is:

[0033] Step 5.1, traverse the bsp nodes and extract all simulation surfaces;

[0034] Step 5.2, use the simulation surface to intersect with the cutting surface to generate constrained edges and points, and then perform constrained Delaunay triangulation to generate common cutting lines;

[0035] Step 5.3, split the common cutting line into two: find the starting endpoint of the common cutting line, and use this endpoint as the initial vertex A1, and the vertex on the next common cutting line connected to it is called B1. Randomly divide the two adjacent faces of the A1B1 edge into two categories, such as red and white, and generate point B1' as the symmetric vertex of B1, then use the face in the red category as the starting face, rotate around point B1, and all the faces it touches are classified as red until it touches the common cutting line B1C1, and change all B1 points in the red category to B1', and so on until it touches the end of the cutting line;

[0036] Step 5.4, splice the common cutting line and the cutting surface to form a simulated cutting surface: copy the points on the cutting surface except the points on the boundary and not on the common cutting line, and form a symmetric surface of the new cutting surface with the newly generated common cutting line and the points on the boundary, and give a label to the simulated cutting surface;

[0037] The specific method of step 7 is:

[0038] Step 7.1, extract all nodes contained in the simulation surface and the spliced ​​surface;

[0039] Step 7.2, delete all bsp nodes in the original bsp tree A that are included in the local area;

[0040] Step 7.3, use Delaunay to re-tetrahedralize all nodes;

[0041] Step 7.4, construct a new bsp tree structure B for the tetrahedron and its contained surfaces;

[0042] Step 7.5, using the BSP subdivision method, restore the re-tetrahedralized feature surface;

[0043] Step 7.6, merge bsp tree A and bsp tree B, and replace the original splicing surface with all new sub-surfaces;

[0044] In step 8, the center of gravity is calculated by taking the average of all node positions contained in the non-tetrahedron. All non-tetrahedrons are convex structures, and because they are convex structures, their calculated center of gravity must be inside the non-tetrahedron.

[0045] In the step 8, because the re-tetrahedralization with the center of gravity is performed, the adjacent tetrahedral structures will not be affected.

[0046] The specific method of step 9 is:

[0047] Step 9.1: Construct an aabb tree structure for the tetrahedron surface.

[0048] Step 9.2: Build a minimum heap based on the edge lengths of all tetrahedrons, extract the edge with the smallest length, and merge it if the edge length is less than the minimum threshold set by the user.

[0049] Step 9.3: If the two nodes V1 and V2 of the edge are not surface nodes, replace all tetrahedrons containing V2 nodes with V1 nodes, and translate the position of V1 nodes to (V1+V2) / 2.

[0050] Step 9.3: If the edge contains a surface node, use the aabb tree to determine whether the distance (V1+V2) / 2 from the surface is within the deviation set by the user. If so, merge them.

[0051] Step 9.4: When there are no edges in the heap, the edge length optimization algorithm ends.

[0052] Step 9.5: Execute the quadrilateral mesh-based Laplace smoothing method on all nodes to complete node optimization.

[0053] The present invention has the following characteristics and beneficial effects:

[0054] By adopting the above technical solution, the present invention can accurately cut the object to be cut while ensuring the quality of the generated tetrahedral mesh. The present invention fully considers the high accuracy requirements of virtual surgery and simplifies the traditional complex tetrahedron cutting process into a cutting operation for triangular surfaces, greatly improving the accuracy of the cutting plane.

[0055] In addition, the present invention focuses on solving the problem of low tetrahedron quality after simulation for the simulation scene of cut objects in virtual surgery. The method mainly performs physical simulation, re-tetrahedronization, BSP subdivision and post-processing on the local cut area to effectively solve the quality problems such as tetrahedron unit inversion that may occur during the simulation process, ensuring the stability and accuracy of the cutting simulation. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use are briefly introduced below. Obviously, the drawings in the following description are only some specific embodiments of the present invention, and technicians familiar with the field can also obtain other related diagrams based on these drawings without creative work.

[0057] Figure 1 is a flow chart of the present invention;

[0058] Figure 2 The results of triangular mesh and tetrahedronization, where (a) is the triangular mesh and (b) is the result after tetrahedronization;

[0059] Figure 3is the cutting plane;

[0060] Figure 4 Schematic diagram of the bounding box construction process, where (a) is the location of the cutting plane, (b) is the bounding box constructed for the cutting plane, and (c) is the cutting plane and the bounding box;

[0061] Figure 5 Schematic diagram of constructing a bounding box to extracting a local area, where (a) is the cut object and the cut surface, (b) is the shape after the bounding box is constructed, (c) is the extracted local area and the remaining cut object shape, and (d) is the extracted local area;

[0062] Figure 6 Schematic diagram of the process of generating a common cutting line from a cutting surface, wherein (a) is the shape of the cutting surface and the object to be cut before the cutting line is formed, and (b) is the shape of the cutting surface and the object to be cut after the cutting line is formed;

[0063] Figure 7 A schematic diagram of a cutting surface generated by a cutting surface and capable of physical simulation;

[0064] Figure 8 Schematic diagram of the common cutting line breaking process, where (a) is the result after the initial division into red and white categories, (b) is the result of the category division after the symmetry point is formed, (c) is the result of the division after the breaking is completed and the result after the symmetry point is generated, and (d) is the rendering result after the common cutting line is broken;

[0065] Fig. 9 Schematic diagram of the physical simulation preprocessing process, where (a) is the object to be cut and the cutting surface, and (b) is the sample after the fixed point is formed;

[0066] Fig.10 It is a schematic diagram of the subdivision of the local cutting area bsp;

[0067] Fig.11 Schematic diagram of the center-of-gravity-based re-tetrahedralization, where (a) shows the shape of the non-tetrahedron before re-tetrahedralization, and (b) shows the shape of the non-tetrahedron ABECD after it is subdivided into four tetrahedrons: ABDG, BDGE, EGDC, and ADCG;

[0068] Fig.12 Schematic diagram of barycentric-based re-triangulation, where (a) is a picture of the triangular surface before subdivision, (b) is the non-triangular surface BDCE generated after subdivision, and (c) is the shape of the non-triangular surface BDCE after re-triangulation based on the barycentric;

[0069] Fig.13 It is a schematic diagram of the final physical simulation results;

[0070] Fig.14 Schematic diagram of the optimization results, where (a) is the tetrahedron shape before optimization and (b) is the tetrahedron shape after optimization. DETAILED DESCRIPTION

[0071] It should be noted that, in the absence of conflict, the embodiments of the present invention and the features in the embodiments may be combined with each other.

[0072] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings.

[0073] The present invention covers any substitution, modification, equivalent method and scheme made on the essence and scope of the present invention as defined by the claims. Further, in order to make the public have a better understanding of the present invention, some specific details are described in detail in the following detailed description of the present invention. Those skilled in the art can fully understand the present invention without the description of these details.

[0074] Example 1

[0075] This embodiment provides a virtual surgery simulation method based on grid adaptive cutting and subdivision. Figure 1 As shown, the specific steps are as follows:

[0076] Step 1: Re-tetrahedralize the input non-tetrahedral mesh;

[0077] The present invention determines whether it is a tetrahedron by judging the input format. Figure 2 As shown, (a) is the input triangular mesh, and (b) is the tetrahedral mesh generated by tetgen.

[0078] Step 2: Use line segments to simulate a scalpel and form a cutting surface;

[0079] The cutting surface formed by two triangles is as follows Figure 3 shown.

[0080] Step 3: construct a suitable bounding box and build a bsp tree A for the input cut object;

[0081] The construction of the bounding box is to extract enough area for normal simulation, so it is necessary to refer to the edge length of the tetrahedron and the edge length of the cutting surface. Figure 4 As shown, (a) is the position of the cutting surface and the cut object, (b) is the generated bounding box area of ​​appropriate size, and (c) is a schematic diagram of the bounding box.

[0082] For the construction of bsp tree A, because it is constructed based on tetrahedralized grid, and tetrahedron itself is a convex polyhedron, it meets the definition of bsp node. Each face of tetrahedron can be regarded as a plane, and its adjacency relationship is equivalent to the adjacency relationship of the face, so the construction of bsp tree can be regarded as a simple assignment operation.

[0083] Step 4: Use the bounding box and winding number calculation to mark the BSP nodes involved in the bounding box, and mark the splicing surface and simulation surface;

[0084] like Figure 5 As shown, (b) is the bounding box constructed using (a). The bounding box is used to calculate the winding number and extract the tetrahedrons contained therein, namely the local simulation area. Each tetrahedron corresponds to a bsp node, and the bsp node can be found through the tetrahedron ID. Figure 5 (c) is the extracted local area and the remaining shape of the cut object, (d) is the local simulation area, which is obtained by traversing all the tetrahedral faces in it. If the face is only adjacent to one tetrahedron, then the face is the local area surface. If the face was previously marked as the cut object surface, then the face is marked as a simulation surface, such as the green area in (d). Otherwise, it is marked as a splicing surface, such as the yellow and red areas in (d), and the face is identified in the bsp tree.

[0085] Step 5: With the simulation surface and cutting surface as input, a constrained Delaunay triangulation method is used to generate a common cutting line, and then the common cutting line and cutting surface are preprocessed by physical simulation to form a surface that can be simulated, and together they are constructed into a simulatable surface.

[0086] The common cutting line is generated by the previously obtained simulated surface and the cutting surface using the constrained Delaunay triangulation method, such as Figure 6 In order to form Figure 7 The surfaces shown can be simulated, requiring the common cutting lines to be disconnected.

[0087] By breaking the common cutting line into two: find the starting endpoint of the common cutting line, and use this endpoint as the initial vertex A. The vertex on the next common cutting line connected to it is called B, such as Figure 8 As shown in (a) in the figure. The two adjacent faces of the AB edge are randomly divided into two categories, such as red and white, and point B' is generated. And B' is the symmetric vertex of B. Then the face in the red category is used as the starting face and rotated around point B. During the rotation process, all the faces it touches are classified as red, until it touches the common cutting line BC, and all the B points in the red category are changed to B', as shown in Figure 8 As shown in (b) in the figure. And so on until the end of the cutting line is reached. Figure 8As shown in (c) in the figure, the common cutting line can be disconnected in this way. Figure 8 (d) shows the example result diagram.

[0088] Then combine the two common cutting lines with the cutting surface to generate a cutting surface that can be used for simulation: copy the points on the cutting surface that are not on the common cutting line except the boundary points, and combine them with the newly generated common cutting line and boundary points to form a new cutting surface that is symmetrical on the other side, and give the generated surface a simulation surface label, such as Figure 7 shown.

[0089] Step 6: Construct a mass-spring model for the simulated surface: mark the movable vertices and the immovable vertices respectively, apply gravity and resistance to the movable vertices, and then use the Newton iteration method to calculate the position of the vertices to simulate the cut shape, and mark the simulated surface as the outer surface of the local cutting area;

[0090] Based on the experience of observing the shape after cutting, the approximate impact range of the object after cutting is determined, and the vertices that have basically no impact are automatically marked by setting the distance from the cutting surface, such as Fig. 9 As shown in (a), the red dots are the vertices that cannot be moved.

[0091] In order to apply appropriate external force to each vertex, the weighted normal of each vertex is calculated and the appropriate external force is applied along the normal. Fig. 9 (b) shows that the blue line is the normal of each vertex.

[0092] Finally, in order to solve the spring model, the following formula was derived:

[0093] 1. According to Hooke's law, we know that vertex x is in a certain period of time. The moving position of is as follows:

[0094]

[0095]

[0096]

[0097] in is the initial side length, is the mass of the current vertex, k is the spring coefficient, is the current side length, Vertex By Location The next position to move to;

[0098] 2. From the formula derivation, it can be seen that the implicit equation in 1 is equivalent to solving the minimum value of the following equation:

[0099]

[0100]

[0101]

[0102] in for The minimum value of

[0103] 3. In 2, solve the target position Transform to solve The minimum value problem of , and then the following inference:

[0104]

[0105]

[0106]

[0107]

[0108] Because of the solution The minimum value of , so its derivative should be 0, and then The position is expanded by Taylor's formula, and the above formula is inferred, where H is the Hessian matrix. Because the final solution equation is a quadratic convex function, the Newton iteration method is used to solve it, such as Fig.13 The results obtained after simulation.

[0109] Step 7: Perform Delaunay re-tetrahedronization and BSP subdivision on the local area to generate a new BSP tree B, and merge BSP trees A and B.

[0110] After Delaunay re-tetrahedronization and BSP subdivision of the local area, the following is generated: Fig.10 In the local area shown, each tetrahedron still meets the definition of a bsp node, and each face of the tetrahedron can be regarded as a plane, so a bsp tree B can be constructed. The bsp tree A and the bsp tree B are based on the adjacency relationship of the original splicing faces, so the bsp tree A and the bsp tree B can be simply merged.

[0111] Step 8: Traverse all BSP nodes and use the center of gravity to re-tetrahedralize the non-tetrahedral nodes.

[0112] Because in the process of re-tetrahedron, there are only two problems that need to be solved, namely:

[0113] 1. The original splicing surface becomes a one-to-many relationship such as Fig.11 As shown in (a).

[0114] 2. The original one becomes a one-to-many relationship. Fig.12 As shown in (b).

[0115] Fig.11 The original splicing plane BCD shown in (a) is split into two planes BED and ECD, so tetrahedron ABCD needs to be re-tetrahedralized, where tetrahedrons FEBD and FEDC belong to the nodes of the original bsp tree B, and tetrahedron ABCD belongs to the node of bsp tree A. Point G is obtained by calculating the centroid of tetrahedron ABCD, and then it is re-tetrahedralized based on the centroid, as shown in Fig.11 In (b), ABCD can be subdivided into the following tetrahedrons: ABCG, ABDG, BEDG, and ECDG.

[0116] Fig.12 The original splicing surface ABC shown in (a) is re-subdivided into two triangular surfaces ABE and AEC due to the local cutting area. Fig.12 In (b), the adjacent faces of the original CB edge become the two ends of BE and EC, and the original triangle BCD becomes a polygon BDEC. Therefore, the polygon BDEC is re-triangulated based on the centroid, as shown in Fig.12 In (c), the center of gravity is calculated to obtain point G, and the original BDEC is re-divided into triangles BEG, EGC, GDC, and BDG. In this way, the center of gravity-based re-tetrahedronization can be performed.

[0117] Step 9, optimizing the tetrahedron model after simulated cutting;

[0118] First, construct an aabb tree structure for the surface, and construct a minimum heap for all edges based on the edge length. Extract the edge with the smallest edge length, and merge it if the edge length is less than the minimum threshold set by the user. If the two nodes V1 and V2 of the edge are not surface nodes, replace all tetrahedrons containing V2 nodes with V1 nodes, and translate the position of V1 nodes to (V1+V2) / 2. If the edge contains a surface node, use the aabb tree to determine whether the distance from (V1+V2) / 2 to the surface is within the deviation set by the user. If it is, merge it. When there is no edge in the heap, the edge length optimization algorithm ends. Execute the Laplace smoothing method based on quadrilateral mesh for all nodes. The final effect is as follows Fig.14 As shown in (a) and (b).

[0119] In summary, the present invention samples the position of the cutting surface based on the time step, constructs a suitable bounding box for the sampling and extracts a suitable area sufficient for physical simulation, and divides the surface of the area into simulation surfaces and splicing surfaces, and forms a simulation surface that can be physically simulated for the simulation surface and the sampled cutting surface, which makes the generation of the cutting surface 100% accurate. A spring model is constructed, a physical simulation is performed, and then it is re-tetrahedronized and BSP subdivided and post-processed, which effectively reduces the generation of pathological units. The virtual surgery simulation method based on grid adaptive cutting and subdivision proposed in this discovery ensures better accuracy and quality of the cut objects.

[0120] The embodiments of the present invention are described in detail above with reference to the accompanying drawings, but the present invention is not limited to the described embodiments. For those skilled in the art, various changes, modifications, substitutions and variations of these embodiments including components are made without departing from the principles and spirit of the present invention, and still fall within the scope of protection of the present invention.

Claims

1. A virtual surgery simulation method based on grid adaptive cutting and subdivision, characterized in that: The steps include: Step 1: tetrahedralize the input object to be cut, and record the position of the scalpel once every time step. Its starting position and the position after moving one time step together constitute the cutting surface on the tetrahedron; Step 2: According to the cutting surface and the maximum side length of the tetrahedron, construct the bounding box and the bsp tree A, mark the bsp nodes involved in the bounding box, and mark the splicing surface and the simulation surface; Step 3: Generate common cutting lines based on the simulated surface and the cutting surface by using the constrained Delaunay triangulation method, pre-process the common cutting lines and cutting surface, and obtain the simulated cutting surface; Step 4: Mark the movable vertices and the immovable vertices respectively, then use Newton's iteration method according to Hooke's law to solve the position of the new vertex, move the vertex, and simulate the shape after cutting; Step 5: Re-tetrahedralize and BSP subdivide the local area based on the simulated shape, generate a new BSP tree B, and merge BSP trees A and B; Step 6: After merging, traverse all BSP nodes, and for non-tetrahedral nodes, use the center of gravity of the non-tetrahedron to tetrahedronize them, complete the simulated cutting, and perform optimization.

2. The virtual surgery simulation method based on grid adaptive cutting and subdivision according to claim 1, characterized in that: The specific process of forming the cutting surface on the tetrahedron in step 1 is: obtain the four endpoints A, B, C, and D that constitute the starting line segment AB and the ending line segment CD, and use the four endpoints to generate a cutting plane composed of two triangles, where the triangles are ABC and ACD.

3. The virtual surgery simulation method based on grid adaptive cutting and subdivision according to claim 2, characterized in that: The specific implementation process of constructing the bounding box and the bsp tree A is as follows: Step 2.1.

1. By traversing the edge lengths of all tetrahedrons, get the maximum edge length len_ A ; Traverse the edge lengths on the cutting plane and get the maximum edge length len_ of the cutting plane B ; Take the maximum edge length of the tetrahedron len_ A and the maximum side length of the cutting plane len_ B The maximum value in len_ C , and then give len_ C Multiply by a coefficient as the offset value of the diagonal of the bounding box; Step 2.1.2, according to the positions of the four vertices of the cutting plane ABCD, calculate the maximum value of these four points in the three-dimensional coordinate axis (x max ,y max , z max ) and minimum value (x min ,y min , z min ); Step 2.1.3, the bounding box is composed of 8 vertices, and its diagonal is (x max + offset,y max + offset, z max + offset) and (x min - offset, y min -offset,z min -offset), so the 8 vertices are (x max +offset,y max + offset, z max + offset), (x max + offset,y min - offset, z max + offset), (x max +offset,y max + offset, z min - offset), (x max + offset,y min - offset, z min - offset)、(x min -offset,y max + offset, z max + offset)、(x min - offset, y max + offset, z min - offset)、(x min -offset,y min - offset, z max + offset)、(x min - offset, y min - offset, z min - offset); Step 2.1.4, traverse the cut tetrahedron, where each triangular surface is regarded as a plane, and each tetrahedron is regarded as a leaf node of a bsp unit, and a bsp tree A is created.

4. The virtual surgery simulation method based on grid adaptive cutting and subdivision according to claim 3 is characterized in that: The BSP nodes involved in the bounding box are marked, and the splicing surface and the simulated surface are marked. The specific implementation process is as follows: Step 2.2.1, use Jacobson algorithm to construct the inner-outer function of the bounding box; Step 2.2.2, traverse all BSP nodes and calculate the centroid of each node; Step 2.2.3, use the inner-outer function to determine whether the centroid is contained in the bounding box. If it is in the bounding box, mark the bsp node. If it is not in the bounding box, do not process it. Step 2.2.4, traverse all marked BSP nodes, if the face in the node is adjacent to an unmarked BSP node, and the face is not the surface of the cut object, then mark the face as a splicing face; Step 2.2.5, traverse all the faces in the marked bsp nodes. If the face is the surface of the cutting object, mark the face as a simulation surface.

5. The virtual surgery simulation method based on grid adaptive cutting and subdivision according to claim 4, characterized in that: The specific implementation process of step 3 is as follows: Step 3.1, traverse the bsp nodes and extract all simulation surfaces; Step 3.2, use the simulation surface to intersect with the cutting surface to generate constrained edges and points, and then perform constrained Delaunay triangulation to generate common cutting lines; Step 3.3, find the starting endpoint of the common cutting line, and use this endpoint as the initial vertex A1, and the vertex on the next common cutting line connected to it is called B1; randomly divide the two adjacent faces of the A1B1 edge into two categories, and generate point B1' as the symmetric vertex of B1, then use the face in the first category as the starting face, rotate around point B1, and all the faces it touches are classified as the first category until it touches the common cutting line B1C1, and change all B1 points in the first category to B1', and the rest are the second category without change, and so on until it touches the end point of the cutting line; Step 3.4, splice the common cutting line and the cutting surface to form a simulated cutting surface: copy the points on the cutting surface except the points on the boundary and not on the common cutting line, and form a symmetric surface of the new cutting surface with the newly generated common cutting line and the points on the boundary, and give it a label of the simulated cutting surface.

6. The virtual surgery simulation method based on grid adaptive cutting and subdivision according to claim 5, characterized in that: The specific implementation process of step 5 is as follows: Step 5.1, extract all nodes contained in the simulation surface and the spliced ​​surface; Step 5.2, delete all bsp nodes in the original bsp tree A that are included in the local area; Step 5.3, use Delaunay to re-tetrahedralize all nodes; Step 5.4, construct a new bsp tree structure B for the tetrahedron and its contained surfaces; Step 5.5, using the BSP subdivision method, restore the re-tetrahedralized feature surface; Step 5.6: Merge bsp tree A and bsp tree B, and replace the original spliced ​​surface with all new sub-surfaces.

7. The virtual surgery simulation method based on grid adaptive cutting and subdivision according to claim 6, characterized in that: In step 6, the center of gravity is calculated by taking the average value of all node positions contained in the non-tetrahedron.

8. The virtual surgery simulation method based on grid adaptive cutting and subdivision according to claim 7, characterized in that: In step 6, the specific implementation process of the optimization is as follows: Step 6.1, constructing an aabb tree structure on the tetrahedron surface; Step 6.2: construct a minimum heap based on the edge lengths of all tetrahedrons, extract the edge with the smallest edge length, and merge it if the edge length is less than the set minimum threshold; Step 6.3, if the two nodes V1 and V2 of the edge are not surface nodes, replace all tetrahedrons containing V2 nodes with V1 nodes, and translate the position of V1 nodes to (V1+V2) / 2; Step 6.4, if the edge contains a surface node, use the aabb tree to determine whether the distance from (V1+V2) / 2 to the surface is within the deviation set by the user. If it is, merge it. Otherwise, extract the next edge with the smallest length and repeat steps 6.3 and 6.

4. Step 6.5: When there are no edges in the heap, the edge length optimization ends; Step 6.6: Execute the Laplace smoothing method based on quadrilateral mesh on all nodes to complete node optimization.

Images