WO2023000560A1

WO2023000560A1 - Reduction trajectory automatic planning method for parallel fracture surgical robot

Info

Publication number: WO2023000560A1
Application number: PCT/CN2021/131067
Authority: WO
Inventors: 孙涛; 李锦龙; 刘传耙; 连宾宾; 宋轶民
Original assignee: 天津大学
Priority date: 2021-07-21
Filing date: 2021-11-17
Publication date: 2023-01-26
Also published as: CN113781495B; CN113781495A

Abstract

Disclosed in the present invention is a reduction trajectory automatic planning method for a parallel fracture surgical robot. The method comprises: on the basis of a CT image of a patient, obtaining a distal bone model, a proximal bone model and an intact-side bone model by means of three-dimensional reconstruction; acquiring a fracture reduction goal by taking the intact-side bone model of the patient as a reduction reference; defining a collision detection threshold, and searching for the minimum distance point between fractured bones on the basis of an octree search algorithm, so as to realize fractured bone collision detection; establishing a personalized fracture muscle model on the basis of an OpenSim standard model, and acquiring a muscle traction force in the reduction process; and designing a trajectory search node and an evaluation function of an A* algorithm by taking no collision, the minimum muscle traction force and the shortest path as goals, thereby realizing trajectory planning. By means of the method in the present invention, the collision between fractured bone blocks is prevented, the excessive traction of a muscle is reduced, and a fracture reduction path is shortened on the basis of an improved A* search algorithm, thereby improving a fracture reduction effect.

Description

Automatic reset trajectory planning method for parallel fracture surgery robot

technical field

The invention relates to an automatic reset trajectory planning method, in particular to an automatic reset trajectory planning method for a parallel fracture surgery robot.

Background technique

Traditional fracture reduction surgery needs to expose the bone tissue through a large incision, and the doctor restores the anatomical position of the fracture under direct vision. Traditional reduction methods are limited by physician experience and intraoperative equipment, and there are risks such as large trauma, susceptibility to infection, and secondary fractures.

With the integration of robotics, computer information technology and orthopedic medicine, fracture reduction surgery based on parallel robots is generally considered to be an advantageous solution to achieve precise and safe fracture reduction. With its advantages of minimal invasiveness and high precision, parallel robots can effectively solve the defects of traditional reduction surgery.

Parallel robot-assisted fracture reduction means that the doctor uses computer-aided software to plan the reset trajectory of the broken bone, and based on the robot kinematics algorithm, maps the reset trajectory of the broken bone to the robot's motion trajectory, and the robot executes the above trajectory to achieve the purpose of fracture reduction . In the above process, fracture reduction trajectory planning is a key link. Its purpose is to find an optimal trajectory, avoid collisions between bone fragments during the reduction process, reduce damage to soft tissues such as muscles, and enable the fractured bone fragments to be safely and accurately to its anatomical location.

The existing fracture reduction trajectory planning methods can be divided into two categories: interactive and automatic. The interactive trajectory planning method means that the doctor uses the mouse or keyboard to adjust the position and posture of the broken bone fragments, and artificially designate the reset trajectory of the bone fragments. This type of method is mainly implemented by doctors. Due to the different experience of doctors, the generated trajectories are very different. In the automatic trajectory planning method, the computer uses relevant algorithms to automatically generate the reset trajectory of the bone block. The existing automatic trajectory planning methods have problems such as low collision detection accuracy, long reset path, and large muscle pulling force, and cannot obtain safe and effective fractures. Reset track.

In summary, the existing fracture reduction trajectory planning methods are difficult to meet the urgent needs of fracture clinical surgery.

Contents of the invention

The purpose of the present invention is to overcome the shortcomings of the prior art, to provide an automatic planning method for the reset trajectory of a parallel fracture surgery robot, to avoid the collision of bone blocks during the reset operation, to reduce the damage to soft tissues such as muscles, and to effectively improve fracture reduction. accuracy and efficiency.

To achieve the above object, the present invention adopts the following technical solutions:

The present invention is oriented towards an automatic reset trajectory planning method for a parallel fracture surgery robot, comprising the following steps:

(1) Obtain the three-dimensional fracture model of the patient, the specific steps are:

(1a) scan and obtain the CT data of the fracture patient, utilize the ITK tool kit to segment the bone fragments from the CT data, and the described bone fragments include distal bone, proximal bone and healthy side bone;

