Dental implantation robot path generation method

Abstract

The invention provides a dental implant robot path generation method. The method specifically comprises the following steps: constructing a tooth space pose vector; constructing a regression model for predicting spatial features of target missing teeth by combining spatial pose features of opposite-side homonymous teeth and upper and lower occlusal teeth, and predicting spatial positions and attitude angles of the target missing teeth; tCP calibration is carried out on the dental implant robot; carrying out TCF calibration on the dental implant robot; constructing a spatial transformation relation between an end effector of a mechanical arm of the robot and an external visual positioning system; constructing a spatial transformation matrix between the digital model and the entity printing model by adopting a point cloud registration technology; an initial obstacle avoidance path of the mechanical arm is constructed based on a fast expansion random tree sampling algorithm, and meanwhile, smoothing and speed regulation are performed on a joint track through polynomial interpolation and B spline curve interpolation in combination with an improved algorithm RRT *. According to the technical scheme, the problems that in the prior art, planting errors are large, and the fitting degree is not high are solved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of path generation, and particularly relates to a tooth implant robot path generation method. BACKGROUND

[0002] Since the implant tooth body is similar in function to the original tooth body, the effect after treatment is obvious and there is no adverse reaction. For the current clinical diseases such as dentition defect and tooth loss, tooth implant repair is the best treatment scheme. However, the precision and treatment effect of the free-hand implantation mode are closely related to the experience and operation environment of the doctor. For example, visual obstruction may cause implantation errors or differences in treatment effects due to specific physiological structures of patients.

[0003] Therefore, there is a need for a tooth implant robot path generation method with small implantation errors and high fit. SUMMARY

[0004] The main purpose of the present application is to provide a tooth implant robot path generation method to solve the problem of large implantation errors and low fit in the prior art.

[0005] To achieve the above-mentioned purpose, the present application provides a tooth implant robot path generation method, which specifically comprises the following steps:

[0006] S1, extracting an axis in a tooth space based on a CBCT image, expressing the spatial relationship between the tooth axis and a standard Cartesian coordinate system, and constructing a tooth space pose vector.

[0007] S2, combining the spatial pose features of the contralateral homonym tooth and the upper and lower occlusal teeth to construct a regression model for predicting the spatial features of the target missing tooth, and predicting the spatial position and attitude angle of the target missing tooth.

[0008] S3, TCP calibration of the tooth implant robot.

[0009] S4, TCF calibration of the tooth implant robot.

[0010] S5, constructing the spatial transformation relationship between the robot manipulator end effector and the external visual positioning system.

[0011] S6, using point cloud registration technology to construct a spatial transformation matrix between the digital model and the physical printing model.

[0012] S7, constructing an initial obstacle avoidance path of the manipulator based on a rapidly-exploring random tree sampling algorithm, combining the improved algorithms RRT-Connect and RRT*, and smoothing and speed regulating the joint trajectory through polynomial interpolation and B-spline curve interpolation.

[0013] Further, step S1 specifically comprises the following steps:

[0014] S1.1, import the CBCT scanning data into a medical image processing platform, perform three-dimensional reconstruction, and obtain a dental voxel model.

[0015] S1.2, select three layers in the dental voxel model, respectively corresponding to the upper segment, middle segment and apical region of the dental root, and mark them as the first, second and third layers. In each layer, obtain a set of dental outer contour boundary points wherein k = 1, 2, 3 represents the layer number, n k represents the number of boundary points extracted in the layer, and the two-dimensional geometric center C (k) of each layer is calculated as follows:

[0016]

[0017] The two-dimensional geometric center of each layer is upgraded to a three-dimensional coordinate according to the corresponding layer height z k to obtain three groups of space points P1, P2 and P3:

[0018] P1 = (x1, y1, z1), P2 = (x2, y2, z2), P3 = (x3, y3, z3).

[0019] S1.3, using the least squares method, fit P1, P2 and P3 into the dental center axis, and the direction vector of the dental center axis is represented as:

[0020]

[0021] wherein (a, b, c) is a unit vector, representing the components of the dental center axis in the X, Y and Z directions.

[0022] S1.4, the spatial relationship between the dental center axis and the standard Cartesian coordinate system is represented as:

[0023]

[0024] wherein α, β and γ represent the direction cosine angles of the X, Y and Z axes, respectively.

[0025] S1.5, the dental space pose vector is represented as:

[0026] ξ = [x, y, z, α, β, γ] T ;

[0027] wherein (x, y, z) is the three-dimensional space coordinate of the tooth.

[0028] Further, step S2 specifically includes the following steps:

[0029] S2.1, assuming that the space pose of the target missing tooth is ξ tar :

[0030] ξ tar = [x tar , y tar , z tar , a tar , b tar , g tar ] T .

[0031] S2.2, define the spatial pose vector of the contralateral homonymous tooth contra ∈ R 6 ; the spatial pose vector of the occlusal tooth occlu ∈ R 6 . S2.3, construct the regression model of the predicted target missing tooth spatial feature tar :

[0032] ξ tar = W1·ξ contra +W2·ξ occlu +b+ε;

[0033] wherein W1, W2∈ R 6×6 , respectively represent the coefficient matrix of the contralateral tooth and the occlusal tooth predicting the contribution to the target tooth; b∈ R 6 is the intercept term; ε∈ R 6 is the residual term.

[0034] The regression model is expanded as follows:

[0035]

[0036] wherein tar,m ξ represents the mth feature component of the target missing tooth.

[0037] Further, step S3 specifically comprises the following steps:

[0038] S3.1, select a fixed spatial reference point, and the mechanical arm of the dental implant robot clamps the tool through the clamping mechanism, repeatedly positions the tool tip at the spatial reference point.

[0039] S3.2, control the end flange of the mechanical arm to contact the spatial reference point in at least four significantly different poses.

[0040] S3.3, record the corresponding mechanical arm end flange coordinates.

[0041] S3.4, construct the calibration equation set and use the least squares method to fit, obtain the tool center point TCP position

[0042]

[0043] wherein, and respectively represent the rotation matrix of the end flange coordinate system relative to the dental implant robot base coordinate system in four different postures of the mechanical arm.

[0044] Further, the step S4 specifically comprises the following steps:

[0045] S4.1, marking the end flange center point in the initial pose of the mechanical arm movement as o point; dragging the mechanical arm from the o point along the X axis by a distance Δx; dragging the mechanical arm along the z axis from the o point by a distance Δz; establishing the following coordinate relationship:

[0046]

[0047] wherein, is the coordinate of the end flange center point in the base coordinate system after moving the distance Δx; is the coordinate of the end flange center point in the base coordinate system in the initial pose; is the coordinate of the tool end point in the base coordinate system after moving the distance Δx; is the coordinate of the tool end point in the base coordinate system in the initial pose; is the coordinate of the end flange center point in the base coordinate system after moving the distance Δz; is the coordinate of the tool end point in the base coordinate system after moving the distance Δz.

[0048] S4.2, constructing the rotation matrix of the tool coordinate system according to the three calibration points, and obtaining the orthogonal basis through vector difference and cross multiplication operation:

[0049]

[0050] o = a x n;

[0051] wherein, is the inverse matrix of the rotation matrix of the end flange center point in the base coordinate system in the initial pose, is the rotation matrix of the end flange center point in the base coordinate system, is the coordinate of the tool end point in the end flange coordinate system in the initial pose; n is the rotation matrix relative to the x axis, a is the rotation matrix relative to the z axis, and o is the rotation matrix relative to the y axis.

[0052] Further, the step S5 specifically comprises the following steps:

[0053] S5.1, since the marker C is fixed to the end flange E, then at time t and time t+1:

[0054]

[0055] where, is the rotation matrix of the end flange coordinate system relative to the end position marker coordinate system.

[0056] S5.2, the transformation matrix between the end flange coordinate system and the end position marker coordinate system is given by the homogeneous transformation relationship:

[0057]

[0058] where, is the transformation matrix of the robot base coordinate system relative to the end flange coordinate system, is the transformation matrix of the optical position tracker coordinate system relative to the robot base coordinate system, is the inverse of the transformation matrix of the optical position tracker coordinate system relative to the end position marker.

[0059] Further, it is obtained that:

[0060]

[0061] After equivalent conversion, the basic expression of hand-eye calibration is obtained as AX = XB:

[0062]

