An object deformation display method and device based on force feedback interaction

Abstract

This invention discloses a method and apparatus for displaying object deformation based on force feedback interaction. When CT data of the object to be interacted with is received, a polygonal bounding volume is constructed using the CT data. The polygonal bounding volume is then tetrahedralized to generate a tetrahedral mesh model. Based on the vertex relationships between tetrahedral units and the unit volume of each tetrahedral unit, the equivalent mass of the mesh particles in each tetrahedral unit is allocated to obtain a soft mesh model. External forces are applied to the soft mesh model in response to control commands. Based on the equivalent mass of the mesh particles and the force corresponding to the applied external force, the particle positions of each tetrahedral unit are iteratively updated to obtain the updated particle positions. The image is rendered based on the updated particle positions, outputting the deformation image corresponding to the object to be interacted with. By combining tetrahedral modeling and iterative updates of particle positions, the accuracy and consistency of force feedback control for interaction between external forces and soft objects are effectively improved.

Description

Technical Field

[0001] This invention relates to the field of robot control technology, and in particular to a method and apparatus for displaying object deformation based on force feedback interaction. Background Technology

[0002] With the development of intelligent simulation and medical robot technology, traditional physical modeling and manual control training modes have pain points such as difficulty in adjusting model parameters, low dynamic characteristic reproducibility, high control practice cost and lack of quantitative feedback, which can no longer meet the needs of high-precision and high-realism intelligent interactive training.

[0003] To address this issue, high-fidelity software modeling can achieve precise and real-time adaptation to robot control. Among these algorithms, Extended Position-Based Dynamics (XPBD) has become one of the mainstream algorithms for real-time software modeling due to its strong real-time performance, high stability, and support for large deformation simulation.

[0004] However, in current solutions that integrate software modeling and robot control, traditional algorithms can only simulate pure elastic deformation, and there is a lack of viscoelastic modeling of soft materials. This results in an inability to accurately reproduce the dynamic mechanical properties of the soft material, leading to insufficient accuracy and consistency in force feedback control of robot-soft material interaction. Summary of the Invention

[0005] This invention provides a method and device for displaying object deformation based on force feedback interaction, which solves the technical problem that traditional algorithms can only simulate pure elastic deformation, lack viscoelastic modeling of soft bodies, and cannot accurately reproduce the dynamic mechanical properties of soft bodies, resulting in insufficient accuracy and consistency of force feedback control in robot-soft body interaction.

[0006] This invention provides a method for displaying object deformation based on force feedback interaction, comprising:

[0007] When the CT data of the object to be interacted with is received, a polygon bounding volume is constructed using the CT data;

[0008] The polygonal bounding volume is tetrahedralized to generate a tetrahedral mesh model; the tetrahedral mesh model is composed of multiple tetrahedral elements;

[0009] Based on the vertex relationships between the tetrahedral elements and the element volume of each tetrahedral element, the equivalent mass of the mesh particles of each tetrahedral element is allocated to obtain the soft mesh model.

[0010] The system responds to control commands by applying external forces to the soft mesh model.

[0011] Based on the equivalent mass of the mesh particles and the force corresponding to the applied external force, the particle position of each tetrahedral unit is iteratively updated to obtain the updated particle position.

[0012] Based on the updated positions of each particle, the screen is rendered, and the deformed screen corresponding to the object to be interacted with is output.

[0013] Optionally, when receiving CT data of the object to be interacted with, constructing a polygon bounding volume using the CT data includes:

[0014] When the CT data of the object to be interacted with is received, the initial patch model corresponding to the CT data is obtained;

[0015] Match the corresponding voxel mesh step size according to the scale of the initial patch model, and construct a voxel mesh based on the voxel mesh step size; the voxel mesh includes multiple mesh elements;

[0016] Select mesh cells whose center position is within the initial patch model to construct an occupied cell set;

[0017] The initial bounding polygon is obtained by extracting and smoothing the outer contour of the occupied unit set;

[0018] The initial bounding polygon is offset outward along the normal direction of each of the mesh cells by the outward offset corresponding to the object to be interacted with, thereby generating a polygon bounding volume.

[0019] Optionally, the step of tetrahedralizing the polygonal bounding volume to generate a tetrahedral mesh model includes:

[0020] Construct a single initial tetrahedron covering the bounding volume of the polygon;

[0021] Each vertex of the polygonal bounding body is inserted into the initial tetrahedron in sequence, and all extra tetrahedrons containing the insertion points on the circumscribed spheres are deleted to form a common surface cavity;

[0022] Connect the insertion point with all the boundary vertices of the common surface cavity to obtain the intermediate tetrahedron;

[0023] After all the bounding volume vertices have been inserted, the initial tetrahedron is removed to obtain a tetrahedral mesh model.

[0024] Optionally, the step of allocating the equivalent mass of mesh particles of each tetrahedral element based on the vertex association relationship between each tetrahedral element and the element volume of each tetrahedral element to obtain a soft mesh model includes:

[0025] Calculate the unit volume of each tetrahedral unit according to the coordinates of its four vertices;

[0026] If any mesh particle is associated with only a single tetrahedral unit, then the equivalent mass of the mesh particle is calculated according to the first equivalent mass calculation formula and allocated to the mesh particle;

[0027] If any mesh particle is associated with at least two of the tetrahedral elements, the second equivalent mass calculation formula calculates the equivalent mass of the mesh particle and assigns it to the mesh particle;

[0028] Once all the mesh particles are assigned equivalent mass, a soft mesh model is obtained.

[0029] Optionally, it also includes:

[0030] After obtaining the soft mesh model, all tetrahedral elements within the soft mesh model are traversed, and the nearest neighbor distance probability density distribution index and dihedral angle quality index are calculated one by one.

[0031] If there are unqualified tetrahedral units that do not meet the preset standards for the nearest neighbor distance probability density distribution index and / or the dihedral angle quality index, and the number of unqualified tetrahedral units is less than the preset unit threshold, then the region to be adjusted where the unqualified tetrahedral unit is located is located, and the polygon bounding volume corresponding to the region to be adjusted is updated according to the new voxel mesh step size.

[0032] Jump to the step of performing tetrahedral subdivision of the polygon bounding volume to generate a tetrahedral mesh model;

[0033] If there are unqualified tetrahedral cells that do not meet the preset standards for the nearest neighbor distance probability density distribution index and / or the dihedral angle quality index, and the number of unqualified tetrahedral cells is not less than the preset cell threshold, then a new voxel mesh step size and a new outward offset are obtained, and the process jumps to the step of constructing a polygon bounding volume using the CT data.

[0034] If the nearest neighbor distance probability density distribution index and dihedral angle quality index of all tetrahedral elements meet the preset compliance conditions, then the response control command is executed to apply external force to the soft mesh model.

[0035] Optionally, the step of iteratively updating the particle position of each tetrahedral unit based on the equivalent mass of the mesh particle and the force corresponding to the applied external force to obtain the updated particle position includes:

[0036] The current velocity of each particle in the tetrahedral unit is updated based on the force corresponding to the external force application action, and the dynamic velocity dissipation coefficient is calculated in combination with the duration of soft deformation.

[0037] The initial velocity of each vertex is attenuated and corrected using the dynamic velocity dissipation coefficients described above to obtain the particle corrected velocity.

[0038] The initial predicted position of each particle in the current frame is predicted based on the particle correction velocity.

[0039] At each simulation time step, a dynamic compliance constant is calculated based on the duration of soft deformation, and the Lagrange operator of the length constraint is adjusted based on the dynamic compliance constant to obtain the updated length constraint.

[0040] The initial predicted position is iteratively corrected under the updated length constraint, mention constraint, and collision penetration prevention constraint to obtain the updated particle position and update the current velocity of the particle.

[0041] Optionally, the formula for calculating the dynamic velocity dissipation coefficient is:

[0042] ;

[0043] in, The dynamic velocity dissipation coefficient, The velocity dissipation time-varying coefficient, t is the preset base dissipation coefficient, and t is the duration of the soft deformation.

[0044] Optionally, the formula for calculating the dynamic compliance constant is:

[0045] ;

[0046] Where k, b, and a are preset KV model constants. Corrects the particle velocity.

[0047] Optionally, the step of rendering the image based on the updated positions of each particle and outputting the deformed image corresponding to the object to be interacted with includes:

[0048] The position of each particle is updated and transformed to screen pixel coordinates;

[0049] The triangular facet formed by the vertices to which the particle update position belongs is rasterized to obtain multiple pixels and associated with the screen pixel coordinates;

[0050] Each pixel is assigned brightness and / or color to generate and output the deformed image corresponding to the object to be interacted with.

[0051] The present invention also provides an object deformation display device based on force feedback interaction, comprising:

[0052] A polygon bounding volume construction module is used to construct a polygon bounding volume using the CT data of the object to be interacted with when the CT data of the object is received.

[0053] The tetrahedral subdivision module is used to perform tetrahedral subdivision on the polygonal bounding volume to generate a tetrahedral mesh model; the tetrahedral mesh model is composed of multiple tetrahedral elements;

[0054] The particle mass allocation module is used to allocate the equivalent mass of the mesh particles of each tetrahedral unit according to the vertex association relationship between each tetrahedral unit and the unit volume of each tetrahedral unit, so as to obtain a soft mesh model.

[0055] An external force application module is used to perform external force application actions on the soft mesh model in response to control commands;

[0056] The position update module is used to iteratively update the particle position of each tetrahedral unit based on the equivalent mass of the mesh particle and the force corresponding to the action of the external force applied, so as to obtain the particle update position;

[0057] The image rendering module is used to render the image based on the updated position of each particle and output the deformation image corresponding to the object to be interacted with.

[0058] As can be seen from the above technical solutions, the present invention has the following advantages:

[0059] Upon receiving CT data of the object to be interacted with, a polygonal bounding volume is constructed using the CT data. This bounding volume is then tetrahedralized to generate a tetrahedral mesh model. Based on the vertex relationships between tetrahedral elements and the element volume of each tetrahedral element, the equivalent mass of the mesh particles for each tetrahedral element is allocated, resulting in a soft mesh model. External forces are applied to the soft mesh model in response to control commands. Based on the equivalent mass of the mesh particles and the forces corresponding to the applied external forces, the particle positions of each tetrahedral element are iteratively updated to obtain the updated particle positions. Finally, the scene is rendered based on the updated particle positions, outputting the deformation display corresponding to the object to be interacted with. This combination of tetrahedral modeling and iterative particle position updates effectively improves the accuracy and consistency of force feedback control for external forces and soft-body interaction. Attached Figure Description

[0060] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0061] Figure 1A flowchart illustrating the steps of a method for displaying object deformation based on force feedback interaction, as provided in an embodiment of the present invention;

[0062] Figure 2 This is a schematic diagram illustrating the generation process of a polygonal bounding volume according to an embodiment of the present invention;

[0063] Figure 3 A schematic diagram illustrating the process of generating a polygonal bounding volume for the soft organs liver and gallbladder, provided in an embodiment of the present invention;

[0064] Figure 4 This is a schematic diagram illustrating the process of generating a polygonal bounding body for the gallbladder, as provided in an embodiment of the present invention.

[0065] Figure 5 A schematic diagram illustrating the generation process of the tetrahedral mesh model provided in an embodiment of the present invention;

[0066] Figure 6 A schematic diagram illustrating the construction of a tetrahedral mesh model of the liver and gallbladder models provided in an embodiment of the present invention;

[0067] Figure 7 This is a detailed schematic diagram of a tetrahedral mesh model provided in an embodiment of the present invention;

[0068] Figure 8 A schematic diagram of particle mass distribution for a tetrahedral mesh model provided in an embodiment of the present invention;

[0069] Figure 9 A logical flowchart for model evaluation provided in this embodiment of the invention;

[0070] Figure 10 This is a schematic diagram showing the comparison of indicators of a liver model tetrahedron provided in an embodiment of the present invention;

[0071] Figure 11 This is a schematic diagram of the nearest neighbor distance probability density distribution provided in an embodiment of the present invention;

[0072] Figure 12 A bar chart comparing dihedral angle indices provided in embodiments of the present invention;

[0073] Figure 13 This is a schematic diagram of the damped oscillation phenomenon of a soft mesh model after a force is applied, provided in an embodiment of the present invention.

[0074] Figure 14 A detailed schematic diagram of a damped oscillation phenomenon provided in an embodiment of the present invention;

[0075] Figure 15 A schematic diagram of an experimental platform provided in an embodiment of the present invention;