(1b) Use the moving cube algorithm encapsulated in the VTK toolkit to perform 3D reconstruction of the segmented bone block, and store the reconstruction result as a binary STL mesh model to obtain the distal bone model, proximal bone model and uninjured bone model ;Each model surface is composed of several triangle faces;

(2) Define the initial pose and target pose of the distal bone reset, and provide the initial and target nodes for subsequent trajectory planning. The specific steps are:

(2a) Use VTK to build a virtual 3D scene, and import the distal bone model, proximal bone model, and uninjured bone model respectively. In the initial import state, the model coordinate systems of each model coincide with the world coordinate system of the 3D virtual scene;

(2b) Define the model coordinate system of the distal bone model as O _i -xi y _i z _i , _i =1,2,...,t, and use the coordinate system O _i _{-xi y i z i} _relative _to the world The pose of the coordinate system represents the pose T _i of the distal bone model, and the change in the value of i corresponds to the change of the pose of the distal bone model; when i=1, the model coordinate system of the distal bone model coincides with the world coordinate system , that is, the pose of the distal bone model in the initial state is the identity matrix T ₁ ; when i=t, the pose T _t of the distal bone model is the desired target pose;

(2c) Taking the XZ plane in the three-dimensional virtual scene as a reference plane, performing mirror transformation on the bone model of the healthy side;

(2d) fix the proximal bone model, use the ICP algorithm to obtain the transformation matrix between the proximal bone model and the mirrored healthy side bone model, so that the proximal bone model and the mirrored healthy side bone model overlap;

(2e) Fix the mirrored bone model of the healthy side, and use the ICP algorithm to obtain the transformation matrix between the distal bone model and the mirrored healthy side bone model, so that the mirrored healthy side bone model and the distal bone model overlap; at this time The distal bone model reaches the reset state, read the pose T _t of the distal bone model, and T _t is the target pose;

(3) Perform interference analysis on the distal bone model and the proximal bone model during the fracture reduction process, the specific steps are:

(3a) respectively import the distal bone model and the proximal bone model into the three-dimensional virtual scene built by VTK;

(3b) Based on the octree search algorithm, calculate the shortest distance l between the proximal bone model and the distal bone model;

(3c) Obtain the triangular faces that constitute each model by intercepting the voxels of the CT image, and then calculate the detection threshold l' for the collision between the distal bone model and the proximal bone model during the reset process;

(3d) Judging whether there is a risk of collision between the distal bone model and the proximal bone model based on the distance l, if l>l', the distal bone and the proximal bone will not collide; if l≤l', the distal bone model There is a risk of collision with the proximal bone;

(4) Utilize the standard musculoskeletal model provided by OpenSim software to obtain the pulling force of each muscle in the reset process, and the specific steps are:

(4a) Calculate the scaling factor s according to the topological relationship between the OpenSim standard bone model and the anatomical landmarks of the patient's bone model after the reset, and then use the scaling factor s to scale the construction parameters of the standard bone model to obtain a bone model that is similar in size to the patient's The standard skeletal muscle model;

(4b) Use the ICP algorithm to obtain the conversion matrix between the patient's bone model after step (2e) and the standard skeletal muscle model obtained through step (4a), and then replace the standard skeletal muscle model with the reset one Patient bone model, and set the fixed connection relationship between the muscle attachment point and the distal bone model and the proximal bone model. If the muscle attachment point is close to the distal bone model, it is stipulated that the muscle attachment point is fixedly connected with the distal bone model; On the contrary, it is stipulated that the muscle attachment point is fixedly connected with the proximal bone model;

(4c) First, use the C++ API interface of _OpenSim software to establish a six-degree-of-freedom joint between the distal bone model and the proximal bone model in the three-dimensional virtual scene. The six-degree-of-freedom joint has _t y _t z _t coordinate axis x _t , y _t , z _t three translation degrees of freedom and three rotation degrees of freedom around the coordinate axis x _t , y _t , z _t ; secondly, establish the distal bone model along the coordinate axis The translational distance d _x , d _y , d _z of x _t , y _t , z _t and the angle α, β, γ of the rotation of the distal bone model around the coordinate axes x _t , y _t , z _t and the position of the distal bone model The relational expression of pose T _i ; finally, take T ₁ as the pose of the distal bone model in the initial state, and T _t as the target pose, to simulate the movement of the distal bone model relative to the proximal bone model; if the distal bone model is known T _i of the bone model, the joint parameters d _x , d _y , d _z , α, β, γ can be obtained; the joint parameters d _x , d _y , d _z , α, β, γ are used as the input of the OpenSim software , then the OpenSim software outputs the pulling force of each muscle corresponding to the pose T _i ;