[0063] where, A and B are the transformation matrices of the mechanical arm system and the vision system between adjacent two poses, respectively, and X is the fixed transformation matrix to be solved; R A , R X and R B are the rotation sub-matrices of the matrix, and b A , b X and b B are the translation sub-matrices of the matrix, respectively.

[0064] S5.3, after the homogeneous matrix is decomposed into rotation and translation parts, it is further decomposed as:

[0065]

[0066] where, t X is the translation vector.

[0067] S5.4, the Tasi-Lenz method is used to model the rotation matrix by axis-angle representation, and R A R X = R X R B , it is obtained that Taking the logarithm of , it is obtained that:

[0068]

[0069] Using the similarity transformation property of matrix logarithms:

[0070]

[0071] Since a three-dimensional rotation matrix is ​​obtained by performing matrix exponentiation operations on the antisymmetric matrix form of the rotation vector:

[0072]

[0073] Among them, [α] × and [β] × It is an antisymmetric matrix generated by rotation vectors α and β.

[0074] S5.5, this problem is a classic absolute orientation problem, solved using Singular Value Decomposition (SVD) to find R. X :

[0075] S5.6, After solving the rotation part, the translation vector t is solved using the following formula. X :

[0076] (R Ai -I)t X =R X t Bi -t Ai .

[0077] S5.7 introduces the LM algorithm for nonlinear optimization. The iterative increment of the LM algorithm satisfies:

[0078] (J T J+λI)Δp=-J T r;

[0079] in, ρ is the parameter vector, r is the rotation part of the rotation vector parameterization, and r is the residual vector composed of rotation and translation residuals; J is the Jacobian matrix of residuals with respect to parameters; I is the identity matrix; Δp is the parameter update amount.

[0080] S5.8, the objective function E(p) is defined as the sum of squared rotation and translation errors:

[0081]

[0082] Where R(ρ)=exp([ρ] × ) represents the rotation matrix corresponding to the rotation vector.

[0083] vee(·) converts an antisymmetric matrix into a three-dimensional vector; They are The corresponding rotation matrix; are corresponding translation vector.

[0084] S5.9, solve the equation (J T J+λI)Δp=-J T r, update the parameters p until the objective function E(p) converges to a preset threshold, to obtain the optimal estimate of the rotation matrix and the displacement vector.

[0085] Further, step S6 specifically comprises the following steps:

[0086] S6.1, downsample the original point cloud, and calculate the expression:

[0087]

[0088] wherein, is the center point in the voxel, p i is the original point in the voxel, and n is the number of voxel points;

[0089] S6.2, use the FPFH algorithm to construct a feature histogram by counting the relative angle relationship between point normal vectors.

[0090] S6.3, based on the initial transformation matrix obtained in the coarse registration, use the ICP algorithm for fine registration, and realize accurate alignment of the point cloud by iteratively minimizing the Euclidean distance between the two sets of corresponding points, and the objective function of ICP is:

[0091]

[0092] wherein, p i is the i-th point of the source point set, q i is the corresponding point of the target point set, R is the rotation matrix, t is the translation vector, and min is the minimum.

[0093] Further, step S7 specifically comprises the following steps:

[0094] S7.1, use the axis-aligned bounding box (AABB) to structurally model each object, and the AABB wraps the object Ω in a minimum cuboid parallel to the coordinate axes, which is defined as:

[0095] AABB(Ω)=[x min ,x max ]×[y min ,y max ]×[z min ,z max ];

[0096] wherein, x min and x max represent the minimum and maximum boundaries of the object Ω on the x-axis, ymin ,y max Let z represent the minimum and maximum boundaries of object Ω on the y-axis, respectively. min ,z max These represent the minimum and maximum boundaries of object Ω on the z-axis, respectively.

[0097] Collision determination is based on whether the AABB coordinates of two objects A and B overlap in all three axes.

[0098]

[0099] S7.2, Path search is based on collision detection results, and the RRT algorithm is first used; RRT randomly samples points q in the continuous configuration space C. rand And find the node q that is closest to the random sampling point in the existing path tree. near :

[0100]

[0101] Where, ‖·‖ represents the Euclidean distance, Q is the set of all nodes in the current tree, and arg min is the position of the minimum value.

[0102] Subsequently, a new node q is generated using the step-size constraint function Steer. new :

[0103]

[0104] Where η is the maximum step size.

[0105] S7.3 introduces the RRT-Connect algorithm, simultaneously starting from the starting point q. start and the endpoint q goal Each grows a random tree and The two trees are expanded alternately and attempts to connect; the connection condition is the number of nodes in the two trees. and Distance not exceeding threshold ∈ conn :

[0106]

[0107] S7.4 introduces the RRT* algorithm, which is used to expand new nodes q. new Then, the path is optimized through a local rewiring mechanism.

[0108] S7.5 uses fifth-order polynomial interpolation and B-spline curve interpolation to smooth the path points and transform them into continuous, dynamically controllable joint angle trajectories.

[0109] Furthermore, step S7.4 specifically includes the following steps:

[0110] S7.4.1, search the neighboring node set Q in the radius r n near :

[0111] Q near = {q e Q | ||q-q new || <= r n};

[0112] Wherein, r n is dynamically adjusted with the sampling number.

[0113] S7.4.2, select the neighboring node q min e Q near with the minimum total cost as the parent node of q new :

[0114]

[0115] Wherein, c(q start ,q) represents the path cost from the starting point to the node q, and Cost(q,q new ) represents the cost from q to q new , usually the path length.

[0116] S7.4.3, perform the rewiring operation on the neighboring node, and if a lower path cost can be obtained through q new , update the parent node and the path cost.

[0117] The present application has the following beneficial effects:

[0118] (1) The personalized initial path is constructed by combining statistical regression and geometric fitting, so that the systematic error caused by unreasonable direction setting is reduced.

[0119] (2) The sampling algorithm such as RRT-Connect / RRT* is adopted to complete the whole obstacle avoidance in the Cartesian space, so as to ensure the safety of the path in the anatomical sense.

[0120] (3) The seven polynomial interpolation is applied in the joint space to perform trajectory smoothing optimization, so as to suppress the end shaking caused by the sudden change of speed and acceleration, and convert the "feasible path" into "highly executable path".

[0121] Therefore, the present application realizes the comprehensive optimization of the prior art in the aspects of path rationality, execution stability and result consistency. BRIEF DESCRIPTION OF DRAWINGS

[0122] ​In order to more clearly illustrate the technical solutions in the specific embodiments of the present application or the prior art, the following will briefly introduce the drawings needed to be used in the specific embodiments or the prior art description. Obviously, the drawings described below are some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor based on these drawings. In the drawings:

[0123] Figure 1 A flow chart of a dental implant robot path generation method of the present application is shown.

[0124] Figure 2 A space mapping relationship diagram is shown.

[0125] Figure 3 A tooth distribution diagram of the present embodiment is shown.

[0126] Figure 4 A diagram of the change of the rotation angle of each axis of the upper teeth is shown.

[0127] Figure 5 A diagram of the change of the rotation angle of each axis of the lower teeth is shown.

[0128] Figure 6 A diagram of the data of the rotation angle of each axis after removing the mean value is shown.

[0129] Figure 7 A two-dimensional obstacle avoidance and path planning diagram of the RRT algorithm is shown.

[0130] Figure 8 A three-dimensional obstacle avoidance and path planning diagram of the RRT algorithm is shown.

[0131] Figure 9 An initial obstacle avoidance path diagram is shown. DETAILED DESCRIPTION

[0132] The technical solutions of the present application will be described in detail below with reference to the drawings. Obviously, the described embodiments are some embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0133] As shown in a dental implant robot path generation method, the method specifically comprises the following steps: Figure 1

[0134] S1, extracting the axis in the tooth space based on the CBCT image, expressing the spatial relationship between the tooth axis and the standard Cartesian coordinate system, and constructing the tooth space pose vector;

[0135] ​S2, combine the spatial pose features of the contralateral homonymic tooth and the upper and lower occlusal teeth to construct a regression model for predicting the spatial features of the target missing tooth, and predict the spatial position and attitude angle of the target missing tooth;

[0136] S3, TCP calibration of the dental implant robot;

[0137] S4, TCF calibration of the dental implant robot;

[0138] S5, construct the spatial transformation relationship between the robot manipulator end effector and the external vision positioning system;

[0139] S6, use point cloud registration technology to construct the spatial transformation matrix between the digital model and the physical printing model;