[0076] Figure 16A schematic diagram of an experimental process provided in an embodiment of the present invention;

[0077] Figure 17 This is a schematic diagram of experimental image and marker recognition provided in an embodiment of the present invention;

[0078] Figure 18 A location trajectory recording curve is provided as an embodiment of the present invention;

[0079] Figure 19 A schematic diagram of key frame and key point changes in the simulation process provided in this embodiment of the invention;

[0080] Figure 20 A schematic diagram of the trajectory of the simulation procedure provided in the embodiments of the present invention;

[0081] Figure 21 This is a structural block diagram of an object deformation display device based on force feedback interaction, provided in an embodiment of the present invention. Detailed Implementation

[0082] This invention provides a method and apparatus for displaying object deformation based on force feedback interaction, which solves the technical problem that traditional algorithms can only simulate pure elastic deformation, lack viscoelastic modeling of soft bodies, and cannot accurately reproduce the dynamic mechanical properties of soft bodies, resulting in insufficient accuracy and consistency of force feedback control in robot-soft body interaction.

[0083] To make the objectives, features, and advantages of this invention more apparent and understandable, the technical solutions of the embodiments of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the embodiments described below are only some embodiments of this invention, and not all embodiments. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0084] Please see Figure 1 , Figure 1 This is a flowchart illustrating the steps of a method for displaying object deformation based on force feedback interaction, as provided in an embodiment of the present invention.

[0085] This invention provides a method for displaying object deformation based on force feedback interaction, comprising the following steps:

[0086] Step 101: When the CT data of the object to be interacted with is received, a polygon bounding volume is constructed using the CT data;

[0087] The object to be interacted with refers to an object with viscoelastic properties, such as animal soft tissue, including but not limited to liver, gallbladder, trachea, and lumbar vertebral nucleus pulposus.

[0088] CT data refers to continuous tomographic images of an interactive object obtained through computed tomography (CT) scans, containing core geometric and physical information such as the object's outline and tissue density distribution.

[0089] A polygon bounding volume refers to a triangular mesh model that is watertight and closed. It is the geometric basis for subsequent tetrahedral subdivision and can completely restore the macroscopic geometric shape of the object to be interacted with.

[0090] In this embodiment, after receiving CT tomographic image data with interactive objects, the contour boundary of the interactive objects is extracted using a medical image segmentation algorithm to generate an initial triangular patch model. Based on the scale of the initial model, the voxel mesh step size is determined and a regular voxel mesh covering the entire model is constructed. All voxel mesh units are traversed to determine their spatial position relationship with the initial model. The set of occupied units located inside the model is selected, the outer surface contour of the occupied unit set is extracted and smoothed, and finally the smoothed contour is uniformly expanded and offset along the patch normal to generate a polygonal bounding body with watertightness.

[0091] Furthermore, during the construction of the polygon bounding volume, tissue density information from CT data is simultaneously incorporated, and different voxel mesh step sizes are set for tissue regions of different densities. Smaller mesh step sizes are used for high-density regions containing fine structures such as blood vessels and bile ducts to preserve geometric details, while relatively larger mesh step sizes are used for homogeneous soft tissue regions to control computational load. At the same time, in the normal outward offset stage, different offset amounts are set for the mechanical properties of different tissues to avoid numerical penetration problems in subsequent collision simulations.

[0092] In one example of the present invention, step 101 may include the following sub-steps:

[0093] When the CT data of the object to be interacted with is received, the initial patch model corresponding to the CT data is obtained;

[0094] Match the corresponding voxel mesh step size according to the scale of the initial patch model, and construct the voxel mesh based on the voxel mesh step size; the voxel mesh consists of multiple mesh elements;

[0095] Select the mesh cells whose center position is within the initial patch model to construct the occupied cell set;

[0096] Extract and smooth the outer contour of the occupied cell set to obtain the initial bounding polygon;

[0097] Offset the initial bounding polygon outwards along the normal direction of each grid cell by the corresponding outward offset of the object to be interacted with, and generate the polygon bounding volume.

[0098] The initial patch model refers to the triangular patch mesh model generated based on CT data through medical image segmentation and surface reconstruction algorithms, which is used to restore the original geometric shape of the object to be interacted with.

[0099] The voxel grid step size refers to the side length of a single cubic grid cell in a voxel grid.

[0100] A voxel mesh refers to a three-dimensional regular cubic mesh that covers the entire range of the initial patch model, constructed based on the voxel mesh step size, and is composed of multiple independent mesh units.

[0101] The occupancy set refers to the collection of all mesh cells whose center position is inside the initial patch model, used to characterize the extent to which the object to be interacted occupies in three-dimensional space.

[0102] The initial bounding polygon refers to the closed triangular mesh obtained by extracting the outer contour of the occupied cell set and smoothing it.

[0103] Outward offset refers to the distance parameter by which the initial enclosing polygon is offset outward along the normal of the face. It is used to reserve a safety margin for subsequent collision simulation and avoid numerical penetration problems.

[0104] In this embodiment, after receiving the CT data of the object to be interacted with, the anatomical contour of the object to be interacted with is first extracted from the CT data using a medical image segmentation algorithm, and a corresponding initial patch model, such as a mesh set, is generated using a surface reconstruction algorithm. Then, based on the overall bounding box size and minimum anatomical structure scale of the initial patch model, the corresponding voxel mesh step size h is matched, and a three-dimensional regular voxel mesh that completely covers the initial patch model is constructed based on this voxel mesh step size. The voxel mesh is composed of multiple equally sized cubic mesh cells. The center position vector of each mesh cell is... for:

[0105] ;

[0106] in, This represents the smallest corner point of the grid cell (the smallest i, j, k point in the cell).

[0107] Please see Figure 2 , Figure 2 This is a schematic diagram illustrating the generation process of a polygonal bounding volume in an embodiment of the present invention. By traversing all mesh cells in the voxel mesh, iterating through each mesh cell, determining whether the center position of the mesh cell is inside the original mesh M, and distinguishing between all inner and outer mesh cells, an occupancy set G is constructed based on the mesh cells within the shape. vSpecifically, a spatial location judgment algorithm is used to verify whether the center position of each mesh cell is within the internal space of the initial patch model. All mesh cells whose center positions are within the initial patch model are then selected, and an occupied cell set G is constructed. v ,like Figure 2 The yellow area.

[0108] Then, extract the outer surface triangular facets of the occupying cell set to form the initial contour. A smoothing algorithm is then used to smooth the initial contour to eliminate the jagged edges caused by voxelization, resulting in the initial bounding polygon as shown below. Figure 2 The area enclosed by the red line is then further defined. Finally, based on the organization type of the object to be interacted with and the subsequent collision simulation requirements, the corresponding outward offset is determined. All triangular faces of the initial enclosing polygon are then offset outward along their normal directions by this outward offset, generating a polygonal enclosing body with watertightness and closure, as shown below. Figure 2 The blue line encloses the region. The above processing provides a foundation for the subsequent I-XPBD framework tetrahedron in several aspects. First, the number and complexity of the mesh can be controlled by the step size h, avoiding uncontrollable computational costs in subsequent iterations caused by inconsistent accuracy of soft organ source data acquired from different formats. Second, many models acquired from CT data exhibit surface breaks; this processing ensures the closure of the enclosing polygons, generating a watertight mesh. Third, the generated enclosing polygon patches have relatively consistent lengths, which is beneficial for establishing stable voxelized topological relationships. If the subsequently established I-XPBD particle units are uniformly distributed in geometric space, the mass of each particle can be directly determined from the spatial density distribution. Finally, the use of a normal outward expansion method avoids numerical intrusion in subsequent interactive simulation processes such as collision detection.

[0109] In practical implementation, Blender open-source software can be used as the main preprocessing software, and the function of generating the watertight polygon bounding volume described above can be implemented by writing a Python plugin. To test the effectiveness of the algorithm, the liver, a typical soft organ, and the gallbladder, a smaller-scale organ, were selected for testing. Figure 3 As shown, Figure 3 (a) in the image represents the original shape of the liver. Figure 3 In (b), the step size h=1 and the outward offset L are... o =1 generates a polygon bounding volume (vertices = 508, faces = 506). Figure 3 In the diagram, (c) represents the step size h = 0.5 and the outward offset L. o =1 generates a polygon bounding volume (vertices = 2050, faces = 2048). Figure 3 In this context, (e) represents the step size h = 0.1 and the outward offset L. o =1 generates a polygon bounding volume (vertices = 51556, faces = 51554). Figure 3 In this context, (f) represents the step size h = 0.5 and the outward offset L. o The polygon bounding volume generated by setting the value to 0.5 (vertices = 2050, faces = 2048). Figure 3 (d) in Figure 3 (a) and Figure 3 Comparing cases (b) in the above examples, it is evident that when the step size h is too large, local details of the original shape are lost. Therefore, if the hardware and software support this, skinning or other methods should be used to improve the detail representation of the software model during simulation. Figure 3 The results in (b), (c), and (d) show that the number of vertices, faces, and quality of the generated polygon bounding volume can be well controlled by the step size h, achieving its functional requirements. (Comparison) Figure 3 (c) and Figure 3 In (f), the outward offset L can be seen. o It affects the size of the polygon bounding volume, but does not directly affect the number of vertices or faces, thus satisfying the requirement of single step size control.

[0110] In addition, to more intuitively demonstrate the modeling effect, a gallbladder model with a more complex shape and relatively smaller scale was selected for testing, and the results are as follows: Figure 4 As shown. Figure 4 (a) in the image represents the original shape of the gallbladder. Figure 4 In (b), the step size h = 0.2 and the outward offset L are... o =0.2 generated polygon bounding volume (vertices = 836, faces = 804). Figure 4 In this context, (c) represents the step size h = 0.1 and the outward offset L. o =0.2 generated polygon bounding volume (vertices = 3336, faces = 3328). Figure 4 In this context, (d) represents the step size h = 0.05 and the outward offset L. o The polygon bounding volume generated by 0.2 (vertices = 13514, faces = 13512). Figure 4 As shown in (b), excessively large step sizes result in disruption of connectivity. For more complex soft organ models, such as those with blood vessels, bile ducts, or other adhesions, it is crucial to ensure that the step size h is less than the minimum thickness of these critical structures. (Comparison) Figure 4 As can be seen from (a), (c), and (d) in the figure, for soft organ shapes with high complexity and small scale, setting the step size h appropriately can also effectively control the number of vertices, the number of faces, and the quality of the polygon bounding body, generating a polygon bounding body that meets the watertightness requirements.

[0111] Step 102: Tetrahedral subdivision is performed on the polygon bounding volume to generate a tetrahedral mesh model; the tetrahedral mesh model is composed of multiple tetrahedral elements;

[0112] A tetrahedral element is the smallest tetrahedral geometry in three-dimensional space, consisting of 4 non-coplanar vertices, 6 edges, and 4 triangular facets. It is the smallest computational unit for software physics simulation.

[0113] A tetrahedral mesh model refers to a voxelized mesh model formed by splicing together multiple tetrahedral elements by sharing vertices, edges, or faces.

[0114] In this embodiment, all vertices of the generated watertight polygon bounding volume are used as the input point set. The Bowyer-Watson algorithm is used to perform 3D Delaunay tetrahedral subdivision. First, an initial bounding tetrahedron that completely covers all input point sets is constructed. Then, the vertices in the input point set are inserted into the current mesh one by one. Bad tetrahedrons containing the insertion points in the circumscribed sphere are deleted to form cavities. All boundary vertices of the insertion points and cavities are connected to reconstruct tetrahedral elements. After all vertices are inserted, the elements related to the initial bounding tetrahedron are deleted. Finally, a tetrahedral mesh model composed of multiple tetrahedral elements is generated, and a topological association index between each tetrahedral element and a vertex is established.

[0115] It should be noted that during the tetrahedral meshing process, for areas prone to degenerate elements, such as slender branch structures and thin soft tissue walls of the object to be interacted with, adaptive densification meshing and real-time mesh quality verification are performed. The dihedral angle parameters are calculated synchronously for each tetrahedral element generated by the meshing. Degenerate tetrahedral elements with dihedral angles exceeding the reasonable range are removed and re-meshed. At the same time, the size differences between adjacent tetrahedral elements are smoothed to avoid numerical oscillations in subsequent simulations caused by abrupt changes in mesh size.

[0116] In one example of the present invention, step 102 may include the following sub-steps:

[0117] Construct a single initial tetrahedron that covers the bounding volume of the polygon;