(5) On the basis of the original execution process of the A* algorithm, the collision detection method in step (3) and the muscle force analysis method in step (4) are combined, and the trajectory nodes and evaluation functions of the A* algorithm are redesigned, Perform fracture reduction trajectory planning, the specific steps are:

(5a) Define the trajectory node N _i (O _xi ,O _yi ,O _zi ,θ _i )(i=1,2,···,t) corresponding to the pose T _i of the distal bone model one-to-one, the trajectory node The components O _xi , O _yi , and O _zi of N _i are the x, y, and z components of the origin O _i of the coordinate system O _i -xi y _i z _i respectively _, and the component θ _i indicates that the distal bone model revolves around the The angle by which the axis n rotates counterclockwise;

(5b) Combining the constraints of no collision, the shortest path, and the smallest muscle pulling force on the fracture reduction trajectory, the evaluation function of the A* algorithm is defined as:

Among them, g(N _i ) is the cost of moving the origin of the distal bone model coordinate system from the initial node N ₁ to node N _i ; h(N _i ) indicates that the origin of the distal bone model coordinate system moves from node N _i to the target node N The cost of _t ; r(N _i )=c ₁ |θ _i -θ _t | is the attitude penalty function, which represents the cost required for the distal bone model coordinate system to rotate from the current node N _i to the target node N _t , and c ₁ is Correlation penalty coefficient, the value of c ₁ should make r(N _i ) have the same order of magnitude as g(N _i ), h(N _i ), m(N _i ); m(N _i )=c ₂ (F _i -F ₁ ) is the muscle force penalty function, F _i (i=1,2,...,t) represents the algebraic sum of the muscle forces when the distal bone model is located at the trajectory node N _i , where F _i (i= 1,2,···,t) use the method in step (5a) to map the node N _i to the pose T _i of the distal bone model, and then use the method in step (4c) to obtain the muscle force Algebraic sum; c ₂ is the correlation penalty coefficient, the value of c ₂ should make m(N _i ) have the same order of magnitude as g(N _i ), h(N _i ), r(N _i );

(5c) Based on the improved A* search algorithm, the trajectory nodes N ₁ (0,0,0,0) and N _t (O _xt ,O _xt ,O _xt ,θ _t ) set in step (5a) are used as the trajectory respectively planning start node and end node, and optimize the reset trajectory node according to the evaluation function f(N _i ) set in step (5b), in the optimization process, use the method in step (3) to select all neighbors of N _c The domain nodes perform collision detection, and the reset trajectory node sequence N ₂ ,···,N _t-1 is optimized, and then the trajectory node N _i (O _xi ,O _yi ,O _zi ,θ _i ) and the distal bone pose T According to the corresponding relationship between _i , the optimized reset trajectory node sequence N ₂ ,···,N _t-1 is converted into the pose of the distal bone model to complete the fracture reduction trajectory planning.

The beneficial effects of the present invention are: the use of the collision detection method avoids the collision between the fractured bone fragments in the reset process, reduces the excessive pulling of the muscles through the OpenSim muscle force simulation, and shortens the fracture reduction path based on the improved A* search algorithm And improve the efficiency of trajectory planning.

Description of drawings

Fig. 1 is the overall flow chart of the automatic planning method for the reset trajectory of the parallel fracture surgery robot of the present invention;

Fig. 2 is a flow chart of fracture model acquisition;

Figure 3 is a flow chart of the definition of the fracture initial and target poses;

Fig. 4 is a schematic diagram of fracture reduction target extraction;

Fig. 5 is a collision detection flowchart;

Fig. 6 is a schematic diagram of CT image voxel size;

Fig. 7 is a schematic diagram of the definition of the collision detection threshold;

Fig. 8 is a schematic diagram of the OpenSim standard skeletal muscle model;