[0140] S7, based on the rapid extended random tree sampling algorithm, construct the initial obstacle avoidance path of the manipulator, and combine the improved algorithms RRT-Connect and RRT* to smooth and speed control the joint trajectory through polynomial interpolation and B-spline curve interpolation, aiming to realize the safe, high-precision and continuous motion control of the dental implant robot manipulator in complex space.

[0141] Specifically, step S1 specifically includes the following steps:

[0142] S1.1, import the CBCT scanning data into the medical image processing platform, perform three-dimensional reconstruction, and obtain the tooth voxel model; import the CBCT scanning data in DICOM format into the medical image processing platform such as Mimics and 3D Slicer. According to the gray value (HU) difference between tooth and bone tissue, set multiple threshold intervals to realize segmentation. Due to the sharp change in density of tooth root apex and crown-root junction area, a single threshold cannot effectively distinguish the boundary, so a multi-layer cross-sectional segmentation strategy is adopted to identify the boundary and extract the point set layer by layer for each tooth. S1.2, select three layers in the tooth voxel model, corresponding to the upper segment (near the crown end), the middle segment and the root apex region, respectively, denoted as the 1st, 2nd and 3rd layers. In each layer, the tooth outer contour boundary point set is obtained where k=1, 2, 3 represents the layer number, n k represents the number of boundary points extracted in the layer, and the two-dimensional geometric center C (k) (i.e. the centroid) of each layer is calculated as follows:

[0143]

[0144] The two-dimensional geometric center of each layer is upgraded to a three-dimensional coordinate according to the corresponding layer height z k to obtain three groups of spatial points P1, P2 and P3:

[0145] P1 = (x1, y1, z1), P2 = (x2, y2, z2), P3 = (x3, y3, z3);

[0146] S1.3, using the least squares method, P1, P2 and P3 are fitted as the dental axis, and the dental axis direction vector is

[0147]

[0148] Where (a, b, c) is a unit vector, representing the components of the dental axis in the X, Y, Z three directions;

[0149] S1.4, the spatial relationship of the dental axis and the standard Cartesian coordinate system is represented as:

[0150]

[0151] Where α, β, γ represent the direction cosine angles of X axis, Y axis and Z axis respectively;

[0152] S1.5, the above three angle parameters (α, β, γ) together constitute a complete set of pose vectors, which are used to describe the orientation characteristics of the dental body in the Cartesian space. Combined with the position coordinates (x, y, z) of the tooth in the CBCT image, the spatial pose vector is constructed, and the dental spatial pose vector is represented as:

[0153] ξ = [x, y, z, α, β, γ] T ;

[0154] Where (x, y, z) is the three-dimensional space coordinates of the tooth.

[0155] The six-degree-of-freedom pose parameters not only constitute the basis input for subsequent implant path generation and guide plate direction control, but also are the key indicators for setting the target pose in robot trajectory planning. In order to facilitate regression modeling and parameter training, the present application standardizes all dental pose parameters, unifies the numerical scale, and effectively eliminates the interference caused by individual differences. In addition, through statistical analysis of different dental Euler angle parameters, the regularity of the spatial distribution of the jawbone anatomy is revealed, providing a solid data support for the implant path adaptive generation method based on anatomical driving. In view of the high correlation of human oral structure in left-right symmetry and upper-lower jaw occlusion coupling, the present application proposes a multiple linear regression model combining the spatial pose characteristics of the contralateral homonymy tooth and the upper and lower occlusion teeth, to predict the spatial position and attitude angle of the missing tooth, so as to realize the preliminary deduction and automatic generation of individualized implant path.

[0156] Specifically, step S2 specifically comprises the following steps:

[0157] ​Variable definition and eigenvector construction:

[0158] S2.1, set the spatial pose of the target missing tooth as ξ tar :

[0159] ξ tar = [x tar , y tar , z tar , α tar , β tar , γ tar ] T ;

[0160] S2.2, define the spatial pose vector of the contralateral homonym tooth as ξ contra ∈R 6 ; the spatial pose vector of the upper and lower occlusal tooth as ξ occlu ∈R 6 ;

[0161] Regression model construction:

[0162] S2.3, construct the regression model for predicting the spatial feature ξ tar of the target missing tooth:

[0163] ξ tar = W1·ξ contra +W2·ξ occlu +b+ε;

[0164] Wherein, W1,W2∈R 6×6 , respectively representing the coefficient matrix of the contralateral tooth and the occlusal tooth for predicting the contribution to the target tooth; b∈R 6 is the intercept term; ε∈R 6 is the residual term; it is assumed to obey zero-mean normal distribution.

[0165] The regression model is expanded as follows (taking the mth feature as an example):

[0166]

[0167] Wherein, ξ tar,m represents the mth feature component of the target missing tooth.

[0168] In order to enhance the robustness and generalization ability of the regression model, the present application adopts the 5-fold cross-validation (5-Fold Cross Validation) strategy to train and evaluate the performance of the model. The coefficient of determination R 2With Mean Squared Error (MSE) as the main evaluation index, the influence of different model structures and feature combinations on the prediction performance is compared, and the regression model with the best fitting accuracy and explanation ability is finally determined. The method can effectively predict the six-degree-of-freedom spatial pose of any missing tooth position, provide quantitative parameter support for implant initial path fitting and robot trajectory planning, and realize intelligent modeling of personalized implant path based on geometric correlation. In order to uniformly describe the position and attitude relationship of each component of the system in three-dimensional space, multiple Cartesian coordinate systems are introduced for geometric modeling: robot base coordinate system {B}, patient model coordinate system {Q}, model marker coordinate system {D}, optical tracking system coordinate system {M}, end marker coordinate system {C}, end flange coordinate system {E}, drill tool coordinate system {T}. Figure 2 The mapping structure and calibration process between each coordinate system are shown.

[0169] In the coordinate architecture of the application, first, the spatial transformation relationship between each joint of the robot and the end flange is established by using the forward kinematics equation of the robot, and the basic motion model is established; then, the tool center point TCP and tool coordinate system TCF calibration method is used to accurately obtain the position and attitude of the tool carried by the end effector in the flange coordinate system; then, the rigid transformation matrix between the optical tracking system and the robot coordinate system is analyzed by the classical hand-eye calibration algorithm, and the unified mapping of the two spatial systems is realized; finally, combined with the three-dimensional point cloud registration technology, the spatial alignment between the model marker and the patient's anatomical structure is realized, and the complete mapping link of the virtual path to the physical execution space is constructed, providing high-precision spatial basis for implant path interpolation, trajectory planning and obstacle avoidance control.

[0170] In order to ensure that the robot has sufficient spatial accuracy and direction control ability when executing the implant path, the actual spatial position and attitude of the surgical tool (such as the implant handset) installed on the end effector must be calibrated. Since the clamping assembly may introduce attitude deviation, the TCP and TCF calibration method is used for modeling and error compensation of the end tool, so as to accurately obtain the pose parameters of the tool in the end flange coordinate system.

[0171] The application adopts four-point method to geometrically model the tool center point (TCP). This method requires that the tool tip is fixed at a certain stationary point in space, and the robot end is controlled to contact the point at multiple different attitudes, and the spatial coordinates of the end flange at each attitude are recorded. According to the geometric constraint equation constructed by these discrete poses, the least squares method is used to solve the three-dimensional displacement vector of TCP in the flange coordinate system, so as to obtain the accurate spatial position.

[0172] Specifically, step S3 specifically includes the following steps:

[0173] S3.1, select a fixed spatial reference point (such as a tip reference), the robotic dental implant machine arm holds the tool (such as a dental drill) through a clamping mechanism, repeatedly positions the tool tip at the spatial reference point;

[0174] S3.2, control the mechanical arm flange to contact the spatial reference point with at least four significantly different poses;

[0175] S3.3, record the corresponding mechanical arm flange coordinates;

[0176] S3.4, construct a calibration equation set and use least squares fitting to obtain the tool center point TCP position

[0177] During the calibration process, the tool is rigidly installed on the mechanical arm flange face through the clamping mechanism, and the relative pose between the tool and the flange remains unchanged. Therefore, the tool system can be regarded as the "seventh joint" of the mechanical arm for unified modeling. To describe the spatial transformation relationship between the tool and the end flange, the invention introduces the representation method of homogeneous transformation matrix, and its transformation relationship can be formally expressed as:

[0178]