[0118] Insert each vertex of the bounding volume of the polygon into the initial tetrahedron in turn, and delete all extra tetrahedrons containing the insertion points of the circumscribed spheres to form a common surface cavity;

[0119] Connect the insertion point to all boundary vertices of the common surface cavity to obtain the intermediate tetrahedron;

[0120] After all the bounding volume vertices have been inserted, the initial tetrahedron is removed to obtain the tetrahedral mesh model.

[0121] In this embodiment, a single initial tetrahedron, whose spatial extent completely covers all vertices of the bounding volume of the polygon, is first constructed to establish the initial mesh for 3D Delaunay meshing. Then, each vertex of the bounding volume of the polygon is sequentially inserted into the current tetrahedral mesh as an insertion point. For each insertion point, all tetrahedral elements within the current mesh are traversed, and it is verified whether the circumsphere of each tetrahedral element contains the insertion point. Any extra tetrahedrons whose circumspheres contain the insertion point are deleted, and the common faces of the remaining valid tetrahedral elements enclose a closed common-face cavity. The insertion point is then connected to all boundary vertices of the common-face cavity one by one, generating multiple intermediate tetrahedrons conforming to the Delaunay criterion within the cavity, completing the mesh update after a single vertex insertion. This process of vertex insertion, invalid element removal, cavity reconstruction, and mesh update is repeated until all vertices of the bounding volume of the polygon are inserted. Finally, the initial tetrahedron and all redundant elements associated with it are removed, resulting in a tetrahedral mesh model that perfectly fits the geometry of the polygon and consists of multiple tetrahedral elements conforming to the Delaunay criterion. Figure 5 As shown, Figure 5 This is a schematic diagram illustrating the generation process of the tetrahedral mesh model in an embodiment of the present invention. First, an initial enclosing triangle (corresponding to an initial tetrahedron in 3D) completely covering all points to be inserted is constructed as the starting carrier for meshing. Then, each vertex to be inserted (e.g., p1, p2…p…) is sequentially… n Insert the element into the current mesh. Each time you insert, first delete all invalid triangles (tetrahedrons) containing the insertion point from the circumcircle (circumsphere in 3D). This forms a cavity enclosed by the common edges (common faces in 3D) of the remaining elements. Then, connect the insertion point to all boundary vertices of the cavity to generate new elements that conform to the Delaunay criterion. This completes a single iteration update. Repeat the above insertion, deletion, and reconstruction process until all vertices have been inserted. Finally, remove the initial enclosing triangle (initial tetrahedron in 3D) and associated redundant elements to obtain the final Delaunay triangulation mesh (tetrahedral mesh model in 3D).

[0122] It should be noted that the tetrahedron in this embodiment can be the Delaunay tetrahedron as the basic configuration of the algorithm. The tetrahedron mesh has better quality and higher stability, avoiding tetrahedrons with "too small angles and too long sides", ensuring that the constraint gradient is within a reasonable range and preventing the values ​​from being too extreme or even overflowing. The particle distribution is more uniform, which is beneficial for the subsequent calculation of the mass weighting coefficients of discrete particles based on the density field. The internal force distribution is more reasonable, so as to avoid force oscillation and sudden changes when performing force feedback interactive calculations.

[0123] To provide a more vivid and intuitive understanding of the construction process of the aforementioned tetrahedral mesh model, liver and gallbladder models can be selected for testing. The test results are as follows: Figure 6 As shown. For quantitative comparison, the outward offset of the liver and gallbladder was fixed at 1 and 0.2, respectively, based on their actual dimensions. Figure 6 (a), (b), and (c) in the figure represent the Delaunay tetrahedrons generated by the liver model with step sizes h=1, 0.5, and 0.1, respectively, with the number of vertices being 214, 319, and 512. Figure 6 In the diagram, (d), (e), and (f) represent the Delaunay tetrahedrons generated from the gallbladder model with step sizes h=0.1, 0.05, and 0.025, respectively, with the number of vertices as follows; the overall quality of the generated tetrahedron models is relatively good, compared to... Figure 6 As can be seen from (a), (b), (c), (d), (e), and (f), the number and quality of Delaunay tetrahedrons are controllable and change accordingly with the change of step size. The tetrahedron distribution is relatively uniform, and the overall basic unit configuration requirements of the XPBD algorithm are met.

[0124] Furthermore, in comparison Figure 6 (b), (c) or comparison Figure 6 In (d) and (e), it was found that although the step size h, as a control parameter, can directly control the number of Delaunay tetrahedra and the model quality, it lacks a quantitative evaluation of the uniformity of tetrahedral particle distribution. In some details, even the uniformity and smoothness are regressed, such as... Figure 7 As shown, observation reveals that a larger number of tetrahedra is not necessarily better. In fact, in slender model features, a step size that is too small can lead to uneven distribution along the model's extension direction. Therefore, when generating Delaunay tetrahedra, the step size h should be within an appropriate range. While achieving the required accuracy, attention should be paid to the meshing of detailed parts. In some cases, choosing a larger step size h can actually balance the uniformity of tetrahedra or tetrahedral particles, while also reducing computational costs for real-time simulation.

[0125] Step 103: Based on the vertex association relationship between each tetrahedral element and the element volume of each tetrahedral element, allocate the equivalent mass of the mesh particles of each tetrahedral element to obtain the soft mesh model.

[0126] Vertex association refers to the topological connection between a single mesh particle and all tetrahedral units containing that vertex, and is used to characterize the spatial topological properties of mesh particles in a voxel model.

[0127] The volume of a single tetrahedral unit refers to the volume occupied by a single tetrahedral unit in three-dimensional space.

[0128] The equivalent mass of a mesh particle refers to the physical mass parameter based on the topological association of mesh particles and the volume distribution of the corresponding tetrahedral elements.

[0129] A soft mesh model refers to a tetrahedral mesh model that has completed the assignment of physical mass parameters and can be directly connected to algorithms to perform real-time viscoelastic physical simulations.

[0130] In this embodiment, by traversing all tetrahedral elements in the tetrahedral mesh model, the element volume of each element is calculated based on the coordinates of its four vertices. A centroid subdivision algorithm is used to divide a single tetrahedral element into four equal subdomains, each corresponding to a vertex of the tetrahedral element. For isolated vertices associated with only a single tetrahedral element, the volume of the corresponding subdomain is used as the equivalent volume of that vertex. For higher-order vertices associated with multiple tetrahedral elements, the volume accumulation of the corresponding subdomains of all associated tetrahedral elements is added together to obtain the equivalent volume of that vertex. The equivalent mass of each mesh particle is calculated by combining the organization density parameters of the object to be interacted with. After completing the mass allocation of all vertices, the overall mass conservation can be verified, resulting in a soft mesh model suitable for simulation.

[0131] Furthermore, during the quality allocation process, tissue density field information from CT data can be combined to set differentiated density parameters for tetrahedral units of different tissue types, thereby achieving differentiated quality allocation for different tissue regions within the same organ.

[0132] In one example of the present invention, step 103 may include the following sub-steps:

[0133] Calculate the element volume of each tetrahedral element based on the coordinates of its four vertices.

[0134] If any mesh particle is associated with only a single tetrahedral element, then the equivalent mass of the mesh particle is calculated according to the first equivalent mass calculation formula and assigned to the mesh particle;

[0135] If any mesh particle is associated with at least two tetrahedral elements, the second equivalent mass calculation formula calculates the equivalent mass of the mesh particle and assigns it to the mesh particle;

[0136] Once all mesh particles are assigned equivalent mass, a soft mesh model is obtained.

[0137] After generating the tetrahedral mesh model, it is also necessary to model or define the physical constants of its tetrahedral elements. These parameters include, but are not limited to, environmental physical constants and software physical constants. Environmental physical constants include gravitational acceleration, iteration time interval / frequency, and the number of STS small steps, etc.; software physical constants include particle mass, maximum particle velocity, energy dissipation coefficient, and constants for various algorithm constraints, etc. The environmental physical constants are generally at conventional values, such as... The soft body physics constants can be determined based on the simulation object and human-computer interaction requirements, among which particle mass (i.e., vertex mass) is a fundamental quantity that needs to be determined first. Given that existing schemes also use uniform distribution or a combination of fixed distribution and particle spatial density distribution for mass allocation, problems such as mass mismatch with geometric sampling and excessive coupling between the mass field and the mesh may exist. Therefore, in this embodiment, the particle's mass is allocated based on its spatial distribution and topological relationships, according to the tetrahedral elements and the topological relationships between elements, based on the volume and density of the tetrahedral elements connected to each particle. Since the tetrahedral mesh of the soft organ has complex topological relationships, a single particle / vertex may be associated with multiple tetrahedral elements. To ensure reasonable particle mass allocation, a method of calculating and accumulating the volume based on the tetrahedral element to which the particle resides is chosen. This ensures that the mass allocation does not overflow, i.e., the total particle mass is consistent with the total mass of the tetrahedral elements.

[0138] ;

[0139] For the mass of a tetrahedral element, Let be the mass of the particle in the density field. The number of tetrahedrons. The number of particles.

[0140] Specifically, for the tetrahedral mesh model generated by tetrahedral subdivision, all tetrahedral elements within the mesh are traversed first. Based on the three-dimensional coordinates of the four vertices of each tetrahedral element, the element volume of the corresponding element is calculated. Simultaneously, the topological association between each mesh particle and all tetrahedral elements containing that vertex is established. Then, all mesh particles are traversed again, and the number of tetrahedral elements associated with each mesh particle is checked. If any mesh particle is associated with only a single tetrahedral element, the first equivalent mass calculation formula is called to calculate the equivalent mass of the mesh particle and complete the allocation. If any mesh particle is associated with at least two tetrahedral elements, the second equivalent mass calculation formula is called to calculate the equivalent mass of the mesh particle and complete the allocation. After the equivalent mass allocation of all mesh particles is completed, the overall mass conservation check is performed, and finally, a soft mesh model that can be directly used for physical simulation is obtained.

[0141] It should be noted that, as Figure 8 Point (a) in the middle The volume of each tetrahedron is... It can be determined by its vertex coordinates calculate:

[0142] ;

[0143] For isolated vertices that are associated with only a single tetrahedral element The formula for calculating the first equivalent mass is:

[0144] ;

[0145] Its equivalent volume for:

[0146] ;

[0147] More commonly, a single vertex may have topological relationships with multiple tetrahedral elements simultaneously (hereinafter referred to as "higher-order vertices"). According to the framework definition of the XPBD algorithm, multiple constraints will generate superimposed corrective gradients on particles. Therefore, the mass of particles at higher-order vertices should be reasonably allocated under multiple constraints; otherwise, it will reduce the convergence speed or even cause oscillations. Based on the algorithm's gradient calculation formula and the Lagrange operator:

[0148] ;

[0149] in, For constraints, The reciprocal of the particle mass, For Lagrange operators, These are constraint constants. It is the iteration time interval.

[0150] It can be seen that when the soft body deforms, the particle mass... reciprocal These are the weighting coefficients of the gradient magnitude. Simultaneously, the weighted total mass of the particles participating in this constraint (i.e., a pair of particles under length constraint or four particles under tetrahedral volume constraint) is an important component of the denominator of the Lagrange operator. In a rigid system, the weighted proportion of the reciprocal of the particle's mass directly affects the particle's gradient. Therefore, numerically, the larger the particle's mass, the smaller its reciprocal, the smaller the effect of the constraint on the particle's position, and the greater its effect on the positions of other particles under the same constraint (hereinafter referred to as "related particles").

[0151] Macroscopically, the larger the mass of a particle, the better it resists the effects of constraints; conversely, its associated particles are more susceptible to the influence of constraints. Furthermore, the more constraints a particle is associated with, the higher its mass should be to improve its stability under multiple constraints. Therefore, a superposition equivalent mass allocation method considering the particle space topology was designed. For example... Figure 8 As shown in (c), if point Formed into a tetrahedron , Formed into a tetrahedron And if the two tetrahedrons are adjacent, then The points then form a tetrahedron. And in the tetrahedron In this context, if the particle is a high-order vertex particle, then its equivalent volume should include... and Two parts. Using the superposition method, it can be mathematically represented as:

[0152] ;

[0153] in, This represents the number of tetrahedra associated with the particle. Let be the vector position of the particle and the vector positions of the associated particles that together form the corresponding tetrahedron. Let the volume of the tetrahedron be, for example Figure 8 As shown in (d) above. Then, combining the following second equivalent mass calculation formula, the equivalent mass of the particle can be calculated as:

[0154] ;