Fig. 9 is a flow chart of muscle force analysis;

Figure 10a is a schematic diagram of the patient's bone model and the standard bone model before replacement;

Figure 10b is a schematic diagram of the patient's bone model and the standard bone model after replacement;

Fig. 11 is a flow chart of fracture reduction track node optimization.

detailed description

The specific implementation manners of the present invention will be described in detail below in conjunction with the accompanying drawings.

Referring to the accompanying drawings, the automatic reset trajectory planning method for parallel fracture surgical robots of the present invention includes the following steps:

(1) Obtain the patient's three-dimensional fracture model and provide a data source for the subsequent steps. The specific steps are:

(1a) Scan and obtain the CT data of the fracture patient, use the ITK (Insight Segmentation and Registration Toolkit) toolkit to segment the bone fragments from the CT data, and the bone fragments include the distal bone (the proximal end of the broken bone), the proximal The end bone (the distal end of the broken bone) and the uninjured bone (the side where the fracture did not occur).

(1b) Use the moving cube algorithm encapsulated in the VTK toolkit to perform 3D reconstruction of the segmented bone block, and store the reconstruction result as a binary STL mesh model to obtain the distal bone model, proximal bone model and uninjured bone model ;Each model surface is composed of several triangle faces.

(2a) Create a virtual 3D scene using VTK, and import the distal bone model, proximal bone model and uninjured bone model respectively. In the initial import state, the model coordinate system of each model coincides with the world coordinate system of the 3D virtual scene.

(2b) Define the model coordinate system of the distal bone model as O _i -xi y _i z _i ( _i =1,2,···,t), using the coordinate system O _i _{-xi y i z i} _relative _to The pose of the world coordinate system represents the pose T _i of the distal bone model, and changes in the value of i correspond to changes in the pose of the distal bone model. When i=1, the model coordinate system of the distal bone model coincides with the world coordinate system, that is, the pose of the distal bone model in the initial state is the identity matrix T ₁ ; when i=t, the pose of the distal bone model T _t is the desired target pose, and the calculation of T _t can be found in the following steps.

(2c) Taking the XZ plane in the three-dimensional virtual scene as a reference plane, performing mirror transformation on the bone model of the healthy side.

(2d) Fix the proximal bone model, and use the ICP algorithm to obtain the transformation matrix between the proximal bone model and the mirrored healthy side bone model, so that the proximal bone model and the mirrored healthy side bone model coincide.

(2e) Fix the mirrored healthy side bone model, and use the ICP algorithm to obtain the transformation matrix between the distal bone model and the mirrored healthy side bone model, so that the mirrored healthy side bone model and the distal bone model overlap. At this time, the distal bone model reaches the reset state, and the pose T _t of the distal bone model is read, and T _t is the target pose.

(3) Perform interference analysis on the distal bone model and the proximal bone model during the fracture reduction process to provide constraints for subsequent trajectory planning and avoid collisions between the distal bone model and the proximal bone model during the reduction process. The specific steps are:

(3a) Import the distal bone model and the proximal bone model respectively into the 3D virtual scene built by VTK.

(3b) Calculate the shortest distance l between the proximal bone model and the distal bone model based on an octree search algorithm.

The calculation process of the detection threshold l' is as follows:

The size of the voxel is determined by the CT scanning accuracy. Assuming that the CT scanning accuracy is a×b×c, as shown in Figure 6, the maximum side length s _max of the triangular surface can be obtained as:

Construct an equilateral triangle with the maximum side length s _max , as shown in Figure 7, in this case define the collision detection threshold l' as:

The accuracy of CT scanning is generally 1mm×1mm×0.625mm, and the threshold of collision occurrence is 0.893mm.

(3d) Judging whether there is a risk of collision between the distal bone model and the proximal bone model based on the distance l, the specific instructions are as follows:

The surface of the three-dimensional model of the distal bone and the proximal bone is composed of several triangular faces. If the bone blocks collide and the distal bone and the proximal bone interfere, there must be two or more intersecting triangular faces. If l>l', there will be no collision between the distal bone and the proximal bone; if l≤l', there is a risk of collision between the distal bone and the proximal bone.

(4) Use the standard musculoskeletal model provided by OpenSim software to obtain the pulling force of each muscle during the reset process, and use it as a constraint for subsequent trajectory planning to avoid damage to muscles and other soft tissues during the reset process. The specific steps are:

(4a) According to the topological relationship between the OpenSim standard bone model (as shown in Figure 8) and the anatomical landmarks of the patient's bone model after the reset, calculate the scaling factor s, and then use the scaling factor s to scale the construction parameters of the standard bone model, Obtain a standard musculoskeletal model that is similar in size to the patient's bone model.

The construction parameters include geometry, mass, center of mass, joint position, muscle attachment points, and muscle parameters, etc. The scaling factor is calculated as follows:

Among them, L _i and L' _i represent the distance between the individualized bone model and a certain pair of anatomical landmarks (such as the tip of the inner and outer malleolus) in the standard bone model after the reset is completed, and k represents the group number of the selected anatomical landmarks .

(4b) Use the ICP algorithm to obtain the conversion matrix between the patient's bone model after step (2e) has been reset (see Figure 10a and the standard skeletal muscle model obtained through step (4a), and then replace the standard skeletal muscle model with The bone model of the patient after the reset is completed, as shown in Figure 10b. And set the fixed relationship between the muscle attachment point and the distal bone model and the proximal bone model. If the muscle attachment point is close to the distal bone model, the muscle attachment point is specified. The attachment point is fixedly connected with the distal bone model; otherwise, the muscle attachment point is fixedly connected with the proximal bone model.

(4c) First, use the C++ API interface of the OpenSim software to establish a six-degree-of-freedom joint between the distal bone model and the proximal bone model in a three-dimensional virtual scene, as shown in Figure 10b, the six-degree-of-freedom joint It has three translation degrees of freedom along the O _t -x _t y _tz _t coordinate axes x _t , y _t , z _t and three rotation degrees of freedom around the coordinate axes x _t , y _t , z _t . Secondly, establish the translational distance d _x , d _y , d _z of the distal bone model along the coordinate axes _{x t} _, y _t _, _{z t} _and the angle α, The relationship between β, γ and the pose T _i of the distal bone model. Finally, taking T ₁ as the pose of the distal bone model in the initial state, and T _t as the target pose, the movement of the distal bone model relative to the proximal bone model is simulated. If the pose T _i of the distal bone model is known, the joint parameters d _x , d _y , d _z , α, β, γ can be obtained; the joint parameters d _x , d _y , d _z , α, β, γ As the input of the OpenSim software, the OpenSim software outputs the pulling force of each muscle corresponding to the pose T _i .

d _x , d _y , d _z , α, β, γ have the following relationship with the distal bone pose T _i :

Among them, ^t T _i is the pose matrix of the coordinate system O _i _{-xi y i z i} _relative _to the coordinate system O _t -x _t y _t z _t , ^t R _i is the pose matrix of ^t T _i , and ^t d _i is The displacement vector of ^t T _i , x=[1 0 0] ^T , y=[0 1 0] ^T , z=[0 0 1] ^T . (5) On the basis of the original execution process of the A* algorithm, the collision detection method of step (3) and the muscle force analysis method of step (4) are combined, and the trajectory nodes and evaluation functions of the A* algorithm are redesigned, Perform fracture reduction trajectory planning, the specific steps are:

(5a) Define the trajectory node N _i (O _xi ,O _yi ,O _zi ,θ _i )(i=1,2,···,t) corresponding to the pose T _i of the distal bone model one-to-one, the trajectory node The components O _xi , O _yi , and O _zi of N _i are the x, y, and z components of the origin O _i of the coordinate system O _i -xi y _i z _i respectively _, and the component θ _i indicates that the distal bone model revolves around the The angle by which the axis n is rotated counterclockwise.

The process of deriving the distal bone pose T _i from the trajectory nodes N _i (O _xi , O _yi , O _zi , θ _i ) is as follows:

Since the axis n is the rotation axis of the distal bone model, the axis n remains unchanged during the rotation of the distal bone model from the initial pose to the target pose, so the direction of the axis n can be calculated by the following formula:

R _t n=n

In the formula, R _t is the attitude matrix of T _t . The component θ _i of the trajectory node N _i is calculated as follows:

Knowing the direction of the rotation axis n and the rotation angle θ _i , it can be obtained from the Rodrigues formula:

R _i ＝cosθ _i I+(1-cosθ _i )nn ^T +sinθ _i n^