[0179] The conversion relationship is expanded by block matrix, and the linear relationship between the position vectors can be obtained:

[0180]

[0181] Since all calibration points are the same point in space (i.e. the calibration tip point remains constant), we have:

[0182]

[0183] Solving the above two equations together, and differentiating the data under different poses, the following linear equation set can be constructed:

[0184]

[0185] Under ideal conditions, the augmented matrix and the coefficient matrix should both have a rank of 3, but due to measurement errors and mechanical arm repositioning errors, it may lead to incompatible solutions. Therefore, the least squares method is used to solve the generalized inverse, and then the optimal approximation solution is obtained:

[0186]

[0187] Since the coefficient matrix is a non-square matrix and the columns are full rank, the invention selects the normal equation method for calculation, which is in the form of: + = (A T A) -1 A T ;

[0188] Substituting A into the above equation, we have:

[0189]

[0190] wherein, and respectively represent the rotation matrix of the end flange coordinate system relative to the dental implant robot base coordinate system in four different poses of the mechanical arm;

[0191] To improve the calibration accuracy, it is necessary to reasonably screen the collected end pose data. Specifically, data sets with small pose change amplitude should be removed to reduce the influence of TCP offset error on the calibration result. The pose information of the mechanical arm end is expressed by Euler angles (RPY) and further converted into a rotation matrix form for pose solving and unified modeling in the homogeneous transformation matrix.

[0192]

[0193] is the rotation matrix of the end flange coordinate system relative to the robot base coordinate system.

[0194] After substituting the obtained multiple sets of end pose data into the above geometric model, the three-dimensional coordinates of the tool center point (TCP) in the end flange coordinate system can be solved, thereby completing the consistency calibration between the end tool and the mechanical arm coordinate system. This calibration process not only provides accurate initial parameters for subsequent spatial trajectory interpolation and pose control of the implant path, but also significantly reduces the path deviation caused by tool installation error, ensuring the spatial consistency and precision consistency between the robot execution path and the planned path.

[0195] After completing the positioning of the end tool center point TCP, in order to further realize accurate control of the end pose, the tool coordinate system TCF needs to be calibrated. TCF calibration aims to establish a complete orientation reference system of the end tool in three-dimensional space, providing accurate coordinate basis for pose interpolation of the implant path and control of the mechanical arm end.

[0196] The present application adopts the Z-X direction three-point method (Three-Point Method along Z-X axes) for TCF calibration. This method moves the mechanical arm end along the Z axis and X axis directions of the actual orientation of the tool respectively, and records three representative space points, from which the three orthogonal axis vectors of the tool coordinate system are derived, and a right-handed coordinate frame is established.

[0197] Specifically, step S4 specifically includes the following steps:

[0198] S4.1, mark the center point of the end flange as o point in the initial pose of the robot arm movement, drag the robot arm along the X axis from o point by a distance Δx, and drag the robot arm along the z axis from o point by a distance Δz;

[0199] Since the tool is rigidly connected to the end flange, its own pose remains unchanged during the dragging process, so the following coordinate relationship is established:

[0200]

[0201] wherein, is the coordinate of the center point of the end flange in the base coordinate system after moving a distance Δx; is the coordinate of the center point of the end flange in the base coordinate system in the initial pose; is the coordinate of the tool end point in the base coordinate system after moving a distance Δx; is the coordinate of the tool end point in the base coordinate system in the initial pose; is the coordinate of the center point of the end flange in the base coordinate system after moving a distance Δz; is the coordinate of the tool end point in the base coordinate system after moving a distance Δz;

[0202] S4.2, construct the rotation matrix of the tool coordinate system according to the three calibration points, and obtain the orthogonal basis through vector difference and cross multiplication operation:

[0203]

[0204] o = a x n;

[0205] wherein, is the inverse matrix of the rotation matrix of the center point of the end flange in the base coordinate system in the initial pose, is the rotation matrix of the center point of the end flange in the base coordinate system, is the coordinate of the tool end point in the end flange coordinate system in the initial pose; n is the rotation matrix relative to the x axis, a is the rotation matrix relative to the z axis, and o is the rotation matrix relative to the y axis.

[0206] Considering the existence of experimental errors, the three vectors obtained by preliminary calculation may not fully satisfy the orthogonal condition, so the z axis vector needs to be corrected to ensure the orthogonality of the base matrix. Finally, the rotation matrix R of the tool coordinate system is constructed through normalization processing, thereby completing the pose calibration of the TCF.

[0207] The method achieves a good balance between precision and practicability, and is particularly suitable for dental implant robot tasks with high requirements for end posture control precision. Combined with the TCP calibration result, the TCF calibration effectively ensures the spatial consistency and motion stability of the robot end effector in the path planning and posture control process.

[0208] To realize the accurate mapping of the target planting path in the space of the mechanical arm, it is necessary to construct the spatial transformation relationship between the end effector of the mechanical arm and the external vision positioning system. Hand-eye calibration, as a common spatial registration technology, aims to solve the rigid transformation matrix between the robot end and the optical tracking system. According to the installation mode of the vision sensor, hand-eye calibration can be divided into "eye on hand" and "eye off hand" two modes. The present application adopts "eye off hand" configuration, that is, a spatial identifier is installed at the end of the mechanical arm, and the pose information thereof is collected in real time by the fixed optical positioning system.

[0209] Specifically, step S5 specifically includes the following steps:

[0210] S5.1, the change of the pose of the end identifier caused by the movement of the mechanical arm can be mapped between the robot coordinate system and the vision coordinate system through a fixed rigid transformation X. Because the identifier C is fixed to the end flange E, at time t and time t+1, there is:

[0211]

[0212] Among them, is the rotation matrix of the end flange coordinate system relative to the end positioning marker coordinate system;

[0213] S5.2, the transformation matrix between the end flange coordinate system and the end positioning marker coordinate system is given by the homogeneous transformation relationship:

[0214]

[0215] Among them, is the transformation matrix of the robot base coordinate system relative to the end flange coordinate system, is the transformation matrix of the optical positioning tracker coordinate system relative to the robot base coordinate system, is the inverse of the transformation matrix of the optical positioning tracker coordinate system relative to the end positioning marker;

[0216] Further, it is obtained that:

[0217]

[0218] After equivalent conversion, the basic expression of hand-eye calibration AX=XB is obtained:

[0219]

[0220] where A and B are the transformation matrices between two adjacent poses of the robot system and vision system respectively, and X is the fixed transformation matrix to be solved. A , R X and R B are the rotation sub-matrices of matrixes; A , b X and b B are the translation sub-matrices of matrixes respectively.

[0221] S5.3, after the homogeneous matrix is decomposed into rotation and translation parts, it is further decomposed into:

[0222]

[0223] where t X is the translation vector.

[0224] S5.4, to reduce the redundancy of solving, the Tasi-Lenz method is used to model the rotation matrix with axis-angle representation to avoid the parameter redundancy problem existing in the traditional rotation matrix. The three-parameter representation is used in the axis-angle form, which significantly reduces the calculation dimension and improves the robustness. The efficient conversion between the axis-angle representation and the rotation matrix can be realized through the Rodrigues formula. From R A R X = R X R B , we get Take the logarithm of , then:

[0225]

[0226] Using the similarity transformation property of matrix logarithm:

[0227]

[0228] Since a three-dimensional rotation matrix is obtained through the matrix exponential operation of the skew-symmetric matrix form of the rotation vector:

[0229]

[0230] where [α] × and [β] × are the skew-symmetric matrices generated by the rotation vectors α and β.

[0231] S5.5, solve R X is the classical absolute orientation problem, which uses singular value decomposition (SVD) to solve R X :

[0232]

[0233] S5.6, After solving the rotation part, the translation vector t is solved using the following formula. X :

[0234] (R Ai -I)t X =R X t Bi -t Ai ;

[0235] The linear equations in steps S5.7 and S5.6 can be estimated using the least squares method. However, considering the influence of noise in actual measurement data, to improve the accuracy and robustness of the solution, the LM algorithm is introduced for nonlinear optimization. The LM algorithm combines the Gauss-Newton method and the gradient descent method, and dynamically adjusts the iteration step size and direction by adjusting the damping factor λ. The iteration increment of the LM algorithm satisfies:

[0236] (J T J+λI)Δp=-J T r;