[0155] It should be noted that the proof that the total mass of discrete particles is equivalent to the total mass of the entire model under the same density field is presented in the following two main steps:

[0156] 1) The sum of the equivalent volumes is consistent:

[0157] First, from the perspective of a tetrahedral element, the volume of a tetrahedral element is the sum of the equivalent volumes of all four constituent particles:

[0158] ;

[0159] in This is the current tetrahedron's index. Let $\mathbf$ be the volume occupied by the corresponding particle in the current tetrahedron (hereinafter referred to as "cell allocation volume"). Then the total volume of all tetrahedrons in this tetrahedral mesh is:

[0160] ;

[0161] in This represents the total number of tetrahedrons. Secondly, from the perspective of particles, based on the method for calculating the equivalent volume of a single particle, the sum of the effective volumes of all particles in the current mesh is:

[0162] ;

[0163] in This represents the total number of all particles. Following the method of equivalent volume distribution, since the volume of each particle is the sum of the volumes of its associated tetrahedral units, the sum of all unit volumes equals both the sum of the particle volumes and the sum of the tetrahedral volumes. That is:

[0164] ;

[0165] It can be proven that the equivalent volume of a tetrahedron is relative to the equivalent volume of a particle:

[0166] ;

[0167] 2) The relationship between the tetrahedron and its associated particles in the density field:

[0168] Since both the tetrahedron and the particles are discretely distributed relative to the density field, the relationship between the tetrahedron and its associated particles needs to be studied in different cases. If the equivalent density of the tetrahedron is determined by the densities of its four associated particles, then:

[0169] ;

[0170] Substituting into the equation, we get:

[0171] ;

[0172] This ensures that the total mass of the particles matches the total mass of the tetrahedral unit cells. For soft organs, the density variation of an individual organ or tissue is within 0.1 g / cm³. 3 Within this structure, uneven distribution is unlikely to occur. Therefore, when constructing the basic tetrahedral mesh, it is necessary to configure it independently according to the type of organ and tissue, and the algorithm must ensure independent iteration for different organs. Under this premise, it can be ensured that the particle equivalent mass of this allocation method is highly close to the mass of the entire soft organ.

[0173] In another example of the present invention, in order to evaluate the quality of the generated soft mesh model, the method may further include the following steps:

[0174] After obtaining the soft mesh model, traverse all tetrahedral elements within the soft mesh model and calculate the nearest neighbor distance probability density distribution index and dihedral angle quality index one by one.

[0175] If there are unqualified tetrahedral elements whose nearest neighbor distance probability density distribution index and / or dihedral angle quality index do not meet the preset standards, and the number of unqualified tetrahedral elements is less than the preset element threshold, then the region to be adjusted where the unqualified tetrahedral elements are located is located, and the polygon bounding volume corresponding to the region to be adjusted is updated according to the new voxel mesh step size.

[0176] Jump to the step of performing tetrahedral subdivision of the polygon bounding volume to generate a tetrahedral mesh model;

[0177] If there are unqualified tetrahedral elements whose nearest neighbor distance probability density distribution index and / or dihedral angle quality index do not meet the preset standards, and the number of unqualified tetrahedral elements is not less than the preset element threshold, then obtain a new voxel mesh step size and a new outward offset, and jump to the step of constructing a polygon bounding volume using CT data.

[0178] If the nearest neighbor distance probability density distribution index and dihedral angle quality index of all tetrahedral elements meet the preset standards, then the response control command is executed to apply external force to the soft mesh model.

[0179] In this embodiment, after obtaining the soft mesh model by completing the vertex equivalent quality allocation, all tetrahedral elements within the soft mesh model are traversed. For each tetrahedral element, the corresponding nearest neighbor distance probability density distribution index and dihedral angle quality index are calculated one by one. The two indices are compared with preset compliance conditions to complete the qualification verification of all tetrahedral elements. If the verification finds that there are unqualified tetrahedral elements whose nearest neighbor distance probability density distribution index and / or dihedral angle quality index do not meet the preset compliance conditions, and the total number of unqualified tetrahedral elements is less than the preset element threshold, then the adjustment region where all unqualified tetrahedral elements are located is located by spatial clustering algorithm. A new voxel mesh step size corresponding to the unqualified index is matched, and the polygon bounding volume fragment corresponding to the adjustment region is locally retopologically re-updated and updated. Then, the step of tetrahedral subdivision of the updated polygon bounding volume and generation of tetrahedral mesh model is executed to complete the regeneration and verification of the locally optimized mesh. If the verification finds that the nearest neighbor distance probability density distribution index and / or dihedral angle quality index do not meet the preset compliance conditions, the adjustment region is further verified. If any tetrahedral elements fail to meet the preset standards for dihedral angle quality, and the total number of such elements is not less than the preset element threshold, then a new voxel mesh step size and a new outward offset are obtained based on the deviation of the non-compliance index. The process then proceeds to the step of constructing a polygon bounding volume using CT data to complete the global reconstruction and re-verification of the mesh. If the verification confirms that the nearest neighbor distance probability density distribution index and dihedral angle quality index of all tetrahedral elements meet the preset standards, then the soft mesh model is deemed qualified, and the subsequent response control command is executed to apply external force to the soft mesh model.

[0180] In practical implementation, taking a virtual endoscopic surgery scenario as an example, due to the complex shapes and irregular geometry of soft organs, such as the spinal nucleus pulposus and gallbladder, it is impossible to model them using uniform parameters. Currently, discretization and particleization remain the main methods. In this embodiment, by introducing the Nearest-Neighbor Distance Probability Density Function (ND-PDF) from finite element analysis, and combining it with the purpose of Delaunay tetrahedral partitioning, the basic voxel tetrahedron is quantitatively and comprehensively evaluated using position density and angular quality parameters as core indicators.

[0181] Based on the characteristics of the tetrahedral model of XPBD soft organs, the nearest neighbor distance probability density function is suitable for evaluating the uniformity of its mesh model. Its specific definition is as follows. For the tetrahedral mesh... The coordinates of a vertex can be represented as:

[0182] ;

[0183] in, This refers to a three-dimensional real-valued position vector. The process involves calculating the distance between each vertex and selecting the nearest neighbor distances to obtain a set. :

[0184] ;

[0185] in The distance between a vertex and its nearest neighbor:

[0186] ;

[0187] Subsequently obtained When performing distribution statistics on the set, for the tetrahedral distribution used in the XPBD algorithm, it is necessary to focus on the overall situation and uniformity of the tetrahedrons, while also paying attention to the extreme cases of the densest / sparsest parts. Therefore, the nearest neighbor distance probability density distribution index considered is as follows:

[0188] (1) Probability of maximum / minimum, mean and median and Because the algorithm is relatively stable, its requirements for tetrahedral geometry are not overly stringent. At most, oscillations can cause convergence instability. Therefore, the primary statistical methods for the distribution of nearest neighbor distance sets should be the mean and median. Generally, the mean better reflects the main characteristics of the set, while the median better reflects some extreme cases; both are important indicators.

[0189] (2) Extreme value probability , , and For the metric used to evaluate the uniformity of spatial distribution, the extreme value probability includes points in the densest and sparsest regions. This embodiment uses the distribution probabilities of the smallest 1% and largest 5% nearest neighbor distances, commonly used in engineering, for statistical analysis. On the one hand, excessively dense vertices are a significant cause of local oscillations in the algorithm, especially in a few extreme cases of high density; therefore, the smallest 1% metric is used. On the other hand, the algorithm allows for large deformations, thus allowing for a small number of extreme cases with relatively large distances. The focus is on statistically analyzing the distribution of the 95% nearest neighbor distance; therefore, the 95% distribution metric is selected.

[0190] (3) Standard Deviation of Distribution Width. As a typical indicator for evaluating the uniformity of the distribution of tetrahedron vertices, the standard deviation of distribution width is chosen as the evaluation indicator. The smaller the standard deviation, the more uniform the tetrahedron model is. It is mainly used for the horizontal comparison of different models.

[0191] (4) Probability density distribution function curve. To gain a more intuitive understanding of the distribution of nearest neighbor distances, direct analysis using the probability density function distribution curve is also an important method. By observing the curve, one can determine whether it is unimodal, long-tailed, or multimodal. In this embodiment, histogram-based estimation (HBE) is used to construct the probability density function. Let the set... The interval is Then the interval is divided into equal-width intervals Each equal-width interval Width is:

[0192] ;

[0193] Secondly, statistics for each Number of distances within Obtain the distribution set Finally, the probability density function is evaluated as follows:

[0194] ;

[0195] in, For this interval The midpoint of the probability density. Using this revised definition, we can essentially guarantee that the integral of the probability density is close to 1, that is:

[0196] ;

[0197] Because the aforementioned method for generating tetrahedral models has parameterized characteristics and ensures that the probability density integral is 1, it is beneficial for lateral comparison of different nearest neighbor distance sets, and is particularly important for comparing the quality of meshes generated by different parameters for the same shape.

[0198] At the same time, it is also necessary to calculate the quality index of the dihedral angle:

[0199] In this embodiment, a dihedral angle specifically refers to the angle formed by two adjacent triangular faces of any side of a tetrahedral element. Let the vector coordinates of the four vertices of the tetrahedron be... Then one of the edges Two adjacent faces are The normal directions of both are:

[0200] ;

[0201] Secondly, the unit vector of this edge is:

[0202] ;

[0203] The perpendicular vector components of the normals of two adjacent triangular faces on that side are:

[0204] ;

[0205] Finally, the tetrahedral edge The dihedral angle between adjacent faces is:

[0206] ;

[0207] Using the above definition and calculation method, we can traverse each tetrahedral element to obtain the size of its six dihedral angles, thus obtaining the set of dihedral angles of the tetrahedral mesh to be evaluated. Based on the characteristics of the algorithm, it is necessary to set statistical indicators for the obtained dihedral angle set. The specific indicators and reasons are as follows.

[0208] (1) Global distribution statistics of dihedral angles. Global distribution statistics of dihedral angles are helpful in evaluating the overall situation of tetrahedral meshes. The main indicators of global distribution statistics include the maximum / minimum value of angles, the average value, and the angle distributions of 1%, 5%, 95%, and 99%. These global statistical indicators are helpful in reflecting the overall angle situation of tetrahedral meshes;

[0209] (2) Unit extreme value index. In addition, as a unit of the network, the maximum / minimum value of the dihedral angle of each tetrahedron is also important. Except for some tetrahedrons located near the surface of the shape, the dihedral angle of the main body should not be too large or too small, otherwise tetrahedron degeneration or even negative volume will occur.

[0210] (3) Pass rate index. For a tetrahedral network, since the number of its dihedral angles is 6 times that of the tetrahedral elements, the large number makes it difficult to analyze them one by one. Therefore, the pass range is selected. This acceptable range is defined based on experience from typical PBD physics simulations [3,11]. Based on this acceptable range, the acceptable rate at the angular level is defined. (Angle Pass Rate) and the pass rate at the tetrahedral unit level (Tet Pass Rate), calculated as follows:

[0211] ;

[0212] in, This represents the number of qualified dihedral angles. Represents the total number of dihedral angles; This represents the number of tetrahedral elements that meet all dihedral requirements. This represents the total number of tetrahedral elements. By defining the pass rate, the dihedral pass rate of tetrahedral meshes with different element counts can be effectively compared, making it an important dihedral pass rate indicator.

[0213] Based on the above indicators, the proposed evaluation method for tetrahedral meshes used for virtual surgical soft organ simulation mainly includes two parts: nearest neighbor distance probability density distribution and dihedral quality. It can evaluate the overall, local, and extreme regions of the tetrahedral mesh from the perspectives of global indicators, extreme values, and HBE distribution, as detailed in Table 1.

[0214] Table 1

[0215]

[0216] Furthermore, based on the geometric and numerical definitions and requirements of the aforementioned evaluation methods, a secondary development plugin for Blender software can be written using Python. This process mainly involves five steps. The main implementation logic of the plugin is as follows: Figure 9 As shown.

[0217] Step 1: Tetrahedral Identification of Voxel Mesh. Before calculating the metrics, it is necessary to identify and extract point groups and their topological relationships using tetrahedrals as the unit. First, the mesh particles are uniformly converted to world coordinates; second, all vertex groups are traversed and tetrahedral elements with 4 points are selected. If there are duplicate points or non-tetrahedral elements in a 4-vertex group, they are skipped; finally, the number of identified elements is counted, including non-tetrahedrals (non_quad) or degenerate tetrahedrals (bad_quad), to indicate possible abnormal elements in the mesh.