In the formula, R _i is the attitude matrix of matrix T _i , n^ represents the anti-symmetric matrix of vector n, and I is the identity matrix.

Therefore, each component of the trajectory node N _i is known, and the distal bone pose can be obtained:

Among them, g(N _i ) is the cost of moving the origin of the distal bone model coordinate system from the initial node N ₁ to node N _i ; h(N _i ) indicates that the origin of the distal bone model coordinate system moves from node N _i to the target node N The cost of _t ; r(N _i )=c ₁ |θ _i -θ _t | is the attitude penalty function, which represents the cost required for the distal bone model coordinate system to rotate from the current node N _i to the target node N _t , and c ₁ is Correlation penalty coefficient, the value of c ₁ should make r(N _i ) have the same order of magnitude as g(N _i ), h(N _i ), m(N _i ); m(N _i )=c ₂ (F _i -F ₁ ) is the muscle force penalty function, F _i (i=1,2,...,t) represents the algebraic sum of the muscle forces when the distal bone model is located at the trajectory node N _i , where F _i (i= 1,2,···,t) use the method in step (5a) to map the node N _i to the pose T _i of the distal bone model, and then use the method in step (4c) to obtain the muscle force Algebraic sum. c ₂ is the correlation penalty coefficient, the value of c ₂ should make m(N _i ) have the same order of magnitude as g(N _i ), h(N _i ), r(N _i ).

As shown in Figure 11, the specific process of optimizing fracture reduction trajectory nodes is as follows:

Step1: Define empty tables open and close, and insert the start node N ₁ of the reset trajectory into the open table.

Step2: Determine whether the open table is empty, if open is empty, end the search; if it is not empty, find the node N _c (O _xc ,O _xc ,O _xc ,θ _c ) with the smallest value of the evaluation function in the open table, record as the current node, then delete the node from the open table and insert it into the close table.

Step3: The neighbor node of the current node N _c in the close table is N _c+1 , and then use the method in step (5a) to map N _c+1 to the pose of the distal bone model, and use the method in step (3) The method performs collision detection on all neighborhood nodes of N _c , and defines the collection of non-collided nodes as a child table;

The relationship between the current node _Nc and the neighbor node _Nc ₊₁ of the current node Nc is as follows:

Step4: Traverse each node N _j in the child table, and calculate the new evaluation function value of node N _j when N _c is the parent node; if the open table contains node N _j , and the new evaluation function value of node N _j is smaller than the old value value of the evaluation function, then set the parent node of N _j to N _c ; if the open table does not contain node N _j , set the parent node of N _j to N _c , and add N _j to the open table.

Step5: Determine whether the end node N _t is included in the open table, if so, terminate the search process, and execute Step6, otherwise return to Step2 and execute again.

Step6: When the target node N _t is searched, it means that the trajectory planning is completed. Starting from the terminal node N _t , the parent node is continuously searched upward until the parent node is the starting node N _1. N ₁ , N ₂ , ... N _t together constitute the trajectory node of fracture reduction, and the above nodes are transformed into The pose of the distal bone model can be used to complete the fracture reduction trajectory planning.

The invention utilizes the collision detection technology to avoid the collision between broken bone blocks in the reset process, reduces the excessive pulling of the muscles through the OpenSim muscle force simulation, shortens the fracture reset path and improves the trajectory planning based on the improved A* search algorithm s efficiency. Therefore, an automatic reset trajectory planning method for a parallel fracture surgery robot of the present invention can meet the needs of clinical fracture reduction surgery.

Claims

The reset trajectory automatic planning method for parallel fracture surgery robot is characterized in that it comprises the following steps:

(1) Obtain the three-dimensional fracture model of the patient, the specific steps are:

(1a) scan and obtain the CT data of the fracture patient, utilize the ITK tool kit to segment the bone fragments from the CT data, and the described bone fragments include distal bone, proximal bone and healthy side bone;

(1b) Use the moving cube algorithm encapsulated in the VTK toolkit to perform 3D reconstruction of the segmented bone block, and store the reconstruction result as a binary STL mesh model to obtain the distal bone model, proximal bone model and uninjured bone model ;Each model surface is composed of several triangle faces;

(2) Define the initial pose and target pose of the distal bone reset, and provide the initial and target nodes for subsequent trajectory planning. The specific steps are:

(2a) Use VTK to build a virtual 3D scene, and import the distal bone model, proximal bone model, and uninjured bone model respectively. In the initial import state, the model coordinate systems of each model coincide with the world coordinate system of the 3D virtual scene;

(2b) Define the model coordinate system of the distal bone model as O i -xi y i z i , i =1,2,...,t, use the coordinate system O i -xi y i z i relative to the world coordinate system The pose of represents the pose T i of the distal bone model, and the change of the value of i corresponds to the change of the pose of the distal bone model; when i=1, the model coordinate system of the distal bone model coincides with the world coordinate system, that is In the initial state, the pose of the distal bone model is the identity matrix T 1 ; when i=t, the pose T t of the distal bone model is the desired target pose;

(2c) Taking the XZ plane in the three-dimensional virtual scene as a reference plane, performing mirror transformation on the bone model of the healthy side;

(2d) fix the proximal bone model, use the ICP algorithm to obtain the transformation matrix between the proximal bone model and the mirrored healthy side bone model, so that the proximal bone model and the mirrored healthy side bone model overlap;

(2e) Fix the mirrored bone model of the healthy side, and use the ICP algorithm to obtain the transformation matrix between the distal bone model and the mirrored healthy side bone model, so that the mirrored healthy side bone model and the distal bone model overlap; at this time The distal bone model reaches the reset state, read the pose T t of the distal bone model, and T t is the target pose;

(3) Perform interference analysis on the distal bone model and the proximal bone model during the fracture reduction process, the specific steps are:

(3a) respectively import the distal bone model and the proximal bone model into the three-dimensional virtual scene built by VTK;

(3b) Based on the octree search algorithm, calculate the shortest distance l between the proximal bone model and the distal bone model;

(3c) Obtain the triangular faces that constitute each model by intercepting the voxels of the CT image, and then calculate the detection threshold l' for the collision between the distal bone model and the proximal bone model during the reset process;

(3d) Judging whether there is a risk of collision between the distal bone model and the proximal bone model based on the distance l, if l>l', the distal bone and the proximal bone will not collide; if l≤l', the distal bone model There is a risk of collision with the proximal bone;

(4) Use the standard musculoskeletal model provided by OpenSim software to obtain the pulling force of each muscle during the reset process. The specific steps are:

(4a) Calculate the scaling factor s according to the topological relationship between the OpenSim standard bone model and the anatomical landmarks of the patient's bone model after reset, and then use the scaling factor s to scale the construction parameters of the standard bone model to obtain a bone model that is similar in size to the patient's The standard skeletal muscle model;

(4b) Use the ICP algorithm to obtain the transformation matrix between the patient's bone model after step (2e) is reset and the standard skeletal muscle model obtained through step (4a), and then replace the standard skeletal muscle model with the reset one Patient bone model, and set the fixed connection relationship between the muscle attachment point and the distal bone model and the proximal bone model. If the muscle attachment point is close to the distal bone model, it is stipulated that the muscle attachment point is fixedly connected with the distal bone model; On the contrary, it is stipulated that the muscle attachment point is fixedly connected with the proximal bone model;

(4c) First, use the C++ API interface of OpenSim software to establish a six-degree-of-freedom joint between the distal bone model and the proximal bone model in the three-dimensional virtual scene. The six-degree-of-freedom joint has t y t z t coordinate axis x t , y t , z t three translation degrees of freedom and three rotation degrees of freedom around the coordinate axis x t , y t , z t ; secondly, establish the distal bone model along the coordinate axis The translational distance d x , d y , d z of x t , y t , z t and the angle α, β, γ of the rotation of the distal bone model around the coordinate axes x t , y t , z t and the position of the distal bone model The relational expression of pose T i ; finally, take T 1 as the pose of the distal bone model in the initial state, and T t as the target pose, to simulate the movement of the distal bone model relative to the proximal bone model; if the distal bone model is known T i of the bone model, the joint parameters d x , d y , d z , α, β, γ can be obtained; the joint parameters d x , d y , d z , α, β, γ are used as the input of the OpenSim software , then the OpenSim software outputs the pulling force of each muscle corresponding to the pose T i ;