[0237] in, ρ is the parameter vector, r is the rotation part of the rotation vector parameterization, and r is the residual vector composed of rotation and translation residuals; J is the Jacobian matrix of residuals with respect to parameters; I is the identity matrix; Δp is the parameter update amount.

[0238] S5.8, the objective function E(p) is defined as the sum of squared rotation and translation errors:

[0239]

[0240] Where R(ρ)=exp([ρ] × ) represents the rotation matrix corresponding to the rotation vector.

[0241] vee(·) converts an antisymmetric matrix into a three-dimensional vector; They are The corresponding rotation matrix (an orthogonal matrix with a determinant of 1); They are The corresponding translation vector.

[0242] Jacobian matrix J i The expression form is:

[0243]

[0244] By iteratively solving the equation (J) T J+λI)Δp=-J T r, update parameter p until the objective function E(p) converges to the preset threshold, and obtain the optimal estimate of the rotation matrix and displacement vector.

[0245] To realize the accurate mapping of the virtual planting path in the real operation space, the application adopts a point cloud registration technology to construct a spatial rigid transformation matrix between the digital model and the physical printing model. The method discretizes the three-dimensional planting model into a point cloud form, combines the posture information of the physical sample obtained by the optical positioning system, realizes high-precision pose alignment of the virtual simulation space and the physical space, and ensures accurate positioning of the regression prediction path in the real environment.

[0246] In view of the large-scale error and local noise existing in the initial data of the point cloud, a "coarse-precision combined" rigid registration process is designed, including downsampling, feature extraction and coarse registration, and fine registration.

[0247] Specifically, step S6 specifically includes the following steps:

[0248] S6.1, voxel grid downsampling

[0249] In order to reduce the calculation complexity of the point cloud and maintain the integrity of the geometric structure, the original point cloud is first downsampled. The specific method is to divide the point cloud into a three-dimensional voxel grid with a fixed side length, and replace each voxel with the geometric center of all points in the voxel to form a sparse point set. The calculation expression is:

[0250]

[0251] Wherein, is the center point in the voxel, p i is the original point in the voxel, and n is the number of voxel points. This process reduces the number of points while retaining local geometric information.

[0252] S6.2, FPFH feature extraction and RANSAC coarse registration

[0253] In order to realize the preliminary coarse registration of the point cloud, the local geometric features of each point need to be described. The FPFH (Fast Point Feature Histogram) algorithm is adopted to construct a feature histogram by counting the relative angle relationship between the point normal vectors. The specific steps are:

[0254] ①Calculate the simplified point feature histogram (SPFH) of the query point and its neighborhood points, encode the angle features α, φ, θ, and represent the relative direction between the normal vectors;

[0255] ②Aggregate the SPFH using weighted neighborhood information to form the final FPFH feature:

[0256]

[0257] Wherein, ω i is the weight of the query point p to the neighborhood point p iEuclidean distance reciprocal weight between them, k is the number of neighborhood points.

[0258] Based on the FPFH descriptor, RANSAC (Random Sample Consensus) is used for coarse registration. This method constructs a candidate registration model by multiple random sampling, eliminates abnormal points, selects an inner point set, and then calculates the rigid transformation matrix. This process enhances the stability and robustness of the initial registration.

[0259] S6.3, ICP fine registration: Based on the initial transformation matrix obtained by coarse registration, ICP (Iterative Closest Point) algorithm is used for fine registration, and the accurate alignment of point cloud is realized by iterative minimization of the Euclidean distance between the two sets of corresponding points. The objective function of ICP is:

[0260]

[0261] Where p i is the i-th point of the source point set, q i is the corresponding point of the target point set, R is the rotation matrix, and t is the translation vector.

[0262] The solving steps are as follows:

[0263] ① Calculate the weighted centroid of the source point set and the target point set, set as and Let the partial derivative of the translation vector be zero, that is:

[0264]

[0265] The translation vector is derived as:

[0266]

[0267] ② Construct the covariance matrix:

[0268]

[0269] ③ Singular value decomposition (SVD) is performed on the covariance matrix H:

[0270] H = UΣV T ;

[0271] ④ The optimal rotation matrix is calculated as:

[0272] R = VU T ;

[0273] ⑤ Substitute the rotation matrix into the translation calculation formula to obtain the complete rigid transformation matrix, realizing the best rigid matching between the two point clouds.

[0274] The registration process combines voxel down-sampling, feature matching based on FPFH and RANSAC coarse registration, and SVD-based ICP fine registration, effectively solving the point cloud data error and noise influence, and realizing the sub-millimeter spatial alignment between the virtual digital model and the physical printed model, providing a reliable geometric basis for subsequent mechanical arm trajectory import and high-precision surgery execution.

[0275] In the limited space in the oral cavity, the dental implant robot path planning needs to overcome many space limitations and complex obstacles to ensure that the mechanical arm movement path not only avoids key anatomical structures, but also meets the real-time response and path executability. Therefore, a collision detection mechanism based on geometric boundaries is designed, and combined with RRT and its improved algorithm RRT-Connect / RRT*, efficient path search and generation are realized.

[0276] Specifically, step S7 specifically includes the following steps:

[0277] S7.1, environment modeling and collision detection: the surgical environment includes patient oral model, mechanical arm body, guide plate device and operating table and other obstacles.

[0278] Structured modeling of each object is carried out by using axis-aligned bounding box (AABB). AABB wraps the object Ω in the smallest cuboid parallel to the coordinate axes, which is defined as:

[0279] AABB(Ω)=[x min ,x max ]×[y min ,y max ]×[z min ,z max ];

[0280] Where x min and x max represent the minimum and maximum boundaries of object Ω on the x-axis, y min and y max represent the minimum and maximum boundaries of object Ω on the y-axis, and z min and z max represent the minimum and maximum boundaries of object Ω on the z-axis.

[0281] Collision determination is based on whether the AABB of two objects A and B has interval overlap Collision in three-axis direction:

[0282]

[0283] If the above conditions are met, it is determined that a collision may occur, and the system will start a more detailed collision detection algorithm to reduce false positives. This method has low computational complexity and is suitable for dynamic detection in real-time environments.

[0284] S7.2, Path search algorithm: RRT with RRT-Connect and RRT*:

[0285] Path search is based on the collision detection results, first using RRT algorithm; RRT randomly samples points q rand in the continuous configuration space C near :

[0286]

[0287] Where ||·|| represents the Euclidean distance, is the set of all nodes in the current tree, and arg min is the position of the minimum value.

[0288] Then, the step limit function Steer is used to generate new nodes q new :

[0289]

[0290] Where η is the maximum step size, ensuring that the path is smooth and the increment is reasonable. After generating new nodes, the collision detection function CollisionFree(q near , q new ) is called to ensure that the path segment is not blocked by obstacles. The path tree is continuously expanded, gradually approaching the target node.

[0291] S7.3, To improve search efficiency and path quality, the RRT-Connect algorithm is introduced, and a random tree is grown from the starting point q start and the ending point q goal : The two trees are alternately expanded and connected; the connection condition is that the distance between the nodes in the two trees and does not exceed the threshold ∈ conn :

[0292]

[0293] Once connected successfully, a collision-free path from the starting point to the ending point is obtained. RRT-Connect significantly reduces search time and path redundancy, and is suitable for high-degree-of-freedom, multi-constrained surgical robot path planning requirements, but the generated path is not globally optimal.

[0294] S7.4, To achieve the progressive optimality of the path, the RRT* algorithm is introduced, and after expanding the new node q new , the path is optimized through a local rewiring mechanism.

[0295] ​S7.5, the path points are smoothed and converted into continuous and dynamically controllable joint angle trajectories by using quintic polynomial interpolation and B-spline curve interpolation. The path generated by RRT* is a discrete node sequence, and there are problems of trajectory discontinuity and velocity mutation. If it is directly executed, it will cause the mechanical arm to vibrate and positioning error. Therefore, the path points need to be smoothed and converted into continuous and dynamically controllable joint angle trajectories. The present application uses the following two interpolation methods:

[0296] (1) Quintic polynomial interpolation

[0297] Quintic polynomial interpolation is widely used in robot trajectory planning due to its good continuity and calculation efficiency. Its trajectory expression is:

[0298] p(t) = a0 + a1t + a2t 2 +a3t 3 +a4t 4 +a5t 5 , t ∈ [t0, t f ];

[0299] Where, p(t) represents the position, the coefficients a i (i = 0, …, 5) are determined by satisfying the boundary conditions of position, velocity and acceleration:

[0300]

[0301] The velocity and acceleration are:

[0302]

[0303] Substitute the boundary conditions into the equation group to form a linear system:

[0304] Ma = b;

[0305] Where, a = [a0, a1, a2, a3, a4, a5] T , b = [p0, v0, a0, p f , v f , a f ] T ;

[0306] The coefficient matrix M is constructed according to t0 and t f , and the coefficients a are obtained after solving, realizing the generation of continuous and smooth trajectory, significantly reducing the impact and vibration in the movement process of the mechanical arm, meeting the demand of oral fine operation.

[0307] (2) B-spline curve interpolation

[0308] B-spline curve (Basis Spline) uses multi-segment polynomial basis function and control point to realize local adjustment and overall smoothness of the trajectory. The B-spline basis function is recursively defined as:

[0309]

[0310] where N i,p (t) is the basis function of order p, and {t i} is the node vector. The cubic B-spline curve (p = 3) trajectory is expressed as:

[0311]

[0312] where P i is the control point coordinate. The trajectory speed and acceleration are respectively:

[0313]

[0314] B-spline curve ensures high-order continuity (at least second order), effectively avoiding the boundary oscillation problem commonly seen in polynomial interpolation, and is suitable for scenarios where path points are dense and motion continuity is high. The disadvantage is that the calculation complexity is large, and the configuration requirements for nodes and weights are high.

[0315] Specifically, step S7.4 specifically includes the following steps:

[0316] S7.4.1, search for a set of adjacent nodes Q n within a radius r near :

[0317] Q near = {q e Q | ||q-q new || < r n};

[0318] where r n is dynamically adjusted with the sampling number to meet the asymptotic optimal condition.

[0319] S7.4.2, select the adjacent node q min e Q near that minimizes the total cost as the parent node of q new :

[0320]

[0321] where c(q start , q) represents the path cost from the starting point to node q, and Cost(q, q new ) represents the cost from q to q new , usually using path length;

[0322] S7.4.3, perform a rewire operation on the neighboring nodes if a lower path cost is obtained by q new If a lower path cost is obtained, the parent node and path cost are updated to achieve local path optimization.

[0323] Although RRT* cannot guarantee absolute optimality in finite time, the path quality gradually tends to global optimality as the number of samples increases, making it suitable for offline planning tasks with high trajectory quality requirements.

[0324] In summary, RRT algorithm is suitable for preliminary path exploration due to its simple structure and rapidity, RRT-Connect accelerates the planning process through bidirectional search, and RRT* achieves progressive optimization of the path through path cost function and rewiring mechanism. According to specific task requirements, dental implant robot path planning can flexibly choose the above algorithms to balance efficiency and path quality.

[0325] To achieve accurate initial estimation and spatial orientation reference of the implantation path, the present application extracts the central axis of the patient's full dentition and performs three-dimensional geometric fitting based on CBCT voxel data. During the entire process, first, the multi-layer cross-sectional images of each tooth in the CBCT image are extracted one by one, and the mesial, distal, labial and lingual four boundary points of the tooth are marked in each cross section. These boundary points constitute the area of the tooth cross section, and then the coordinates of these points are used to calculate the position of the center point of each cross section. Through this method, we obtain the center point sequence of each tooth in three-dimensional space, and further use these center points for fitting the central axis in subsequent analysis, see Table 1 and Table 2.

[0326] Table 1 Upper tooth voxel information table

[0327]

[0328] Table 2 Lower tooth voxel information table

[0329]

[0330] After extracting the cross-sectional center points of all teeth, the least squares fitting method (Least Squares Fitting) is used to model a three-dimensional straight line for the center point sequence of each tooth, and the central axis equation of each tooth is obtained. This method ensures the minimum error of the fitted central axis and can accurately represent the geometric shape and spatial position of the tooth in space. The fitting results of all teeth are finally summarized in the form of a set of spatial vectors (see Table 3), which represent the spatial position and attitude information of each tooth, providing a reference basis for the subsequent implantation path direction.

[0331] Table 3 Tooth fitted axis table

[0332]

[0333] A three-dimensional spatial distribution map of the tooth midline was generated using the Mimics and 3-matic 3D visualization platforms, such as... Figure 3 As shown, this illustrates the tilt direction and spatial distribution of the patient's upper and lower jaw dentition. Figure 3 The location of the missing areas is clearly revealed, especially the specific positions of LR6 and LR7 in the lower right region, which is of great reference value for subsequent path planning and implant placement selection. Visualizing the three-dimensional distribution of the patient's teeth allows for a direct observation of the continuity of the dentition and its spatial distribution patterns. Specifically, the figure shows the tilt direction of the patient's dentition, and through intuitive spatial distribution analysis, it helps confirm whether the arrangement relationship between the teeth conforms to physiological expectations.

[0334] Figure 3 The analysis also revealed the occlusal relationship between the patient's upper and lower dentitions, clearly reflecting the spatial location of the missing teeth. This analysis provides essential data support for subsequent path planning, including direction setting, implant depth determination, and obstacle avoidance space modeling. For example, by displaying the missing LR6 and LR7 areas, it can be verified whether these areas provide sufficient distance for obstacle avoidance, ensuring the feasibility of the implant path planning. In three-dimensional space, accurate positioning of the missing teeth helps ensure that the implant placement is completely consistent with the predetermined path, avoiding any potential errors or path conflicts.

[0335] To ensure the accuracy of implant pathway planning and its spatial adaptability during the planning process, a detailed statistical analysis of the spatial pose of the patient's teeth was conducted. While 3D display data alone can show the spatial distribution of teeth and the location of missing teeth, this method fails to deeply analyze the specific spatial characteristics of the teeth, potentially leading to flaws in the planning results. Therefore, a more comprehensive analysis of spatial pose and position data was employed to reveal the distribution characteristics and interrelationships of teeth in three-dimensional space.

[0336] When analyzing the position and pose of the upper and lower teeth, the changes in the rotation angles of each axis were first observed. Figure 4 and Figure 5 The figure shows the rotation angle changes of the upper and lower teeth along the X, Y, and Z axes. As can be seen from the figure, there is a strong negative correlation between the rotation angle changes along the X-axis, meaning that the rotation trends of the upper and lower teeth on this axis are opposite. Specifically, the rotation angle of the lower teeth is smaller, a phenomenon that can be attributed to the occlusal relationship between the upper and lower teeth. During occlusion, the mandible, as the primary force-generating tissue, has teeth that are more concentrated on the inner side of the oral cavity, resulting in a smaller X-axis rotation angle for the lower teeth.

[0337] A negative correlation was also observed in the change of Y-axis rotation angle, indicating that the rotation angles of the upper and lower teeth have similar trends in the Y-axis. However, a significant positive correlation was found in the change of Z-axis rotation angle, indicating that the rotation angles of the upper and lower teeth tend to be consistent and have large numerical differences. Overall, the upper and lower teeth showed a strong correlation in the inclination angle, although in some specific tooth positions (such as L6 and R1), the rotation angle change was not significant, which may be related to the biting habit or bone density conditions.

[0338] To make these data features more intuitive, mean removal processing was performed on the rotation angle data. Figure 6 The mean-removed upper and lower tooth rotation angle changes are shown. The mean-removed data shows that the rotation angles of the upper and lower teeth on the X-axis and Y-axis are negatively correlated, while on the Z-axis they show a positive correlation, further verifying the correlation between the upper and lower tooth positions. The mean-removed upper and lower tooth distribution shows a clear positive correlation, especially in the X-axis and Y-axis direction position data. Comparing the mean-removed X and Y position distribution maps, it is observed that the contralateral tooth positions exhibit mirror symmetry characteristics, indicating that the upper and lower teeth have regularity and symmetry in spatial position, which provides an important basis for subsequent implant path regression prediction.

[0339] To accurately predict the missing tooth position, the invention selects the contralateral teeth and occlusal teeth of each tooth position as the key elements for regression analysis, and obtains the spatial position and attitude of the missing tooth through regression analysis. Based on this method, the spatial position and attitude data of the upper and lower teeth are used to perform regression analysis on all tooth positions, and finally a set of regression equations are obtained. The calculation results of these regression equations accurately predict the spatial position of each tooth position, especially for the accurate positioning of the missing tooth positions (LR6 and LR7).