[0218] Step 2: Nearest Neighbor Distance Calculation. Utilizing the built-in KDTree acceleration structure, quickly search for the nearest neighbor vertex of each vertex. First, insert all vertices into the KDTree; second, traverse each vertex, obtaining its two nearest points (the vertex itself at a distance of 0 and its nearest neighbor vertex), with the second nearest vertex being the nearest neighbor vertex; third, calculate and record the nearest neighbor distance for each vertex and obtain the specific value of the set; finally, transform the set into a probability density function based on the sampling results, divide it into equal-width intervals, normalize the interval frequencies, calculate and record the probability value sets.

[0219] Step 3: Statistics and calculation of key indicators for nearest neighbor distance. Based on the identified probability value groups, calculate the corresponding indicators, including the probability of the mean, the probability of the median, the probability of the extreme values, and the standard deviation (Std), and form the node values ​​of the NND-PDF curve.

[0220] Step 4: Dihedral Angle Distribution Statistics. First, based on the data identified in Step 1, traverse the edges of the tetrahedral element and calculate the unit normal of a pair of adjacent faces; second, perform inverse cosine calculation on this normal pair to obtain the dihedral angle value and perform interval verification (to avoid floating-point calculation errors that may cause calculation results that violate the definition); finally, refer to the data structure of Wavefront OBJf to establish a tetrahedral element-vertex index-angle value dihedral angle data group.

[0221] Step 5: Statistics and Calculation of Key Dihedral Angle Indicators. Based on the statistically collected dihedral angle data, calculate the corresponding indicators at two levels (angle level and tetrahedral unit level), including maximum / minimum values, average values, 1%, 5%, 95%, and 99% angle distributions, as well as the maximum / minimum value of the dihedral angle for each tetrahedral unit. Finally, statistically calculate the pass rate at the angle level and the pass rate at the tetrahedral unit level.

[0222] Step 6: Generate a text block using the Scripting workspace in the form of a Text Editor, and output the key statistical indicators mentioned above.

[0223] Furthermore, to verify the effectiveness of the evaluation method, this embodiment uses a tetrahedral mesh geometric distortion control experiment for testing. A typical soft organ generated tetrahedral model from virtual surgery is selected as the control group; and the local mesh of this model is geometrically distorted to form the experimental group. Simultaneously, the evaluation method of this embodiment is used to calculate the indicators for both groups of models, and the significance of key indicators is compared. On the one hand, the tetrahedral quality of the generated model can be verified by using existing methods to evaluate and analyze it; on the other hand, the effectiveness of the tetrahedral mesh evaluation method proposed in this embodiment is verified by comparing the undistorted and distorted experimental groups.

[0224] The experiment selected a typical liver model tetrahedron as the test object, such as... Figure 10 As shown in (a) of the table, a moderate number of vertices was chosen for ease of distortion processing and comparison; the specific parameters are shown in Table 2. The experimental group performed local distortion processing on the test model, including local density adjustment (Experiment I) and local angle adjustment (Experiment II), as shown in Table 2. Figure 10 As shown in (b) of 10 and (c) of 10. The control test experiment was conducted using the rapidly developed plugin described above. The test results and nearest neighbor probability distance density distribution index are shown in Table 2. Figure 11 The dihedral angle distribution index is shown in Table 3.

[0225] Table 2

[0226]

[0227] Table 2 shows the comparative experimental statistics obtained using the nearest neighbor distance probability density evaluation method. All three models identified 317 vertices and 1177 tetrahedral elements, with no abnormal tetrahedral elements detected. This indicates that the tetrahedral meshes generated using the Delaunay tetrahedral subdivision method do not contain degenerate or substandard tetrahedral elements, proving the effectiveness of the method. Since only a small number of points in the distorted mesh group were modified, the overall distribution range of the three models is similar. The standard deviation is around 0.3. Divided into 60 segments, the widths of the equal-width intervals are also quite similar, at 0.022, 0.028, and 0.026 respectively. The average values ​​of the standard mesh in the control group and the mesh in Experiment II (which only underwent angular distortion) are... With median All values ​​were above 90%, indicating that the nearest neighbor distances of the mesh are relatively concentrated, with most vertices being very close together. However, in Experiment I, where positional distortion was performed, the width of the equal-width intervals increased to some extent, resulting in a minimum distance of 0.07, with the average value... With median All showed a significant decrease, reaching only 78.71%. From Figure 11 The nearest neighbor distance probability density distribution also shows that Figure 11 The control group in (a) and Figure 11 In group II of angular distortion (c), the nearest neighbor distances are relatively concentrated in the ranges of 0.6–0.8 and 1.6–1.8, while... Figure 11 In group I of positional distortion (b), the distribution is more uniform and lacks concentration; secondly, from Figure 11As can be seen from the comparison of the three in (d), both experimental groups I and II exhibit varying degrees of distribution within the small distance range of 0-0.4. Particularly in group II, the extremely small distribution of positional distortion can lead to uneven length constraints and even numerical overflow, indicating poor model quality. In summary, this paper proposes using nearest-neighbor distance probability density to identify local positional distortion in tetrahedral meshes, enabling objective quantitative evaluation of the positional and distance quality of the tetrahedral mesh. This method is highly sensitive to local anomalies and can effectively locate a small number of positionally distorted regions within the tetrahedral mesh. Through probability density integral normalization, the obtained index exhibits scale consistency and can be used for cross-sectional comparative analysis between different mesh models.

[0228] Table 3

[0229]

[0230] Table 4

[0231]

[0232] Table 5

[0233]

[0234] Using the aforementioned method for dihedral angle distribution statistics, the full-mesh dihedral angle index results for the experimental mesh object are shown in Table 3, with a comparison of their specific values. Figure 12 The bar chart shows the dihedral angle pass rate and the maximum / minimum dihedral angle statistics, respectively. The algorithm plugin identified 7062 dihedral angles in all three models, corresponding to the number of tetrahedral elements in the aforementioned nearest neighbor distance probability density distribution statistics, indicating accuracy in identifying tetrahedral elements and dihedral angles. Comparing the dihedral angle statistics of the three models, since distortion processing only addresses a few local points, the overall dihedral angle probability density distribution is consistent, with the maximum angle approaching 180 degrees and the average around 110 degrees. The pass rate statistics in Table 4 show that, whether considering dihedral angles across the entire mesh or dihedral angles at the tetrahedral level, most dihedral angles fall within the pass range, with a pass rate exceeding 90%, even reaching 95%. This indicates that the tetrahedral generation method presented in this paper is of high quality and the mesh dihedral angle distribution is reasonable. In Experiment I, due to position distortion, a small number of dihedral angle anomalies were also generated, resulting in a lower dihedral angle pass rate. Meanwhile, since distortion processing was performed in Experiment I and Experiment II, such as Figure 12As shown, with most indicators remaining consistent, the minimum dihedral angle, a crucial indicator, exhibited significant anomalies. The minimum angle in the control group was 0.8 degrees, but in the position distortion experiment group I, it was 0.5 degrees, and in the angle distortion experiment group II, it was even more severe at 0.2 degrees. This indicates that the corresponding tetrahedral elements may have reached the edge of degenerating into planes, which is clearly detrimental to the stable simulation of the subsequent XPBD algorithm. In summary, the distortion test results demonstrate that the proposed method of using dihedral angle statistics to evaluate tetrahedral mesh quality can effectively locate local angle distortions in tetrahedral angles. The intuitive comparison of key indicators via bar charts allows for quick identification of extreme dihedral angle cases and facilitates the comparison of tetrahedral quality across different meshes.

[0235] Step 104: Apply external forces to the soft mesh model in response to control commands;

[0236] Control commands refer to digital instructions sent by users through force feedback operating devices, interactive handles, and other input terminals to control virtual instruments to perform corresponding operations, including but not limited to instrument position adjustment, clamp opening and closing, traction, puncture, and other operation commands.

[0237] The application of external force refers to the process of mapping and loading the external forces such as the mechanical interaction force and environmental force corresponding to the control command onto the corresponding mesh particles of the soft mesh model. It is the physical input link that drives the soft model to produce deformation.

[0238] In this embodiment, control commands sent by the user through the interactive terminal are received in real time. The virtual instrument type, pose change, operation type, and operation range corresponding to the control commands are parsed. The contact area and contact vertices between the virtual surgical instrument and the soft mesh model are determined by a spatial collision detection algorithm. Based on the operation type and instrument pose change, the corresponding contact force, clamping force, tensile force, and compressive force are calculated. At the same time, external environmental forces such as gravity and fixed boundary constraints are applied to the soft mesh model. All external forces are mapped to the corresponding mesh particles of the soft mesh model according to their interaction relationship to complete the application of external forces.

[0239] Furthermore, the control commands can be pre-identified to determine the user's operation type and range. Collision detection pre-calculation and external force pre-allocation can then be performed on the predicted contact area, reducing computational latency during the interaction process. Additionally, based on the force feedback device's sensing threshold, minute forces below the threshold can be filtered to avoid interference from ineffective external forces in the simulation process.

[0240] Step 105: Based on the equivalent mass of the mesh particles and the force corresponding to the action of the applied external force, the particle position of each tetrahedral element is iteratively updated to obtain the updated particle position.

[0241] Particle update position refers to the final spatial position of the mesh particles locked after the full-process iterative calculation of a single simulation time step is completed. It is the physical data source for software deformation visualization rendering.

[0242] In this embodiment, the equivalent mass of the mesh particle and the applied force are used as inputs. The improved extended position-based dynamics (I-XPBD) algorithm framework is used to perform iterative calculations. First, the motion velocity of the mesh particle is updated by the external force and the equivalent mass. The dynamic velocity dissipation coefficient is calculated in combination with the soft deformation duration to perform attenuation correction on the velocity. The initial position of the mesh particle is predicted based on the corrected velocity. A single simulation time step is divided into multiple small time steps (STS). In each small step, Gauss-Seidel iterative solutions for length constraints, volume constraints, and collision penetration resistance constraints are performed in sequence according to a preset priority. During the solution process, the dynamic compliance constant is calculated in combination with the soft deformation duration to simulate viscoelastic properties. After completing the iterative solution of all STS small steps, the final spatial position of the mesh particle is locked as the particle update position.

[0243] In one example of the present invention, step 105 may include the following sub-steps:

[0244] The current velocity of each tetrahedral element is updated based on the force corresponding to the action of the external force application, and the dynamic velocity dissipation coefficient is calculated in combination with the duration of soft deformation.

[0245] The initial velocity of each vertex is attenuated and corrected using the dynamic velocity dissipation coefficients to obtain the particle corrected velocity;

[0246] The initial predicted position of each particle in the current frame is predicted based on the corrected velocity of each particle.

[0247] The dynamic compliance constant is calculated based on the duration of soft deformation at each simulation time step, and the Lagrange operator of the length constraint is adjusted based on the dynamic compliance constant to obtain the updated length constraint;

[0248] The initial predicted position is iteratively corrected under the constraints of length update, mention constraint, and collision penetration prevention constraint to obtain the updated position of the particle and update the current velocity of the particle.

[0249] Optionally, the formula for calculating the dynamic velocity dissipation coefficient is:

[0250] ;

[0251] in, The dynamic velocity dissipation coefficient, The velocity dissipation time-varying coefficient, t is the preset base dissipation coefficient, and t is the duration of the soft deformation.

[0252] Optionally, the formula for calculating the dynamic compliance constant is:

[0253] ;

[0254] Where k, b, and a are preset KV model constants. Corrects the particle velocity.

[0255] Force refers to all external forces applied to the vertices of a soft mesh model, including environmental gravity, contact forces, clamping forces, tensile and compressive forces of virtual instruments, and reaction forces of fixed boundary constraints.

[0256] The current particle velocity refers to the velocity of the mesh particles locked after the previous frame of physics simulation iteration. If it is the first simulation, it is zero.

[0257] The duration of soft deformation refers to the cumulative time from the start of significant deformation in the soft mesh model to the current simulation frame.

[0258] The dynamic velocity dissipation coefficient refers to the velocity decay coefficient that changes dynamically with the duration of soft deformation. It is used to simulate the energy dissipation characteristics during viscoelastic damped oscillation and to control the decay rate of the oscillation amplitude.

[0259] The initial vertex velocity refers to the velocity updated solely based on external forces and the current particle velocity, without any dissipation correction. The particle-corrected velocity refers to the velocity after being decayed by a dynamic velocity dissipation coefficient, used to predict the initial vertex position in the current frame. The initial predicted position refers to the initial spatial position of the vertex in the current frame predicted based on the particle-corrected velocity and the particle position of the previous frame.