(5) On the basis of the original execution process of the A* algorithm, the collision detection method in step (3) and the muscle force analysis method in step (4) are combined, and the trajectory nodes and evaluation functions of the A* algorithm are redesigned, Perform fracture reduction trajectory planning, the specific steps are:

(5a) Define the trajectory node N i (O xi , O yi , O zi ,θ i ) (i=1,2,…,t) corresponding to the pose T i of the distal bone model one-to-one, and the trajectory node N i The components O xi , O yi , O zi are the x, y, and z components of the origin O i of the coordinate system O i -xi y i z i respectively , and the component θ i represents the axis of the distal bone model around the distal bone model n the angle of counterclockwise rotation;

(5b) Combining the constraints of no collision, the shortest path, and the smallest muscle pulling force on the fracture reduction trajectory, the evaluation function of the A* algorithm is defined as:

Among them, g(N i ) is the cost of moving the origin of the distal bone model coordinate system from the initial node N 1 to node N i ; h(N i ) indicates that the origin of the distal bone model coordinate system moves from node N i to the target node N The cost of t ; r(N i )=c 1 |θ i -θ t | is the attitude penalty function, which represents the cost required for the distal bone model coordinate system to rotate from the current node N i to the target node N t , and c 1 is Correlation penalty coefficient, the value of c 1 should make r(N i ) have the same order of magnitude as g(N i ), h(N i ), m(N i ); m(N i )=c 2 (F i -F 1 ) is the muscle force penalty function, F i (i=1,2,...,t) represents the algebraic sum of the muscle forces when the distal bone model is located at the trajectory node N i , wherein, F i (i=1, 2,...,t) is the algebraic sum of each muscle force obtained by using the method in step (5a) to map the node N i to the pose T i of the distal bone model, and then using the method in step (4c); c 2 is the correlation penalty coefficient, the value of c 2 should make m(N i ) have the same order of magnitude as g(N i ), h(N i ), r(N i );

(5c) Based on the improved A* search algorithm, the trajectory nodes N 1 (0,0,0,0) and N t (O xt ,O xt ,O xt ,θ t ) set in step (5a) are used as the trajectory respectively planning start node and end node, and optimize the reset trajectory node according to the evaluation function f(N i ) set in step (5b), in the optimization process, use the method in step (3) to select all neighbors of N c The domain nodes perform collision detection, and the reset trajectory node sequence N 2 ,…,N t-1 is optimized, and then the trajectory node N i (O xi , O yi , O zi ,θ i ) and the distal bone pose T i According to the corresponding relationship between them, the optimized reset trajectory node sequence N 2 ,...,N t-1 is converted into the pose of the distal bone model to complete the fracture reduction trajectory planning.
According to claim 1, the automatic reset trajectory planning method for parallel fracture surgical robots is characterized in that: the specific process of optimizing fracture reset trajectory nodes is as follows:

Step1: Define empty tables open and close, and insert the starting node N 1 of the reset trajectory into the open table;

Step2: Determine whether the open table is empty, if open is empty, end the search; if it is not empty, find the node N c (O xc ,O xc ,O xc ,θ c ) with the smallest value of the evaluation function in the open table, record is the current node, then delete the node from the open table and insert it into the close table;

Step3: The neighbor node of the current node N c in the close table is N c+1 , and then use the method in step (5a) to map N c+1 to the pose of the distal bone model, and use the method in step (3) The method performs collision detection on all neighborhood nodes of N c , and defines the collection of non-collided nodes as a child table;

The relationship between the current node Nc and the neighbor node Nc +1 of the current node Nc is as follows:

Step4: Traverse each node N j in the child table, and calculate the new evaluation function value of node N j when N c is the parent node; if the open table contains node N j , and the new evaluation function value of node N j is smaller than the old value value of the evaluation function, then set the parent node of N j to N c ; if the open table does not contain node N j , set the parent node of N j to N c , and add N j to the open table;

Step5: Determine whether the end node N t is included in the open table, and if so, terminate the search process and execute Step6, otherwise return to Step2 and re-execute;

Step6: When the target node N t is searched, it means that the trajectory planning is completed. Starting from the end node N t , the parent node is continuously searched upward until the parent node is the starting node N 1. N 1 , N 2 , ... N t together constitute the trajectory nodes of fracture reduction, and the above nodes are converted into the pose of the distal bone model to complete the trajectory planning of fracture reduction.