[0340] Specifically, the regression equation considers multiple factors such as the symmetry relationship between the upper and lower tooth positions, the position of the occlusal surface, and the distribution pattern of the surrounding teeth. Each regression equation is based on the spatial relationship and statistical data of adjacent teeth to calculate the accurate position of the missing tooth.

[0341]

[0342] The results of the regression analysis are used to calculate the position and attitude of the implant tooth through five regression equations. Table 4 shows the implant path data predicted by the regression equations, which becomes an important basis for customizing the implant scheme for the missing tooth positions LR6 and LR7. The specific calculation results are as follows:

[0343] Table 4 Implant path based on regression equation

[0344]

[0345] Table 4 lists the spatial coordinates and attitude angles of LR6 and LR7, where x, y, and z represent the spatial coordinates, and a, b, and g represent the rotational attitude of the tooth. Through these calculation results, the direction and depth of the implant path can be accurately determined to ensure the precise positioning and occlusal docking of the implant.

[0346] These regression equations provide key data support for subsequent implant path planning. In particular, when performing multi-implant bridge implantation, by selecting the average value of the rotation angle in the regression equation, the attitude of the multiple implants can be ensured to be parallel, thereby ensuring the precise arrangement and functionality of the implants. In addition, all predicted implant path data will eventually be converted into image space data for subsequent surgical operation and implant simulation.

[0347] Through this series of regression analysis and data verification, the present application successfully provides reliable prediction results for the implant paths of LR6 and LR7, and provides a mathematical basis for the precise implementation of implant surgery.

[0348] To ensure that the oral implant robot safely and stably reaches the target implant site in a narrow and complex anatomical structure surgical area, the overall trajectory planning is divided into two stages: first, a collision-free path that meets the anatomical constraints is generated in Cartesian space, and second, the path is smoothed and dynamically executable in joint space. The first stage focuses on path accessibility and obstacle avoidance reliability, and the second stage emphasizes the continuity and controllability of dynamic parameters such as velocity and acceleration.

[0349] Obstacle avoidance analysis and path generation results:

[0350] During the Cartesian space path search process, a geometric space collision detection strategy is used: an axis-aligned bounding box (AABB) is constructed for each link of the robot arm, the support table, and the patient's oral tissue. During the path expansion stage, the bounding box is first used for quick screening, and only when the bounding boxes intersect is detailed collision detection triggered, thereby significantly reducing the computational load while ensuring detection accuracy.

[0351] In terms of path search algorithms, three typical sampling-based planning methods, RRT, RRT-Connect, and RRT*, were simulated in two and three dimensions for comparison, as seen in Figure 7 and Figure 8The results in Table 5 show that RRT-Connect exhibits extremely short planning time (0.002 s and 0.003 s) in both 2D and 3D environments, thanks to its bidirectional tree expansion strategy, making it suitable for online fast response in surgical procedures; RRT* takes significantly more time (11.67 s in 2D and 6.66 s in 3D), but has the shortest path length (12.67 mm in 2D and 16.62 mm in 3D), making it gradually optimal and more suitable for offline generation of high-quality reference trajectories; RRT performs between the two, with a success rate of 100%, but is not superior in terms of time consumption and path optimization. To balance efficiency and optimality, RRT-Connect is preferred in the process to quickly obtain a feasible path, and RRT* is used for offline path quality verification and optimization when time permits, to balance efficiency and optimality. Figure 9 Further illustrate the morphology of the initial obstacle avoidance path: the first collision-free path obtained by a sampling-based planner (such as RRT*) in Cartesian space is composed of a series of discrete nodes and their connections, the starting point and target point are highlighted in different colors or markers, and the intermediate nodes are connected in polyline segments, clearly reflecting the search process of the algorithm under the constraint of obstacle distribution; the main obstacles in the robot working environment and their geometric boundaries are also depicted, and the path remains outside the envelope without penetration. It can be seen that the trajectory as a whole presents a "zigzag" trend, with many turns and large curvature changes, which is a typical "feasible but not smooth" original solution generated by random sampling, laying the foundation for subsequent smoothing optimization.

[0352] Table 5 Comparison of three algorithms in 2D and 3D environments

[0353]

[0354] The original sampling data of the three planting experiments are listed in Table 6; the corresponding top / root endpoint position deviation and attitude angle deviation are listed in Table 7.

[0355] Table 6 Point positions obtained by optical positioning

[0356]

[0357] Table 7 Analysis of experimental results and planning path deviation

[0358]

[0359] The present application is based on CBCT data, and constructs a closed-loop dental implant robot system covering "image reconstruction-path planning-multi-coordinate system registration-robot execution-effect verification". Compared with the traditional freehand operation and semi-automatic guide system, the system has obvious advantages in the spatial positioning and attitude control of the cavity preparation. The experimental results show that the linear deviation of the top and root end points of the implant path is (1.03±0.34) mm and (1.06±0.23) mm respectively, and the attitude deviation is (3.47°±1.99°). Compared with the freehand operation error reported in the existing literature, the position deviation is (1.43±0.49) mm and (2.20±0.79) mm, and the angle deviation is 6.78°±3.31°, the present application improves the position accuracy by about 28%~52%, and the attitude deviation is nearly halved. Compared with the >1.5mm position deviation and >5° attitude error commonly occurring in the conventional semi-automatic guide system, the present application also shows better precision control ability.

[0360] Of course, the above description is not a limitation of the present application, and the present application is not limited to the above examples. Changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present application should also be within the scope of the present application.

Claims

1. A dental implant robot path generation method characterized by, Specifically comprising the following steps: S1, extracting the axis in the dental space based on the CBCT image, expressing the spatial relationship between the dental axis and the standard Cartesian coordinate system, and constructing a dental spatial pose vector; S2, combining the spatial pose characteristics of the contralateral homonym tooth and the upper and lower occlusal teeth to construct a regression model for predicting the spatial characteristics of the target missing tooth, and predicting the spatial position and attitude angle of the target missing tooth; S3, TCP calibration of the dental implant robot; S4, TCF calibration of the dental implant robot; S5, constructing the spatial transformation relationship between the robot mechanical arm end effector and the external visual positioning system; S6, using point cloud registration technology to construct the spatial transformation matrix between the digital model and the physical printing model; S7, constructing the initial obstacle avoidance path of the mechanical arm based on the rapidly-exploring random tree sampling algorithm, combining the improved algorithms RRT-Connect and RRT*, and smoothing and speed regulating the joint trajectory through polynomial interpolation and B-spline curve interpolation.

2. The dental implant robot path generation method of claim 1, wherein, Step S1 specifically comprises the following steps: S1.1, import the CBCT scan data into the medical image processing platform, perform three-dimensional reconstruction, and obtain the dental voxel model; S1.2, select three layers in the dental voxel model, respectively corresponding to the upper segment of the root, the middle segment and the apical region, denoted as the first, second and third layers, and on each layer, obtain the dental outer contour boundary point set wherein k = 1, 2, 3 represents the layer number, n k represents the number of boundary points extracted in the layer, and the two-dimensional geometric center C (k) of each layer is calculated as follows: The two-dimensional geometric center of each layer is calculated according to the corresponding layer height z k The three-dimensional coordinates are obtained by lifting the dimension, and three groups of spatial points P1, P2 and P3 are obtained: P1=(x1,y1,z1),P2=(x2,y2,z2),P3=(x3,y3,z3); S1.3, using least squares, fit P1, P2, and P3 to a dental axis, a dental axis direction vector is represented as: Where (a, b, c) is a unit vector representing the components of the dental axis in the X, Y, and Z directions; S1.4, the spatial relationship between the dental axis and the standard Cartesian coordinate system is expressed as: Where α, β, γ represent the direction cosine angles of the X, Y, and Z axes, respectively; S1.5, the dental spatial pose vector is expressed as: ξ = [x, y, z, a, β, γ] T ; Where (x, y, z) is the three-dimensional spatial coordinate of the tooth.

3. The dental implant robot path generation method of claim 1, wherein, Step S2 specifically comprises the following steps: S2.1, set the spatial pose of the target missing tooth as ξ tar : ξ tar = [x tar , y tar , z tar , a tar , b tar , g tar ] T ; S2.2, defining the spatial pose vector ξ of the contralateral homonymous tooth contra ∈ R 6 ; the spatial pose vector ξ of the upper and lower occlusal tooth occlu ∈ R 6 ; S2.3, construct a prediction target missing tooth spatial feature ξ tar of the regression model: ξ tar = W1 · ξ contra + W2 · ξ occlu + b + ε; wherein W1, W2 ∈ R 6×6 , respectively, denote the coefficient matrixes predicting the contribution of the opposite teeth and the occlusal teeth to the target tooth; b ∈ R 6 is the intercept term; and ε ∈ R 6 is the residual term. The regression model is expanded as follows: wherein, ξ tar,m represents the mth feature component of the target missing tooth site.