[0260] The simulation time step refers to the time interval of a single complete iteration of a physical simulation.

[0261] The dynamic compliance constant refers to the constraint compliance parameter that changes dynamically with the duration of soft deformation. It is used to adjust the rigidity of length constraints and simulate the stress relaxation and creep characteristics of viscoelasticity.

[0262] Length constraints refer to the core geometric constraints that maintain the side length of tetrahedral elements and limit excessive stretching or compression deformation of the soft body. Volume constraints refer to the geometric constraints that maintain the initial volume of tetrahedral elements and ensure the incompressibility of the soft organ. Collision penetration prevention constraints refer to the interaction constraints that prevent virtual instrument colliders from penetrating the soft body's surface particles, used to ensure the stability of interactive simulation.

[0263] The particle update position refers to the final spatial position of the current frame vertex after all constraint iterations have been completed.

[0264] In this embodiment, based on the force corresponding to the action of the applied external force, the current velocity of each tetrahedral unit particle is first updated to obtain the initial velocity of the vertex. At the same time, the corresponding dynamic velocity dissipation coefficient is calculated in combination with the soft deformation duration. The initial velocity of the vertex is attenuated and corrected using the dynamic velocity dissipation coefficient to obtain the particle corrected velocity. Then, based on the particle corrected velocity and the particle position of the previous frame, the initial predicted position of each vertex in the current frame is predicted. Subsequently, a single simulation time step is divided into multiple STS sub-steps. In each STS sub-step, the dynamic compliance constant is first calculated in combination with the soft deformation duration. Based on the dynamic compliance constant, the Lagrange operator of the length constraint is adjusted to obtain the updated length constraint with viscoelastic properties. Then, according to the preset priority, the initial predicted position is corrected by Gauss-Seidel iteration under the updated length constraint, volume constraint, and collision penetration prevention constraint to obtain the middle position after each STS sub-step. After completing all STS sub-step iterations, the particle update position is locked. Finally, the particle current velocity is updated based on the particle update position and the particle position of the previous frame to prepare for the next frame iteration.

[0265] like Figure 13 As shown, this diagram illustrates the damped oscillation phenomenon of a soft mesh model after a force is applied in an embodiment of the present invention. This phenomenon is the behavior that occurs after a soft organ or tissue is subjected to traction, and it mainly consists of three stages. Stage 1 is as follows: Figure 13 As shown in (a), in the initial stage of release, a large amount of the soft potential energy is converted into kinetic energy, thus the soft body produces a relatively large amplitude oscillation; stage 2 is as follows Figure 13 As shown in (b), the oscillating behavior consumes some kinetic energy over time, thus damping becomes dominant, leading to a decrease in oscillation amplitude. In the final stage (3), the soft body tends to stabilize over time. The dominant internal cause of this behavior is the change in stress in the biomaterial over time, manifested as the process of soft body oscillation over time, referred to as damped oscillation. It is worth noting that this behavior differs from the common damped behavior of elastic springs; it is a nonlinear deformation behavior, unique to biomaterials, and highly correlated with the time dimension. Damped oscillation phenomena have quantifiable indicators such as the maximum amplitude ratio, observable amplitude, number of cycles, and settling time.

[0266] like Figure 14 As shown in (a), after the soft organ is released by a traction operation, the clamping point first undergoes a maximum-amplitude asymmetric reciprocating motion. The maximum positive distance of this reciprocating motion is denoted by [missing information - likely a typo]. This indicates that the return trip has one maximum reverse travel amplitude, expressed as... The initial oscillation of the damped oscillation behavior of soft organs is an important basis for doctors to observe and judge the degree of softness and hardness. Therefore, this article uses the maximum amplitude ratio of this stage (abbreviated as...). (Amplitude Ratio, specifically referring to the amplitude ratio of the elastic-viscosity behavior of soft bodies in this paper) is an important parameter for comparing soft body simulation algorithms with real organ tissues. Using a ratio approach can reduce the influence of different scales on the parameter; that is, the maximum amplitude is:

[0267] ;

[0268] Observability of the simulation process is an important quantitative indicator, and the observable reciprocating amplitude of damped oscillation behavior is one of them. On the one hand, using observability as a metric is an important quantitative method for most real-time simulation algorithms for human-computer interaction. Furthermore, defining the resolution of observability can further quantify the definition of "stability" of the software simulation algorithm. On the other hand, based on this definition, quantitative indicators can be provided for the algorithm's accuracy requirements, which can further provide effective judgment support for algorithm acceleration and optimization, and also provide a basis for subsequent comparative experiments of the algorithm. For example... Figure 14 As shown in (b) of the figure, the minimum observable oscillation distance and the maximum positive distance The ratio is the lowest observable amplitude ratio:

[0269] ;

[0270] When the damped oscillation amplitude ratio of the soft body is less than this value, the soft body enters a steady state. Furthermore, the time required for the soft organ to go from deformation release to a steady state is called the settling time (ST), and the number of reciprocations during this period is called the observable reciprocation count.

[0271] In this embodiment, a combination of volume constraints and length constraints is used to simulate the damped oscillation phenomenon of a soft mesh model under applied forces. The driving force of elastoviscous dynamics is the internal stress of the soft body. Under continuous internal stress, the deformation and other parameters of the soft body undergo dynamic changes. The stiffness of the soft body gradually changes with the deformation time, the magnitude of the applied force, the degree of deformation, and the direction. According to the form of Newton's laws of motion, the effects of force and "constraints" can be considered equivalent, and the equivalent internal force is:

[0272] ;

[0273] in, For the quality matrix, The position of the particle. To constrain the equivalent internal forces of the particles, Let be the potential energy function.

[0274] Without considering external forces, The terms are formally identical to those used in force calculations; therefore, based on the constraint equations of XPBD, terms equivalent to forces can be derived, thus leading to the expression for the equivalent internal force. Length constraints are the primary constraint; therefore, in this embodiment, internal force calculations mainly focus on length constraints. If the position vectors of the two particles in the length constraint are respectively... and The original length between them is Then, for this pair of particles, the length constraint can be expressed as:

[0275] ;

[0276] When the distance between the two particles changes, according to the XPBD algorithm, the position constraint gradients that the two particles need to generate are:

[0277] ;

[0278] in, The reciprocal of the particle's mass. For the Lagrange operator in the XPBD algorithm, The absolute value of the distance between the two particles. As a compliance component, particles are allowed to deviate from the constraints to a certain extent. It is the compliance constraint constant, which can be regarded as the reciprocal of the softness or hardness of the software. When the value is 0, the compliance component is 0, representing that the object is a pure rigid body, and the particles should follow the... The direction is restored as quickly as possible.

[0279] Taking the second derivative of the above position constraint gradient formula with respect to time, we can use the difference approximation derivative to obtain the internal force expression of the XPBD length constraint:

[0280] ;

[0281] It can be seen that the internal forces generated by the length constraint are related to the Lagrange operator. This is related to the position vector difference between the two particles. According to elasticity, if an object is infinitely close to a rigid body, when it is stretched and deformed, it will immediately return to its original shape, and the internal force of the restored object will tend to infinity, that is, there is... Combining the formulas, we can obtain:

[0282] ;

[0283] in The sum of the masses of the particles. The factor set to simplify the expression is defined as follows:

[0284] ;

[0285] The above simultaneous formulas are divided into two parts: rigid body terms. and software items In the form of. When When it approaches 0, The soft term also tends towards 0, at which point the algorithm fully exhibits rigid body characteristics, i.e., determined by the rigid body term. The internal forces are determined. Conversely, when length constraints simulate the hardness, softness, and viscoelasticity of non-rigid objects, the rigid body terms remain constant, while the soft body properties mainly depend on... The software component is effective, and and They have a linear relationship. The larger the size, the softer the constraint, and the greater the internal force. The smaller, the opposite. When = 0, internal force The maximum internal force, equal to the rigid body limit, is:

[0286] ;

[0287] Therefore, for software with length constraints as the primary constraint, in order to achieve viscoelastic software simulation without sacrificing excessive computational cost, the compliance constant of the constraints can be dynamically modified. To achieve dynamic control Instead of directly calculating each constraint in real time. This allows for the efficient simulation of elastic-viscosity soft body deformation.

[0288] In position-constrained gradients, the rate of recovery of length-constrained deformation is the reciprocal of the particle's mass. Lagrange operators and gradient The decision is made. The time dimension is the primary dimension in software simulation. The viscoelasticity of software is reflected in the changes in stress and strain over time. According to classical mechanics, time is not the driving force behind soft deformation. Therefore, in real-time simulation, a time-stiffness relationship constraint is not directly established; instead, a stress-time relationship is established to simulate the elastic-viscosity of the software. Secondly, in the stress and strain dimensions, elastic-viscosity software exhibits a state where stress decreases and strain increases over time. These two aspects are highly coupled in the algorithm and are highly dependent on time. To reflect elastic-viscosity in XPBD simulation software, it is possible to establish... Dynamically changing in relation to time This is essentially achieved through a Kelvin-Voigt model (KV) or a Generalized Maxwell model. Since software algorithms are designed for applications with complex shapes and operations, aiming for the most concise and frequent computations possible, the mathematical form of the simplified and efficient KV model is chosen as a crucial reference, combined with the XPBD algorithm. From the relationship, we can obtain:

[0289] ;

[0290] in , , These are KV model constants. For time, For velocity. Considering the small deformation scale and low operation speed in endoscopic surgery, the soft body deformation process is a quasi-static process. Therefore, the velocity term is ignored when calculating the internal force. Differentiating the above equation, we get:

[0291] ;

[0292] In the viscoelastic dynamic simulation of soft materials, the compliance constraint constant It exhibits a hyperbolic relationship with time. It's worth noting that the elastic-viscosity simulation here does not consider the plastic deformation of the soft body. That is, regardless of the amount of deformation, when the external force is removed, the length constraint will always cause the particle to eventually return to its original static length, and the plastic deformation effect can be directly and linearly adjusted according to the amount of deformation. That's all.

[0293] Based on the aforementioned analysis of damped oscillation behavior, the oscillation amplitude is large in the initial stage, with a relatively small velocity dissipation ratio during the initial deformation. However, as time changes, the velocity dissipation ratio increases, the oscillation phenomenon gradually weakens, and the deformation tends to stabilize. Based on this, and according to the mathematical form of the KV model and the formula for calculating equivalent internal forces, under the condition that the energy dissipation process is stable, the velocity dissipation coefficient... It should have a linear positive correlation with time, that is:

[0294] ;

[0295] in, The velocity dissipation time-varying coefficient, The basic dissipation coefficient can be obtained from actual physical experiments. Therefore, the optimized pre-iteration velocity (i.e., the particle correction velocity) should be:

[0296] ;

[0297] Step 106: Render the scene based on the updated position of each particle, and output the deformation scene corresponding to the object to be interacted with.

[0298] Deformation visuals refer to visual images that recreate the deformation state of an interactive object under external forces, matching a specific viewpoint and visual characteristics, and are used to provide users with real-time feedback on the effects of virtual operations. This specific viewpoint can be the viewpoint used in clinical endoscopic surgery.

[0299] In this embodiment, the current spatial position of all vertices of the surface rendering mesh is calculated by linear interpolation of particle update position. Simultaneously, the normal vector, tangent space and other rendering auxiliary data of the surface rendering mesh are updated. Based on the optical parameters and pose data of the clinical endoscope lens, the model-view-projection transformation is performed to convert the three-dimensional vertex coordinates into two-dimensional coordinates on the screen. After rasterization and soft tissue-specific pixel shading processing, visual enhancement effects such as highlighting of the stress area and visualization of stress distribution are superimposed. Finally, the deformation image of the interactive object that matches the visual experience of clinical endoscopic surgery is output.

[0300] In one example of the present invention, step 106 may include the following sub-steps:

[0301] The position of each particle is updated and transformed to screen pixel coordinates;

[0302] Rasterize the triangular facet formed by the vertices of the particle update position to obtain multiple pixels and associate them with the screen pixel coordinates;

[0303] Each pixel is assigned brightness and / or color to generate and output the deformed image corresponding to the object to be interacted with.

[0304] Screen pixel coordinates refer to the pixel position coordinates of a vertex in three-dimensional space after model-view-projection transformation and mapped onto a two-dimensional display screen.

[0305] A triangular facet is the smallest two-dimensional geometric unit that constitutes the watertight surface rendering mesh. It is enclosed by the vertices to which the particle update position belongs and is the basic unit of rasterization processing.

[0306] A pixel is the smallest image display unit on a two-dimensional display screen. After rasterization, it corresponds one-to-one with the screen pixel coordinates and is the object of processing for assigning brightness and color.

[0307] In this embodiment, after obtaining the particle update positions of all surface-rendered mesh particles calculated by the I-XPBD algorithm in the current frame, model transformation, view transformation, and projection transformation are performed sequentially based on the lens optical parameters and pose data of the current screen display perspective. This converts the three-dimensional particle update position of each vertex into the corresponding two-dimensional screen pixel coordinates. Then, all triangular faces of the surface-rendered mesh are traversed, and rasterization is performed on each triangular face composed of vertices to which the particle update position belongs. The screen pixel range covered by the triangular face is determined, and each pixel within the range is associated with the corresponding screen pixel coordinates. Attributes such as vertex normals and texture coordinates are assigned to each pixel through linear interpolation. Subsequently, the brightness of each pixel is assigned by combining the position and intensity of the endoscope cold light source with the vertex normal associated with the pixel. At the same time, the color of each pixel is assigned by combining the organization type, material characteristics, and stress state of the object to be interacted with. The brightness and color values ​​are then fused to generate the deformation image corresponding to the object to be interacted with, and finally output to the display terminal.

[0308] In practical implementation, to verify the damped oscillation behavior of soft organs after being subjected to stretching and other operations, and the characteristics of this behavior during actual deformation, the quantitative parameters of elasticity and viscosity in the damped oscillation behavior are measured. This allows for the quantitative comparison of optimized soft algorithms. A corresponding soft deformation experimental platform needs to be built, and corresponding soft experimental data needs to be designed and collected. For example... Figure 15 As shown, the experimental platform uses an industrial-grade robot as the operating actuator to simulate the traction operations performed by doctors on soft organs such as the liver and gallbladder, while ensuring sufficient positional accuracy and rigidity to avoid additional errors or vibrations. Secondly, a high-precision six-dimensional force sensor (ME) is used to measure the robot's operating force on the soft material, further controlling the robot's tension and compression under specific experimental settings. Thirdly, an RGBD binocular camera capable of detecting depth and point cloud information is used to effectively collect and locate the dynamic changes of key points / local areas of the soft material during the experiment. Finally, a robot teaching controller is used for position control, and a computer records various measurement information in real time. The key indicators of the relevant experimental instruments and equipment are shown in Table 6.

[0309] Table 6

[0310]

[0311] During the verification process, fresh pig liver was selected as the experimental subject and fixed to the experimental table using a custom-designed clamp. To ensure that the robot closely approximates actual traction operations, the grippers used in actual laparoscopic surgery were modified by adding a pair of symmetrical springs to adjust the clamping force. A pin device was designed to control the opening and closing of the grippers and was assembled onto the end effector of a robotic arm equipped with a high-precision force sensor. This enabled the robotic arm to grip, pull, and release the experimental pig liver. The experimental process is as follows: Figure 16As shown, to accurately capture the location and trajectory of key liver positions, white markers were added to the clamping area. A solid-color backplate was installed on the experimental table, and the RGB-D camera was kept recording while saving video images (svo2 format). A robot and clamping device were used to grip one end of the pig liver material. The robot was controlled by a teaching controller to stretch horizontally, and the magnitude of the pulling force was read in real time. Based on force measurement experiments and research in real surgical tasks, except for surgeries involving rigid objects such as orthopedics, the operating force of the clamping device during surgery is usually maintained below 2N, and even below 0.5N for delicate operations (such as those involving nerve tissue). Therefore, the experimental pulling force was set to 1N. After the pulling force at the robot's end reached 1±0.05N, the robot's pulling was stopped, and the force sensor reading was kept stable at this level. In the stable state, the clamping device was operated to release the clamps, causing the liver to produce damped oscillations.

[0312] Based on the ZED2 camera's API, a Python program was written to read the image file. First, it reads the image file and generates a depth map (MEASURE.DEPTH) and a point cloud (MEASURE.XYZRGBA). Second, based on actual scale calibration, a distance transformation relationship is established to uniformly convert the recorded pixel coordinates to the actual scale. Third, image HSV (Hue-Saturation-Value) spatial target recognition is used to search for circular markers, providing key parameter adjustment functions (mainly local search radius, HSV tolerance, and minimum connected component area threshold). Finally, contour filtering, ellipse fitting, and center point calculation are performed to achieve frame-by-frame tracking of marker coordinates in the video stream. The experimental images and marker recognition are shown below. Figure 17 As shown.

[0313] The experiment involved capturing images of 10 different pig liver locations. The trajectory curves of the marked points on the projection plane perpendicular to the camera's orientation were recorded as follows: Figure 18 As shown in (a), the curves of different colors represent the trajectories recorded in different groups. By averaging the positions of the 10 groups of markers, the cumulative average trajectory can be obtained, as shown in (a). Figure 18As shown in (b) in the experiment. Based on the dynamic images collected in the experiment and the trajectory of the pig liver markers, the "damped oscillation" phenomenon generated by the actual pig liver after stretching and deformation is obvious, with the following characteristics: (1) The trajectories of the liver deformation in multiple groups are very similar, showing the universality of the damped oscillation behavior; (2) The oscillation amplitude changes greatly with time. The initial amplitude is large, and even the amplitude ratio of some experimental groups exceeds 50% for the first time, showing a large instability overall. However, as time changes, the amplitude decreases rapidly, which is a typical manifestation of viscoelasticity; (3) Comparing the deformation of each group of experiments, it was found that under the same tension (1.0N), the amount of plastic deformation generated by the soft organ is very small. It can be seen that under the relatively small operating force of the endoscope, the soft deformation is mainly dominated by elasticity and viscoelasticity, verifying the correctness of the rigid body term and the soft body term; (4) Based on the quantitative calculation of the position tracking trajectory, the average maximum amplitude ratio of the actual damped oscillation can be preliminarily obtained as 26.1%.

[0314] In addition, the damped oscillation parameters can be determined by analyzing and labeling observable indicators of the liver soft tissue deformation process through professional personnel or large-scale model identification. The labeling includes: a. whether a "damped oscillation" phenomenon occurred; b. the keyframe time at which the soft tissue has entered a stable state from the observation angle; c. the number of soft tissue oscillations from the observation angle. Each set of images can be viewed repeatedly and paused at keyframes, and keyframes considered stable from the observation angle are marked and recorded. Finally, the results of each measurement are averaged and statistically analyzed. The experimental data are shown in Table 7.

[0315] Table 7

[0316]

[0317] Table 7 shows that the first row contains the serial numbers of all participants or large models. For index a, each participant or large model observed damped oscillations in all 10 cases, demonstrating the prevalence of this phenomenon. For index b, most participants or large models observed that the soft body had reached "steady" state and ceased oscillations within 1.1-1.3 seconds after release, with an average of 1.21 seconds. The standard deviation of the data for each group was 0.022, indicating that the time to steady-state was relatively similar for the pig liver experimental case selected in Section 3.6.1. Furthermore, most participants or large models observed 2-4 observable vibrations before reaching observable stability, with an average of 2.98 vibrations. This means that for the 10 cases of the same scale, the damped oscillation behavior was similar, with an observable number of reciprocating motions of 3. This aligns with... Figure 18 The recorded trajectory in (b) is the same. Finally, the average amplitude ratio of the keyframes that the personnel or large model considers to be stable is calculated for the cases that the personnel or large model considers to have entered a stable state. The range was 4.94%, with a maximum of 5.37%. It can be assumed that if the amplitude ratio of the software's reciprocating motion is less than 5%, then from the perspective of personnel or large models, the software has essentially reached a stable state. Personnel should be surgeons with actual experience in laparoscopic and percutaneous endoscopic surgeries, and who clearly understand the concept of damped oscillations and are proficient in related operation and recording methods.

[0318] In practical implementation, a liver model of the same scale (including a tetrahedral model and the patch models required for skinning) can also be built, and simulations can be performed on models based on the optimized dynamic parameter I-XPBD algorithm and two sets of typical XPBD algorithms with fixed parameters. Some keyframe images from the simulation results are shown below. Figure 19 As shown, the indicator length is the maximum amplitude during the oscillation to steady state. This represents the ratio of the maximum amplitude of the initial oscillation. Figure 19 (a) in the image is a keyframe of the I-XPBD algorithm after skinning; Figure 19 (b) in the image is a keyframe displaying the tetrahedral wireframe; the red dots represent the recorded keypoint trajectories. Figure 19 (c) in the figure is a typical XPBD algorithm keyframe and keypoint record with fixed parameters, the length constraint compliance constant is 100, which is the control group I; there is no velocity damping; Figure 19 (d) represents the case with 10% velocity damping, the most typical XPBD algorithm, and serves as control group II. To further compare the real-time simulation with the actual experiment, the trajectory points from the actual experiment were averaged and smoothed before being compared with the three sets of simulated trajectory points, as shown below. Figure 20 As shown, the black curve represents the average trajectory obtained from 10 sets of actual experimental records, the blue curve is the trajectory recorded by the I-XPBD algorithm simulation software, the green trajectory is the trajectory recorded by the typical undamped XPBD algorithm simulation software, and the yellow trajectory is the velocity dissipation coefficient. The trajectory recorded by the typical XPBD algorithm simulation software is 0.1. The specific maximum amplitude ratio and the number of reciprocations with amplitude greater than 5% are shown in Table 8.

[0319] Table 8

[0320]

[0321] Depend on Figure 19 and Figure 20 It can be seen that due to the differences between the geometric coordinate records in the simulation process and the coordinate records of the actual data, the smoothness of the two types of data is inconsistent, but the overall trend is consistent. The software driven by the optimized I-XPBD algorithm will not exhibit the large, unnatural oscillations seen in undamped algorithms, such as... Figure 19 As shown in (c), it will not stabilize as quickly as a typical XPBD algorithm, such as Figure 19As shown in (d) in the diagram, the oscillations begin initially, but stabilize as constraints and damping increase over time. This behavior closely matches the viscoelastic damped oscillation behavior described earlier in this chapter. Specifically, the I-XPBD maximum amplitude ratio is 32.3%, very close to the maximum amplitude of the actual pig liver trajectory at the same scale, differing by only 6.2%. The observable number of reciprocations is also very close to 3 compared to the actual number. Conversely, control group I exhibits excessively strong oscillation behavior, with twice the actual number of reciprocations, while control group II quickly reaches the observable stable stage. In the pig liver case study designed in this chapter, the statistical analysis of key indicators and the comparison of the dynamic deformation process show that, compared to the fixed-parameter XPBD algorithm, the optimized algorithm more closely approximates the actual pig liver situation.

[0322] It should be understood that although the steps in the flowcharts of the above embodiments are shown sequentially according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated in this embodiment, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in the flowcharts of the above embodiments may include multiple steps or multiple stages. These steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least some of the steps or stages in other steps.

[0323] The object deformation display device based on force feedback interaction provided in the embodiments of the present invention will be described below. The object deformation display device based on force feedback interaction described below can be referred to in correspondence with the object deformation display method based on force feedback interaction described above.

[0324] Please see Figure 21 This invention also provides an object deformation display device based on force feedback interaction, comprising:

[0325] The polygon bounding volume construction module 2101 is used to construct a polygon bounding volume using CT data when CT data of an object to be interacted with is received.

[0326] Tetrahedral subdivision module 2102 is used to tetrahedral subdivide a polygon bounding volume to generate a tetrahedral mesh model; the tetrahedral mesh model is composed of multiple tetrahedral elements.

[0327] The particle mass allocation module 2103 is used to allocate the equivalent mass of mesh particles of each tetrahedral unit according to the vertex association relationship between each tetrahedral unit and the unit volume of each tetrahedral unit, so as to obtain a soft mesh model.

[0328] The external force application module 2104 is used to perform external force application actions on the soft mesh model in response to control commands;

[0329] The position update module 2105 is used to iteratively update the particle position of each tetrahedral element based on the equivalent mass of the mesh particles and the force corresponding to the action of the external force applied, so as to obtain the updated particle position.

[0330] The image rendering module 2106 is used to render the image based on the updated position of each particle and output the deformation image corresponding to the object to be interacted with.

[0331] Optionally, the polygon bounding volume building module 2101 is specifically used for:

[0332] When the CT data of the object to be interacted with is received, the initial patch model corresponding to the CT data is obtained;

[0333] Match the corresponding voxel mesh step size according to the scale of the initial patch model, and construct the voxel mesh based on the voxel mesh step size; the voxel mesh consists of multiple mesh elements;

[0334] Select the mesh cells whose center position is within the initial patch model to construct the occupied cell set;

[0335] Extract and smooth the outer contour of the occupied cell set to obtain the initial bounding polygon;

[0336] Offset the initial bounding polygon outwards along the normal direction of each grid cell by the corresponding outward offset of the object to be interacted with, and generate the polygon bounding volume.

[0337] Optionally, the tetrahedral partitioning module 2102 is specifically used for:

[0338] Construct a single initial tetrahedron that covers the bounding volume of the polygon;

[0339] Insert each vertex of the bounding volume of the polygon into the initial tetrahedron in turn, and delete all extra tetrahedrons containing the insertion points of the circumscribed spheres to form a common surface cavity;

[0340] Connect the insertion point to all boundary vertices of the common surface cavity to obtain the intermediate tetrahedron;

[0341] After all the bounding volume vertices have been inserted, the initial tetrahedron is removed to obtain the tetrahedral mesh model.

[0342] Optionally, the particle mass allocation module 2103 is specifically used for:

[0343] Calculate the element volume of each tetrahedral element based on the coordinates of its four vertices.

[0344] If any mesh particle is associated with only a single tetrahedral element, then the equivalent mass of the mesh particle is calculated according to the first equivalent mass calculation formula and assigned to the mesh particle;

[0345] If any mesh particle is associated with at least two tetrahedral elements, the second equivalent mass calculation formula calculates the equivalent mass of the mesh particle and assigns it to the mesh particle;

[0346] Once all mesh particles are assigned equivalent mass, a soft mesh model is obtained.

[0347] Optionally, the device also includes a model evaluation module, specifically used for:

[0348] After obtaining the soft mesh model, traverse all tetrahedral elements within the soft mesh model and calculate the nearest neighbor distance probability density distribution index and dihedral angle quality index one by one.

[0349] If there are unqualified tetrahedral elements whose nearest neighbor distance probability density distribution index and / or dihedral angle quality index do not meet the preset standards, and the number of unqualified tetrahedral elements is less than the preset element threshold, then the region to be adjusted where the unqualified tetrahedral elements are located is located, and the polygon bounding volume corresponding to the region to be adjusted is updated according to the new voxel mesh step size.

[0350] Jump to the step of performing tetrahedral subdivision of the polygon bounding volume to generate a tetrahedral mesh model;

[0351] If there are unqualified tetrahedral elements whose nearest neighbor distance probability density distribution index and / or dihedral angle quality index do not meet the preset standards, and the number of unqualified tetrahedral elements is not less than the preset element threshold, then obtain a new voxel mesh step size and a new outward offset, and jump to the step of constructing a polygon bounding volume using CT data.

[0352] If the nearest neighbor distance probability density distribution index and dihedral angle quality index of all tetrahedral elements meet the preset standards, then the response control command is executed to apply external force to the soft mesh model.

[0353] Optionally, the location update module 2105 is specifically used for:

[0354] The current velocity of each tetrahedral element is updated based on the force corresponding to the action of the external force application, and the dynamic velocity dissipation coefficient is calculated in combination with the duration of soft deformation.

[0355] The initial velocity of each vertex is attenuated and corrected using the dynamic velocity dissipation coefficients to obtain the particle corrected velocity;

[0356] The initial predicted position of each particle in the current frame is predicted based on the corrected velocity of each particle.

[0357] The dynamic compliance constant is calculated based on the duration of soft deformation at each simulation time step, and the Lagrange operator of the length constraint is adjusted based on the dynamic compliance constant to obtain the updated length constraint;

[0358] The initial predicted position is iteratively corrected under the constraints of length update, mention constraint, and collision penetration prevention constraint to obtain the updated position of the particle and update the current velocity of the particle.

[0359] Optionally, the formula for calculating the dynamic velocity dissipation coefficient is:

[0360] ;

[0361] in, The dynamic velocity dissipation coefficient, The velocity dissipation time-varying coefficient, t is the preset base dissipation coefficient, and t is the duration of the soft deformation.

[0362] Optionally, the formula for calculating the dynamic compliance constant is:

[0363] ;

[0364] Where k, b, and a are preset KV model constants. Corrects the particle velocity.

[0365] Optionally, the image rendering module 2106 is specifically used for:

[0366] The position of each particle is updated and transformed to screen pixel coordinates;

[0367] Rasterize the triangular facet formed by the vertices of the particle update position to obtain multiple pixels and associate them with the screen pixel coordinates;

[0368] Each pixel is assigned brightness and / or color to generate and output the deformed image corresponding to the object to be interacted with.

[0369] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working process of the above-described device and module can be referred to the corresponding process in the foregoing method embodiments, and will not be repeated here.

[0370] In the several embodiments provided by this invention, it should be understood that the disclosed apparatus and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of modules is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple modules or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between apparatuses or modules may be electrical, mechanical, or other forms.

[0371] The modules described as separate components may or may not be physically separate. The components shown as modules may or may not be physical modules; that is, they may be located in one place or distributed across multiple network modules. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs.

[0372] Furthermore, the functional modules in the various embodiments of the present invention can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module. The integrated modules described above can be implemented in hardware or as software functional modules.

[0373] The above-described embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for displaying deformation of an object based on force feedback interaction, characterized in that, include: When the CT data of the object to be interacted with is received, a polygon bounding volume is constructed using the CT data; The polygonal bounding volume is tetrahedralized to generate a tetrahedral mesh model; The tetrahedral mesh model is composed of multiple tetrahedral elements; Based on the vertex relationships between the tetrahedral elements and the element volume of each tetrahedral element, the equivalent mass of the mesh particles of each tetrahedral element is allocated to obtain the soft mesh model. The system responds to control commands by applying external forces to the soft mesh model. Based on the equivalent mass of the mesh particles and the force corresponding to the applied external force, the particle position of each tetrahedral unit is iteratively updated to obtain the updated particle position. Based on the updated positions of each particle, the screen is rendered, and the deformed screen corresponding to the object to be interacted with is output.

2. The method of claim 1, wherein, When the CT data of the object to be interacted with is received, constructing a polygon bounding volume using the CT data includes: When the CT data of the object to be interacted with is received, the initial patch model corresponding to the CT data is obtained; Match the corresponding voxel mesh step size according to the scale of the initial patch model, and construct a voxel mesh based on the voxel mesh step size; the voxel mesh includes multiple mesh elements; Select mesh cells whose center position is within the initial patch model to construct an occupied cell set; The initial bounding polygon is obtained by extracting and smoothing the outer contour of the occupied unit set; The initial bounding polygon is offset outward along the normal direction of each of the mesh cells by the outward offset corresponding to the object to be interacted with, thereby generating a polygon bounding volume.

3. The object deformation display method based on force feedback interaction according to claim 1, characterized in that, The step of tetrahedralizing the polygonal bounding volume to generate a tetrahedral mesh model includes: Construct a single initial tetrahedron covering the bounding volume of the polygon; Each vertex of the polygonal bounding body is inserted into the initial tetrahedron in sequence, and all extra tetrahedrons containing the insertion points on the circumscribed spheres are deleted to form a common surface cavity; Connect the insertion point with all the boundary vertices of the common surface cavity to obtain the intermediate tetrahedron; After all the bounding volume vertices have been inserted, the initial tetrahedron is removed to obtain a tetrahedral mesh model.

4. The object deformation display method based on force feedback interaction according to claim 1, characterized in that, The process of allocating the equivalent mass of mesh particles for each tetrahedral element based on the vertex relationships between the tetrahedral elements and the element volume of each tetrahedral element to obtain a soft mesh model includes: Calculate the unit volume of each tetrahedral unit according to the coordinates of its four vertices; If any mesh particle is associated with only a single tetrahedral unit, then the equivalent mass of the mesh particle is calculated according to the first equivalent mass calculation formula and allocated to the mesh particle; If any mesh particle is associated with at least two of the tetrahedral elements, the second equivalent mass calculation formula calculates the equivalent mass of the mesh particle and assigns it to the mesh particle; Once all the mesh particles are assigned equivalent mass, a soft mesh model is obtained.

5. The object deformation display method based on force feedback interaction according to any one of claims 1 to 4, characterized in that, Also includes: After obtaining the soft mesh model, all tetrahedral elements within the soft mesh model are traversed, and the nearest neighbor distance probability density distribution index and dihedral angle quality index are calculated one by one. If there are unqualified tetrahedral units that do not meet the preset standards for the nearest neighbor distance probability density distribution index and / or the dihedral angle quality index, and the number of unqualified tetrahedral units is less than the preset unit threshold, then the region to be adjusted where the unqualified tetrahedral unit is located is located, and the polygon bounding volume corresponding to the region to be adjusted is updated according to the new voxel mesh step size. Jump to the step of performing tetrahedral subdivision of the polygon bounding volume to generate a tetrahedral mesh model; If there are unqualified tetrahedral cells that do not meet the preset standards for the nearest neighbor distance probability density distribution index and / or the dihedral angle quality index, and the number of unqualified tetrahedral cells is not less than the preset cell threshold, then a new voxel mesh step size and a new outward offset are obtained, and the process jumps to the step of constructing a polygon bounding volume using the CT data. If the nearest neighbor distance probability density distribution index and dihedral angle quality index of all tetrahedral elements meet the preset compliance conditions, then the response control command is executed to apply external force to the soft mesh model.

6. The object deformation display method based on force feedback interaction according to claim 1, characterized in that, The step of iteratively updating the particle position of each tetrahedral unit based on the equivalent mass of the mesh particles and the force corresponding to the applied external force, to obtain the updated particle position, includes: The current velocity of each particle in the tetrahedral unit is updated based on the force corresponding to the external force application action, and the dynamic velocity dissipation coefficient is calculated in combination with the duration of soft deformation. The initial velocity of each vertex is attenuated and corrected using the dynamic velocity dissipation coefficients described above to obtain the particle corrected velocity. The initial predicted position of each particle in the current frame is predicted based on the particle correction velocity. At each simulation time step, a dynamic compliance constant is calculated based on the duration of soft deformation, and the Lagrange operator of the length constraint is adjusted based on the dynamic compliance constant to obtain the updated length constraint. The initial predicted position is iteratively corrected under the updated length constraint, mention constraint, and collision penetration prevention constraint to obtain the updated particle position and update the current velocity of the particle.

7. The object deformation display method based on force feedback interaction according to claim 6, characterized in that, The formula for calculating the dynamic velocity dissipation coefficient is as follows: ； in, The dynamic velocity dissipation coefficient, The velocity dissipation time-varying coefficient, t is the preset base dissipation coefficient, and t is the duration of the soft deformation.

8. The object deformation display method based on force feedback interaction according to claim 6, characterized in that, The formula for calculating the dynamic compliance constant is as follows: ； Where k, b, and a are preset KV model constants. Corrects the particle velocity.

9. The object deformation display method based on force feedback interaction according to claim 1, characterized in that, The step of rendering the image based on the updated positions of each particle and outputting the deformed image corresponding to the object to be interacted with includes: The position of each particle is updated and transformed to screen pixel coordinates; The triangular facet formed by the vertices to which the particle update position belongs is rasterized to obtain multiple pixels and associated with the screen pixel coordinates; Each pixel is assigned brightness and / or color to generate and output the deformed image corresponding to the object to be interacted with.

10. A force feedback interaction-based object deformation display device, characterized in that, include: A polygon bounding volume construction module is used to construct a polygon bounding volume using the CT data of the object to be interacted with when the CT data of the object is received. The tetrahedral subdivision module is used to perform tetrahedral subdivision on the polygonal bounding volume to generate a tetrahedral mesh model; the tetrahedral mesh model is composed of multiple tetrahedral elements; The particle mass allocation module is used to allocate the equivalent mass of the mesh particles of each tetrahedral unit according to the vertex association relationship between each tetrahedral unit and the unit volume of each tetrahedral unit, so as to obtain a soft mesh model. An external force application module is used to perform external force application actions on the soft mesh model in response to control commands; The position update module is used to iteratively update the particle position of each tetrahedral unit based on the equivalent mass of the mesh particle and the force corresponding to the action of the external force applied, so as to obtain the particle update position; The image rendering module is used to render the image based on the updated position of each particle and output the deformation image corresponding to the object to be interacted with.

Images