4. The dental implant robot path generation method of claim 1, wherein, Step S3 specifically comprises the following steps: S3.1, select a fixed spatial reference point, and the mechanical arm of the dental implant robot holds the tool through the clamping mechanism, repeatedly positions the tool tip at the spatial reference point; S3.2, control the mechanical arm end flange to contact the spatial reference point with at least four significantly different attitudes; S3.3, record the corresponding mechanical arm end flange coordinates; S3.4, construct calibration equation set and solve for tool center point TCP position using least squares fit wherein, and respectively represent the rotation matrix of the end flange coordinate system relative to the dental implant robot base coordinate system in four different postures of the mechanical arm.

5. The dental implant robot path generation method of claim 1, wherein, Step S4 specifically comprises the following steps: S4.1, mark the center point of the end flange at the initial pose of the mechanical arm movement as o point; drag the mechanical arm along the X axis from o point by a distance Δx; drag the mechanical arm along the z axis from o point by a distance Δz; establish the following coordinate relationship: wherein, is the coordinate of the center point of the end flange in the base coordinate system after moving a distance Δx; is the coordinate of the center point of the end flange in the base coordinate system in the initial pose; is the coordinate of the tool end point in the base coordinate system after moving a distance Δx; is the coordinate of the tool end point in the base coordinate system in the initial pose; is the coordinate of the center point of the end flange in the base coordinate system after moving a distance Δz; is the coordinate of the tool end point in the base coordinate system after moving a distance Δz; S4.2, construct the rotation matrix of the tool coordinate system according to the three calibration points, and obtain the orthogonal basis through vector difference and cross multiplication operation: o=a×n; wherein, is the inverse of the rotation matrix of the center point of the end flange in the base coordinate system under the initial pose, is the rotation matrix of the center point of the end flange in the base coordinate system, is the coordinate of the tool end point in the end flange coordinate system under the initial pose; n is the rotation matrix relative to the x-axis, a is the rotation matrix relative to the z-axis, and o is the rotation matrix relative to the y-axis.

6. The dental implant robot path generation method of claim 1, wherein, Step S5 specifically comprises the following steps: S5.1, since the marker C is fixed to the end flange E, there is: wherein, is the transformation matrix of the end flange coordinate system with respect to the end positioning marker coordinate system; S5.2, the transformation matrix between the end flange coordinate system and the end positioning marker coordinate system is given by the homogeneous transformation relationship: wherein, is a transformation matrix of the robot base coordinate system with respect to the end flange coordinate system, is a transformation matrix of the optical positioning tracker coordinate system with respect to the robot base coordinate system, is an inverse of a transformation matrix of the optical positioning tracker coordinate system with respect to the end positioning marker. Further obtained: After equivalent conversion, the basic expression of hand-eye calibration AX=XB is obtained: Wherein, A and B are transformation matrices of the robot system and the vision system between two adjacent poses, respectively, and X is the fixed transformation matrix to be solved; R A , R X , and R B are rotation sub-matrices of the matrix; b A , b X , and b B are translation sub-matrices of the matrix, respectively. S5.3, after decomposing the homogeneous matrix into rotation and translation parts, further decompose it into: where t X is the translation vector; S5.4, the rotation matrix is modeled using the axis-angle representation by R A R X = R X R B , we get Taking the logarithm, we get ​ Using the similarity transformation property of matrix logarithm: Since a three-dimensional rotation matrix is obtained by matrix exponential operation through the skew-symmetric matrix form of the rotation vector: where [a] and [β] are skew-symmetric matrices generated from the rotation vectors a and β; and × and [β] × are skew-symmetric matrices generated from the rotation vectors a and β; and S5.5, this problem is the classical absolute orientation problem, solve R using singular value decomposition (SVD) X : S5.6, after solving the rotation part, solve the translation vector t using the following formula X : (R Ai -I)t X = R X t Bi -t Ai ; S5.7, the LM algorithm is introduced for nonlinear optimization, and the LM algorithm iteration increment satisfies: (J T J+λI)Δp=J T r; (J T J+λI)Δp=-J T r; wherein, is the parameter vector, p is the rotation part of the rotation vector parameterization, r is the residual vector composed of the rotation and translation residuals; J is the Jacobian matrix of the residuals with respect to the parameters; I is the identity matrix; Ap is the parameter update; S5.8, the objective function E(p) is defined as the sum of squares of rotation and translation errors: where R(p) = exp([p] x) and p = [px, py, pz]T. × denotes the rotation matrix corresponding to the rotation vector, vee(·) is the conversion from skew-symmetric matrix to three-dimensional vector; are the rotation matrices corresponding to the rotation vectors, are the rotation matrices corresponding to the rotation vectors, are the translation vectors corresponding to the rotation vectors, are the translation vectors corresponding to the rotation vectors. S5.9, solve the equation (J T J + λI) Δp = -J T r, update the parameters p until the objective function E(p) converges to a preset threshold, obtaining the optimal estimation of the rotation matrix and the displacement vector.

7. The dental implant robot path generation method of claim 1, wherein, Step S6 specifically comprises the following steps: S6.1, the original point cloud is down-sampled, and the expression is calculated as: wherein, is the center point within the voxel, p i is the original point within the voxel, n is the number of voxel points; S6.2, the FPFH algorithm is used to construct a feature histogram by statistically analyzing the relative angle relationship between point normal vectors; S6.3, on the basis of the obtained initial transformation matrix of coarse registration, the ICP algorithm is used for fine registration, and the accurate alignment of the point cloud is realized by iteratively minimizing the Euclidean distance between the two sets of corresponding points, and the objective function of ICP is: where p i is the i-th point of the source point set, q i is the corresponding point of the target point set, R is a rotation matrix, t is a translation vector, and min is the minimum.

8. The dental implant robot path generation method of claim 1, wherein, Step S7 specifically comprises the following steps: S7.1, the axis-aligned bounding box (AABB) is used for structured modeling of each object, the object Ω is wrapped in the smallest cuboid parallel to the coordinate axis, and it is defined as: AABB(Ω) = [x min ,x max ] x [y min ,y max ] x [z min ,z max ]; where x min and x max denote the minimum and maximum bounds of the object Ω on the x-axis, respectively, y min and y max denote the minimum and maximum bounds of the object Ω on the y-axis, respectively, and z min and z max denote the minimum and maximum bounds of the object Ω on the z-axis, respectively. Collision determination is based on whether the AABB of two objects A and B exists in the interval overlap in three-axis directions: S7.2, path search is based on collision detection results, first using RRT algorithm; RRT randomly samples points q in continuous configuration space C rand , and finds the nearest node q to the randomly sampled point in the existing path tree near : Where ‖·‖ represents the Euclidean distance, Q is the set of all nodes in the current tree, and arg min is the position of the minimum value; Subsequently, a new node q is generated using a step limit function Steer new : Where η is the maximum step length; S7.3, introduce RRT-Connect algorithm, while growing from start q start and goal q goal each grow a random tree and alternate expanding and trying to connect; connection condition is that nodes in both trees and are within threshold ∈ conn : S7.4, introduce RRT* algorithm, RRT* extends a new node q new after, optimize the path by local re-wiring mechanism; S7.5, the path points are smoothed and converted into continuous and dynamically controllable joint angle trajectories by using quintic polynomial interpolation and B-spline curve interpolation.

9. The dental implant robot path generation method of claim 8, wherein, Step S7.4 specifically comprises the following steps: S7.4.1, search for a set of neighboring nodes Q within radius r n near :​ Q near = { q e Q | || q - q new || < r n}; wherein r n Dynamic adjustment with sampling number; S7.4.2, select the neighboring node q that minimizes the total cost min ∈ Q near as the parent of q new : where c(q start , q) denotes the cost of the path from the start to node q, and Cost(q, q new ) denotes the cost of the path from q to q new , typically the path length. S7.4.3, perform a re-wiring operation to the neighboring nodes, if a lower path cost is obtained by q new The parent node and path cost are updated.

Images