Automatic registration method for neurosurgery facial navigation point cloud

Abstract

The invention relates to an automatic registration method for a neurosurgical facial navigation point cloud. The method comprises the following steps: segmenting and extracting a head point cloud to obtain a front facial point cloud as a target point cloud; a track point cloud obtained by preprocessing a navigation track of face navigation is used as a source point cloud; performing automatic registration from the source point cloud to the target point cloud; the method comprises the following steps: firstly, through coarse registration, obtaining rotation translation homogeneous matrixes of N conditions with minimum Euclidean distance values of origin points of coordinates of two systems in various conditions that a source point cloud coordinate system is aligned with a target point cloud coordinate system; through fine registration, taking the N rotation translation homogeneous matrixes as initial transformation matrixes of source point clouds, and respectively performing M registration of point-to-point, point-to-surface and surface-to-surface relationships between the source point clouds and target point clouds and combinations of every two of the three relationships to obtain M * N rigid body rotation matrixes and Euclidean distances; and screening out a rigid body rotation matrix with the best robustness as an automatic registration result to be output. According to the invention, the efficiency and accuracy of facial navigation point cloud registration are improved.

Description

Technical Field

[0001] This invention relates to the fields of data processing technology, point cloud registration technology, and medical image processing, and particularly to an automatic registration method for point clouds used in neurosurgical facial navigation. Background Technology

[0002] To achieve matching between surgical navigation equipment and corresponding medical images of the operating area inside or outside the patient during neurosurgery, preoperative registration is required. Registration typically uses, for example, surface point clouds of the facial operating area and trajectory point clouds within that area obtained through the navigation equipment.

[0003] Imaging techniques such as magnetic resonance imaging (MRI), computed tomography (CT), and ultrasound typically acquire images of the entire head, used for procedures like facial reconstruction surgery or facial feature analysis. Clinicians and researchers are more focused on data from the anterior face region. This is because the anterior face contains a rich array of facial muscles and important facial features, such as the eye sockets, bridge of the nose, and lips, which are crucial for facial aesthetics and functional restoration. Complete head point cloud data contains a significant amount of additional information beyond the anterior face. This redundant data not only increases data processing time but can also lead to unnecessary waste of computational resources and even affect the efficiency and accuracy of data analysis.

[0004] Furthermore, the trajectories obtained through navigation devices are typically drawn manually by doctors. Since the movement speed during manual drawing is not constant, the point spacing in the trajectory point cloud is not equal. Moreover, due to the low data output frequency and large point spacing, further processing of the obtained trajectory point cloud is necessary to use more corresponding points for matching in the registration process to obtain more accurate registration results.

[0005] Furthermore, the Iterative Closest Point (ICP) algorithm currently used for 3D point cloud matching is prone to getting trapped in local optima and stopping iteration before finding a global optimum. ICP adds a three-point coarse registration step to avoid local optima, but this step cannot be performed automatically and requires manual selection of points on the target point cloud, making the process complex. Moreover, it suffers from low registration success rates for irregularly shaped trajectories. Summary of the Invention

[0006] Based on the above analysis, the present invention aims to disclose an automatic registration method for neurosurgical facial navigation point clouds; and to solve the problem of automatic registration between neurosurgical facial navigation point clouds and target point clouds.

[0007] This invention discloses an automatic registration method for neurosurgical facial navigation point clouds, comprising:

[0008] Step S1: Segment and extract the head point cloud to obtain the front head point cloud as the target point cloud;

[0009] Step S2: The trajectory point cloud obtained by preprocessing the navigation trajectory of the facial navigation is used as the source point cloud;

[0010] Step S3: Perform automatic registration from the source point cloud to the target point cloud;

[0011] In the automatic registration process, coarse registration is first used to obtain the rotation and translation homogeneous matrices for N cases where the Euclidean distance between the origins of the two coordinate systems is minimized among various cases where the source point cloud coordinate system and the target point cloud coordinate system are aligned. Then, fine registration is used to take the N rotation and translation homogeneous matrices from the coarse registration as the initial transformation matrices of the source point cloud. M registrations are then performed between the source point cloud and the target point cloud, including point-to-point, point-to-plane, and plane-to-plane relationships, as well as pairwise combinations of the above three relationships, to obtain M*N rigid body rotation matrices and Euclidean distances. The rigid body rotation matrix with the best robustness is then selected as the automatic registration result output.

[0012] Further, step S1 includes:

[0013] Step S101: Solve for the curvature values ​​of the points in the head point cloud to obtain the curvature distribution data of the head point cloud; filter out the curvature maxima in the head point cloud to obtain the set of curvature maxima.

[0014] Step S102: Perform clustering operation on the set of points with maximum curvature to obtain 6 point clusters including the ears, nose, eye sockets and lips;

[0015] Step S103: Based on the centroid coordinates of the ear and nose point clusters, construct a section to divide the head point cloud, divide the head point cloud into two parts, and extract the point cloud data containing the point clusters to obtain the front head point cloud.

[0016] The cut surface is perpendicular to the XOY plane of the head coordinate system, and intersects the centroid of the isosceles triangle formed by the projections of the centroids of the two ears and the nose on the XOY plane, and is parallel to the line connecting the centroids of the two ears.

[0017] Furthermore, the method for obtaining the curvature distribution data of the head point cloud includes:

[0018] Step S101-1: Select the neighboring points around each point in the head point cloud; establish a set N of neighboring points. p ;

[0019] Step S101-2: In the neighborhood point set N p The neighborhood points are used to fit a local plane;

[0020] Step S101-3: Construct the covariance matrix C of the neighborhood points using the coordinates of the neighborhood points in the local plane;

[0021] Step S101-4: Calculate the eigenvalues ​​and eigenvectors of the covariance matrix C;

[0022] Step S101-5: Sort the eigenvalues ​​of the covariance matrix C of each neighboring point, extract the maximum and minimum values ​​of the eigenvalues, calculate the Gaussian curvature of each neighboring point, and obtain the curvature distribution data of the head point cloud.

[0023] Further, in step S2, the navigation trajectory is the navigation trajectory obtained by the electromagnetic navigation registration device through electromagnetic positioning when the surgeon operates on the anterior part using surgical instruments including an electromagnetic navigation registration device; the preprocessing is equidistant interpolation processing.

[0024] Furthermore, the equidistant interpolation process in step S2 includes:

[0025] Step S201: Perform interpolation initialization;

[0026] Determine the number of interpolation points N and the interpolation point spacing intv; set the trajectory between each pair of adjacent points of the input navigation trajectory as an interpolation segment, divide the trajectory of M points into M-1 segments, set the trajectory start point as the segment interpolation start point of the first interpolation segment, and the trajectory points adjacent to the trajectory start point as the segment end point of the first interpolation segment; initialize the segment number i = 1;

[0027] Step S202: Determine whether segment number i = M is true; if yes, then the equidistant interpolation has reached the end of the trajectory and the interpolation is stopped; if no, then calculate the segment direction vector X = {X(1), X(2), X(3)} of the current interpolation segment with segment number i; calculate the Euclidean distance inter1 between the segment interpolation start point and the segment end point of the current interpolation segment.

[0028] Step S203: Determine whether inter1≥intv is true. If true, proceed to step S204; otherwise, proceed to step S207. Intv is the interpolation point spacing.

[0029] Step S204: In the segment direction of the current interpolation segment, perform intra-segment equidistant interpolation successively from the segment interpolation start point until the segment interpolation end point;

[0030] Step S205: Determine whether the segment interpolation endpoint coincides with the segment endpoint of the current interpolation segment; if yes, determine the segment interpolation starting point of the next segment as the segment starting point; if no, search for a point in the next segment that is at a distance of intv-dist1 from the current segment endpoint as the segment interpolation starting point of the next segment; dist1 is the distance from the segment interpolation endpoint of the current segment to the current segment endpoint; the current segment endpoint is the segment starting point of the next segment.

[0031] Step S206: After determining the starting point of the next segment interpolation, set i = i + 1 and return to step S202;

[0032] Step S207: Do not interpolate the current segment; search for a point in the next segment that is at a distance of intv-dist2 from the end point of the current segment as the starting point for segment interpolation of the next segment; dist2 is the length of the current segment; the end point of the current segment is the starting point of the next segment.

[0033] Step S208: After determining the starting point of the next segment interpolation, set i = i + 1 and return to step S202;

[0034] Interpolated data is stored in a storage array of length N. The first position of the storage array stores the starting point of the input navigation trajectory. Starting from the second position, the interpolated points generated in step S204 are stored sequentially until the storage array is full, thus forming a source point cloud.

[0035] Further, step S3 includes:

[0036] Step S301: Perform coarse registration; perform principal component analysis on the target point cloud and the source point cloud respectively to obtain the coordinate systems of the target point cloud and the source point cloud. Using the coordinate system of the target point cloud as the reference, adjust the pointing direction of each axis of the source point cloud coordinate system in turn to align it with the coordinate system of the target point cloud; calculate each alignment case and output the rotation and translation homogeneous matrix of the N cases with the smallest Euclidean distance between the origin of the source point cloud and the target point cloud during alignment.

[0037] Step S302: Perform fine registration; use the N homogeneous rotation and translation matrices from the coarse registration as the initial transformation matrices of the source point cloud, and perform M types of registration between the source point cloud and the target point cloud, including point-to-point, point-to-surface, and surface-to-surface relationships, as well as pairwise combinations of the above three relationships, to obtain M*N rigid body rotation matrices and Euclidean distances. Select the rigid body rotation matrix with the best robustness as the fine registration result output.

[0038] Furthermore, the coarse registration process in step S301 includes:

[0039] Step S301-1: Perform principal component analysis on the target point cloud and the source point cloud respectively to obtain the coordinate systems of the target point cloud and the source point cloud;

[0040] Step S301-2: Using the coordinate system of the target point cloud as a reference, adjust the pointing directions of the +X, -X, +Y, -Y, +Z and -Z axes of the source point cloud coordinate system in sequence, and form various adjustment situations through permutation and combination;

[0041] Step S301-3: Calculate the rotation and translation homogeneous matrices corresponding to various adjustment cases using the three-point registration method; and eliminate cases where the coordinate systems cannot be aligned after transformation, and select the adjustment cases that can be aligned.

[0042] Step S301-4: Input the source point cloud and the target point cloud, select the rotation and translation homogeneous matrices corresponding to the alignment adjustment cases as the initial transformation, use KD-tree to find the corresponding points of the source point cloud and the target point cloud, calculate the rigid body transformation matrix by SVD decomposition, and update the position of the iterative trajectory point cloud.

[0043] Step S301-5: After satisfying the iteration limit, record the Euclidean distance value (fitness_score) between the origin of the source point cloud and the target point cloud in various alignment adjustment situations.

[0044] Step S301-6: Sort the Euclidean distance values ​​(fitness_score) between the origin and the source point cloud and the target point cloud, and output the rotation and translation homogeneous matrices for the N cases with the smallest fitness_score values.

[0045] Furthermore, the methods for determining principal component elements include:

[0046] 1) Point cloud data centralization; First, calculate the average value of the point cloud, and then subtract the average value from each point in the point cloud;

[0047] 2) Construct the covariance matrix; use the centered dataset to calculate the covariance matrix C;

[0048] 3) Calculate the eigenvalues ​​and eigenvectors; solve for the eigenvalues ​​and corresponding eigenvectors of the covariance matrix C to obtain the principal component vectors;

[0049] 4) Sort the feature values; arrange the feature vectors in descending order of feature value; the direction corresponding to the largest feature value is the direction of the greatest data change;

[0050] 5) Determine the coordinate system; select the eigenvectors corresponding to the first 3 largest eigenvalues ​​as the 3 principal component vectors of the point cloud.

[0051] Furthermore, the fine registration process in step S302 includes:

[0052] Step S302-1: Using the rotation and translation homogeneous matrices of the N cases output in step S301, perform initial transformations on the source point cloud to obtain N initial source point clouds.

[0053] Step S302-2: Using N initial source point clouds, perform a total of 7 registrations between the source point cloud and the target point cloud, including point-to-point, point-to-surface, surface-to-surface, point-to-point followed by surface-to-surface, point-to-surface followed by surface-to-surface, surface-to-surface followed by point-to-point, and surface-to-surface followed by point-to-surface, to obtain 7*N rigid body rotation matrices and Euclidean distances.

[0054] Step S302-3: Select the rigid body rotation matrix with the best robustness as the fine registration result output.

[0055] Furthermore, in step S302-2, during the registration process between the source point cloud and the target point cloud, first as a point-to-point and then as a surface-to-surface operation,

[0056] First, point-to-point registration is performed. Using the obtained rigid body rotation matrix R and translation matrix t, the initial source point cloud is transformed to obtain point cloud data, which is then used as the initial point cloud for subsequent face-to-face registration. This process yields the registered rigid body rotation matrix and the final Euclidean distance value, fitness_score.

[0057] In step S302-2, during the registration process between the source point cloud and the target point cloud, first between points and surfaces, then between surfaces,

[0058] First, point-to-surface registration is performed. The obtained rigid body rotation matrix R and translation matrix t are used to transform the initial source point cloud to obtain point cloud data, which is then used as the initial point cloud for subsequent face-to-surface registration. This process yields the registered rigid body rotation matrix and the final Euclidean distance value, fitness_score.

[0059] In step S302-2, during the registration process between the source point cloud and the target point cloud, first between surfaces and then between points,...

[0060] First, perform face-to-face registration. Then, use the obtained rigid body rotation matrix R and translation matrix t to transform the initial source point cloud to obtain point cloud data, which is used as the initial point cloud for subsequent point-to-point registration. This process yields the registered rigid body rotation matrix and the final Euclidean distance value, fitness_score.

[0061] In step S302-2, during the registration process between the source point cloud and the target point cloud, first between surfaces and then between points and surfaces,

[0062] First, perform face-to-face registration. Then, use the obtained rigid body rotation matrix R and translation matrix t to transform the initial source point cloud to obtain point cloud data, which is used as the initial point cloud for subsequent point-to-face registration. This process yields the registered rigid body rotation matrix and the final Euclidean distance value, fitness_score.

[0063] The present invention can achieve the following beneficial effects:

[0064] 1. This invention provides an efficient, accurate, and reliable frontal point cloud extraction scheme by combining curvature analysis and clustering. Curvature analysis can capture subtle changes on the head surface, thereby enabling accurate identification and extraction of the frontal feature point set in the head point cloud. The automated calculation process reduces manual intervention and improves data processing efficiency. Furthermore, data verification ensures the accuracy of the segmentation operation.

[0065] 2. This invention achieves uniform spatial distribution and enhanced density of point cloud data by performing equidistant interpolation on the trajectory point cloud. This improvement effectively solves the problem of uneven point spacing caused by inconsistent operating speeds of doctors in traditional hand-drawn trajectory point clouds, thereby improving the overall quality of the point cloud data. Uniform and dense point cloud data not only provides richer geometric information but also significantly improves the accuracy and stability of subsequent registration processes.

[0066] 3. The automatic point cloud registration of this invention can automatically process 3D point cloud data, reduce manual intervention, and improve the efficiency and accuracy of registration. Furthermore, by combining multiple registration strategies, the robustness of the registration results can be improved, and the applicability of the registration method can be expanded. Attached Figure Description

[0067] The accompanying drawings are for illustrative purposes only and are not intended to limit the invention. Throughout the drawings, the same reference numerals denote the same parts.

[0068] Figure 1 This is a flowchart of the automatic registration method for neurosurgical facial navigation point clouds in an embodiment of the present invention. Detailed Implementation

[0069] Preferred embodiments of the present invention will now be described in detail with reference to the accompanying drawings, which form part of this application and, together with the embodiments of the present invention, serve to illustrate the principles of the present invention.

[0070] One embodiment of the present invention discloses an automatic registration method for neurosurgical facial navigation point clouds, such as... Figure 1 As shown, it includes:

[0071] Step S1: Segment and extract the head point cloud to obtain the front head point cloud as the target point cloud;

[0072] Step S2: The trajectory point cloud obtained by preprocessing the navigation trajectory of the facial navigation is used as the source point cloud;

[0073] Step S3: Perform automatic registration from the source point cloud to the target point cloud;

[0074] In the automatic registration process, coarse registration is first used to obtain the rotation and translation homogeneous matrices for N cases where the Euclidean distance between the origins of the two coordinate systems is minimized among various cases where the source point cloud coordinate system and the target point cloud coordinate system are aligned. Then, fine registration is used to take the N rotation and translation homogeneous matrices from the coarse registration as the initial transformation matrices of the source point cloud. M registrations are then performed between the source point cloud and the target point cloud, including point-to-point, point-to-plane, and plane-to-plane relationships, as well as pairwise combinations of the above three relationships, to obtain M*N rigid body rotation matrices and Euclidean distances. The rigid body rotation matrix with the best robustness is then selected as the automatic registration result output.

[0075] Specifically, in step S1, the head point cloud is obtained by extracting the surface contours of the stacked DICOM data to obtain head point cloud data composed of surface contour data.

[0076] The head coordinate system follows the right-hand screw rule, and each coordinate axis is defined as follows: the positive X-axis points from the left side of the patient's head to the right side; the positive Y-axis points from the back of the head to the front of the head; and the positive Z-axis points from the lower jaw to the top of the head.

[0077] DICOM data is obtained by scanning the human head using medical devices, including MRI and CT scans.

[0078] In step S1, the head point cloud is segmented and extracted to obtain the frontal point cloud as the target point cloud; including:

[0079] Step S101: Solve for the curvature values ​​of the points in the head point cloud to obtain the curvature distribution data of the head point cloud; filter out the curvature maxima in the head point cloud to obtain the set of curvature maxima.

[0080] Step S102: Perform clustering operation on the set of points with maximum curvature to obtain 6 point clusters including the ears, nose, eye sockets and lips;

[0081] Step S103: Based on the centroid coordinates of the ear and nose point clusters, construct a section to divide the head point cloud, divide the head point cloud into two parts, and extract the point cloud data containing the point clusters to obtain the front head point cloud.

[0082] The cut surface is perpendicular to the XOY plane of the head coordinate system, and intersects the centroid of the isosceles triangle formed by the projections of the centroids of the two ears and the nose on the XOY plane, and is parallel to the line connecting the centroids of the two ears.

[0083] In step S101, the curvature values ​​of the points in the head point cloud are calculated to obtain the curvature distribution data of the head point cloud, including:

[0084] 1) Select the neighboring points around each point in the head point cloud; establish a set N of neighboring points. p ;

[0085] The neighboring points are at point P i Points within a fixed radius 'a' around the cloud, where the value of 'a' can be determined based on the specific point cloud density, and in this embodiment, the value of 'a' is 2 mm.

[0086] 2) In the neighborhood point set N p The neighborhood points are used to fit a local plane;

[0087] The fitted plane equation is: ax + by + cz + d = 0; where a, b, and c are the fitting coefficients, and x, y, and z are the plane coordinates.

[0088] The preferred fitting method is the least squares method, which aims to minimize the sum of squares of the distances from all neighboring points to the fitting plane.

[0089] 3) Construct the covariance matrix C of the neighborhood points using the coordinates of the neighborhood points in the local plane;

[0090]

[0091] Where, p i Let be the coordinates of the neighboring point i in the local plane. Let be the average coordinates of the neighboring points, and n be the number of neighboring points.

[0092] 4) Calculate the eigenvalues ​​and eigenvectors of the covariance matrix C;

[0093] det(C-λ i I) = 0;

[0094] (C-λ i,k I)V i,k =0;

[0095] Where, λ i,k V represents the eigenvalues ​​of the covariance matrix of neighborhood point i. k Let be the eigenvector of the covariance matrix of neighborhood point i, k = 1, 2, 3.

[0096] 5) Sort the eigenvalues ​​of the covariance matrix C of each neighboring point, extract the maximum and minimum eigenvalues, calculate the Gaussian curvature of each neighboring point, and obtain the curvature distribution data of the head point cloud;

[0097] Gaussian curvature K of neighborhood point i i :

[0098] K i =λ i,min ·λ i,max ;

[0099] Where, λ i,min λ is the smallest eigenvalue of the covariance matrix of neighborhood point i.i,max Let be the largest eigenvalue of the covariance matrix of the neighborhood point i.

[0100] In step S101, the curvature maxima in the head point cloud are selected to obtain a set of curvature maxima.

[0101] 1) Calculate the mean curvature μ of the head point cloud:

[0102]

[0103] 2) Calculate the variance σ of the curvature values ​​of the head point cloud. 2 ;

[0104]

[0105] 3) Select points with curvature values ​​greater than 8 times the variance to form a set of curvature maxima;

[0106] Point set N q Medium curvature value K i ≥8σ 2 +μ.

[0107] Step S102 includes:

[0108] Step S102-1: Randomly select 6 points C1, C2, C3, C4, C5, and C6 as the initial cluster centers;

[0109] Step S102-2: Calculate the distance from each point in the set of curvature maxima to the center of each cluster, and assign the point to the nearest cluster C. j ;

[0110] The point of maximum curvature p is at a distance from the point cluster C j The distance is:

[0111]

[0112] Among them, (x p ,y p ,z p Let ) be the coordinates of the point p, which is the maximum curvature. For point cluster C j Center coordinates, j = 1, 2, ..., 6;

[0113] Step S102-3: Calculate the average value of the points within each cluster, and use it as the new cluster center;

[0114] New cluster center C′ j for:

[0115]

[0116] Step S102-4: Repeat steps S102-2 and S102-3 until the change in cluster center is less than the threshold, and output the 6 clusters after clustering.

[0117] The threshold is set to 0.0002 mm, the point spacing of the head point cloud is about 0.2 mm, and the iteration stops when the movement distance of the cluster center is less than one-thousandth of the point spacing.

[0118] After obtaining 6 clusters from the clustering, the clusters are divided according to their centroids. The two clusters with centroids close to the outer contour and symmetrical relationship are identified as the ears, the cluster with the largest mean curvature is identified as the nose, the two clusters with centroids far from the outer contour and symmetrical relationship are identified as the eye sockets, the remaining clusters are the lips, and the area where the clusters are clustered is identified as the front face.

[0119] Step S103 includes:

[0120] Step S103-1: Take the coordinate axis that is closer to the centroid of the ear and nose point clusters as the Z-axis of the head coordinate system; project the centroid of the ear and nose point clusters onto the XOY plane of the coordinate system to obtain an isosceles triangle.

[0121] The centroid of the nose and the centroids of both ears lie roughly in the same plane. Taking the closest coordinate axis among these three points as the z-axis, this means that the analysis can be performed on the XOY plane formed by the X and Y axes. Based on the symmetry of the ears, the projections of the three points onto the XOY plane can form an isosceles triangle.

[0122] Step S103-2: Draw a straight line on the XOY plane that passes through the centroid of the triangle and is parallel to the line connecting the centroids of the two ears; and draw a cross-section perpendicular to the XOY plane, which intersects the XOY plane at the straight line; the cross-section divides the head point cloud into two parts; the part containing the front face where the point clusters are gathered is taken out for subsequent operations.

[0123] Preferably, step S1 further includes a verification step:

[0124] 1) Filter out the first verification point with the largest Y-axis in the complete head point cloud, and the second verification point with the largest Y-axis in the segmented front head point cloud;

[0125] 2) If the difference between the coordinates of the first verification point and the second verification point is within the error range, then the segmentation of the front point cloud is accurate; otherwise, it is inaccurate.

[0126] The error range is determined based on the accuracy of the point cloud.

[0127] Step S1 involves segmenting and extracting the head point cloud to obtain the frontal point cloud as the target point cloud. This target point cloud contains most of the feature information of the patient's head, and the trajectory is drawn only on the front part. Using only the frontal point cloud for subsequent operations can reduce the amount of computation and also reduce the influence of other irrelevant points on the registration.

[0128] Specifically, in step S2, the navigation trajectory is the navigation trajectory obtained by the electromagnetic navigation registration device through electromagnetic positioning when the surgeon operates on the anterior part using surgical instruments including an electromagnetic navigation registration device; the preprocessing is equidistant interpolation processing.

[0129] The equidistant interpolation process in step S2 includes:

[0130] Step S201: Perform interpolation initialization;

[0131] Determine the number of interpolation points N and the interpolation point spacing intv; set the trajectory between each pair of adjacent points of the input navigation trajectory as an interpolation segment, divide the trajectory of M points into M-1 segments, set the trajectory start point as the segment interpolation start point of the first interpolation segment, and the trajectory points adjacent to the trajectory start point as the segment end point of the first interpolation segment; initialize the segment number i = 1;

[0132] Wherein, the interpolation point spacing intv = dist / (N-1); dist is the total distance after accumulating the Euclidean distances between adjacent points of the input trajectory; the interpolation point spacing is not greater than 0.2mm.

[0133] Step S202: Determine whether segment number i = M is true; if yes, then the equidistant interpolation has reached the end of the trajectory and the interpolation is stopped; if no, then calculate the segment direction vector X = {X(1), X(2), X(3)} of the current interpolation segment with segment number i; calculate the Euclidean distance inter1 between the segment interpolation start point and the segment end point of the current interpolation segment.

[0134] Specifically, the segment direction vector X = nextpt – curpt; where curpt is the coordinate vector of the starting point of the current interpolation segment, and nextpt is the coordinate vector of the ending point of the current interpolation segment;

[0135] Step S203: Determine whether inter1≥intv is true. If true, proceed to step S204; otherwise, proceed to step S207. Intv is the interpolation point spacing.

[0136] Step S204: In the segment direction of the current interpolation segment, perform intra-segment equidistant interpolation successively from the segment interpolation start point until the segment interpolation end point;

[0137] Step S205: Determine whether the segment interpolation endpoint coincides with the segment endpoint of the current interpolation segment; if yes, determine the segment interpolation starting point of the next segment as the segment starting point; if no, search for a point in the next segment that is at a distance of intv-dist1 from the current segment endpoint as the segment interpolation starting point of the next segment; dist1 is the distance from the segment interpolation endpoint of the current segment to the current segment endpoint; the current segment endpoint is the segment starting point of the next segment.

[0138] Step S206: After determining the starting point of the next segment interpolation, set i = i + 1 and return to step S202;

[0139] Step S207: Do not interpolate the current segment; search for a point in the next segment that is at a distance of intv-dist2 from the end point of the current segment as the starting point for segment interpolation of the next segment; dist2 is the length of the current segment; the end point of the current segment is the starting point of the next segment.

[0140] Step S208: After determining the starting point of the next segment interpolation, set i = i + 1 and return to step S202;

[0141] Interpolated data is stored in a storage array of length N. The first position of the storage array stores the starting point of the input navigation trajectory. Starting from the second position, the interpolated points generated in step S204 are stored sequentially until the storage array is full, thus forming a source point cloud.

[0142] Specifically, in step S204, the equidistant interpolation includes:

[0143] 1) Find a point on the segment direction vector X that is at a distance of the interpolation point spacing intv from the segment interpolation starting point, and use it as the second interpolation point for segment interpolation;

[0144] 2) Starting from the second interpolation point, determine whether the distance between the current interpolation point and the end point of this segment is less than or equal to the interpolation point spacing intv; if not, proceed to 3); if yes, proceed to 4).

[0145] 3) Find the interpolation point on the segment direction vector X; find the point on the segment direction vector X that is at a distance of the interpolation point spacing intv from the current interpolation point, and use it as the next interpolation point for segment interpolation; use the next interpolation point as the current interpolation point and return to step 2) to continue the distance judgment;

[0146] Find the interpolation point k on the segment direction vector X i The formula is:

[0147] k i =[intv*l+k i-1 (1),intv*m+k i-1 (2), intv*n+k i-1 (3)]

[0148] Where, l=x(1),m=x(2),n=x(3); x={x(1),x(2),x(3)} are the three direction vectors after normalizing the segment direction vector X={X(1),X(2),X(3)}; the normalization method is: x=X / norm(X); i=2,3,…,{k1(1),k1(2),k1(3)} are the coordinates of the three directions of the segment interpolation starting point.

[0149] 4) Determine whether the distance between the current interpolation point and the end point of the current segment is equal to or less than the interpolation point spacing intv; if equal, then the end point of the segment is taken as the end point of the segment interpolation; if less than, then the current interpolation point is taken as the end point of the segment interpolation.

[0150] By performing equidistant interpolation on the navigation trajectory of the facial navigation in step S2, the uniformity and density of the trajectory point cloud, which serves as the source point cloud, are achieved. This provides a high-quality input point cloud for the registration in step S3, avoiding the adverse effects on registration caused by inaccurate calculation of normal and tangent vectors due to the unevenness of the trajectory point cloud.

[0151] Specifically, step S3 includes:

[0152] Step S301: Perform coarse registration; perform principal component analysis on the target point cloud and the source point cloud respectively to obtain the coordinate systems of the target point cloud and the source point cloud. Using the coordinate system of the target point cloud as the reference, adjust the pointing direction of each axis of the source point cloud coordinate system in turn to align it with the coordinate system of the target point cloud; calculate each alignment case and output the rotation and translation homogeneous matrix of the N cases with the smallest Euclidean distance between the origin of the source point cloud and the target point cloud during alignment.

[0153] Step S302: Perform fine registration; use the N homogeneous rotation and translation matrices from the coarse registration as the initial transformation matrices of the source point cloud, and perform M types of registration between the source point cloud and the target point cloud, including point-to-point, point-to-surface, and surface-to-surface relationships, as well as pairwise combinations of the above three relationships, to obtain M*N rigid body rotation matrices and Euclidean distances. Select the rigid body rotation matrix with the best robustness as the fine registration result output.

[0154] Specifically, the coarse registration process in step S301 includes:

[0155] Step S301-1: Perform principal component analysis on the target point cloud and the source point cloud respectively to obtain the coordinate systems of the target point cloud and the source point cloud;

[0156] Step S301-2: Using the coordinate system of the target point cloud as a reference, adjust the pointing directions of the +X, -X, +Y, -Y, +Z and -Z axes of the source point cloud coordinate system in sequence, and form various adjustment situations through permutation and combination;

[0157] Step S301-3: Calculate the rotation and translation homogeneous matrices corresponding to various adjustment cases using the three-point registration method; and eliminate cases where the coordinate systems cannot be aligned after transformation, and select the adjustment cases that can be aligned.

[0158] Step S301-4: Input the source point cloud and the target point cloud, select the rotation and translation homogeneous matrices corresponding to the alignment adjustment cases as the initial transformation, use KD-tree to find the corresponding points of the source point cloud and the target point cloud, calculate the rigid body transformation matrix by SVD decomposition, and update the position of the iterative trajectory point cloud.

[0159] Step S301-5: After satisfying the iteration limit, record the Euclidean distance value (fitness_score) between the origin of the source point cloud and the target point cloud in various alignment adjustment situations.

[0160] Step S301-6: Sort the Euclidean distance values ​​(fitness_score) between the origin and the source point cloud and the target point cloud, and output the rotation and translation homogeneous matrices for the N cases with the smallest fitness_score values.

[0161] Specifically, step S301-1 includes:

[0162] Step S301-1-1: Perform principal component analysis on the target point cloud and the source point cloud to obtain their respective three principal component vectors;

[0163] Specifically, including:

[0164] 1) Point cloud data centralization; First, calculate the average value of the point cloud, and then subtract the average value from each point in the point cloud;

[0165] Let the point cloud dataset be X = {X1, ..., X...} i ,…,X N}, where each point X i Both are three-dimensional vectors X i =(x i ,y i ,z i );

[0166] Calculate the average position of all points

[0167]

[0168] Then, the data is centralized to obtain centralized data X′={X′1,…,X′ i ,…X′ N}

[0169]

[0170] 2) Construct the covariance matrix; use the centered dataset to calculate the covariance matrix C;

[0171] The covariance matrix can be represented as

[0172] 3) Calculate the eigenvalues ​​and eigenvectors; solve for the eigenvalues ​​λ of the covariance matrix C. i and the corresponding feature vector v i This yields the principal component vectors;

[0173] Cv = λv;

[0174] 4) Sort the feature values; arrange the feature vectors in descending order of feature value; the direction corresponding to the largest feature value is the direction of the greatest data change;

[0175] 5) Determine the coordinate system; select the eigenvectors corresponding to the first 3 largest eigenvalues ​​as the 3 principal component vectors of the point cloud.

[0176] Step S301-1-2: Calculate the centroids of the target point cloud and the source point cloud;

[0177] The formula for the centroid of a point cloud with N points is:

[0178]

[0179] Step S301-1-3: Based on the directions of the three principal component vectors, and with the centroid of each point cloud as the origin, determine the coordinate systems of the target point cloud and the source point cloud respectively using the right-hand screw rule.

[0180] Specifically, in step S301-2, the coordinate system of the target point cloud is used as a reference, and the pointing directions of the +X, -X, +Y, -Y, +Z and -Z axes of the source point cloud coordinate system are adjusted in sequence, forming 2*2*2*6=48 adjustment situations through permutation and combination.

[0181] Specifically, step S301-3 includes:

[0182] Step S301-3-1: Determine three points in each of the target point cloud and source point cloud coordinate systems as three-point registration points based on the origin coordinates and unit direction vector of the target point cloud and source point cloud coordinate systems.

[0183] The unit direction vectors of the three coordinate axes of the target point cloud are (x1, y1, z1), (x2, y2, z2), and (x3, y3, z3), respectively, and its origin coordinates are (O1, O2, O3). Therefore, the three points used for three-point registration are:

[0184] q1(O1+x1,O2+y1,O3+z1),

[0185] q2(O1+x2,O2+y2,O3+z2),

[0186] q3(O1+x3,O2+y3,O3+z3);

[0187] The unit direction vectors of the three coordinate axes of the source point cloud are (x′1, y′1, z′1), (x′2, y′2, z′2), and (x′3, y′3, z′3), respectively, and the coordinates of its origin are (O′1, O′2, O′3). Therefore, the three points used for three-point registration are:

[0188] p1(O′1+x′1,O′2+y′1,O′3+z′1),

[0189] p2(O′1+x′2,O′2+y′2,O′3+z′2),

[0190] p3(O′1+x′3,O′2+y′3,O′3+z′3).

[0191] Step S301-3-2: Perform three-point registration on three points in each of the target point cloud and source point cloud coordinate systems to obtain the rigid body transformation matrix R and translation matrix t;

[0192] The three-point registration method includes:

[0193] 1) Calculate the centroids of the source point cloud and the target point cloud respectively;

[0194] Centroid of Source Cloud

[0195] centroid of the target point cloud

[0196] 2) Centralized point set;

[0197] Centralized source cloud midpoint

[0198] Centralized target point cloud midpoint

[0199] 3) Construct the rigid body rotation matrix; In order to find the rigid body rotation matrix R, two 3x3 matrices A and B need to be constructed and the singular value (SVD) solution needs to be obtained.

[0200] Matrices A and B are composed of the centered source point set and target point set, respectively;

[0201]

[0202] The outer product matrix H of A and B:

[0203] H = BA T ;

[0204] Perform singular value decomposition on H:

[0205] H = USV T ;

[0206] The rigid body rotation matrix R can be calculated using U and V (W is a 3x3 identity matrix):

[0207] R = UWV T ;

[0208] 4) Calculate the translation matrix t;

[0209]

[0210] Step S301-3-3: Perform rotation and translation transformations on the three points in the source point cloud coordinate system according to the rigid body transformation matrix R and translation matrix t to obtain new positions of the three points; calculate the Euclidean distance between these three points and the corresponding three points in the target point cloud coordinate system; if the calculated distance value is greater than the distance threshold, it is determined that it cannot be aligned and is eliminated, and the alignable adjustment cases are selected.

[0211] The distance threshold is 1e-2. In this embodiment, the input trajectory point cloud and target point cloud are used to filter out 24 alignable adjustment cases, and their corresponding transformation matrices are used for initialization transformation.

[0212] In step S301-4, the rotation and translation homogeneous matrices corresponding to the 24 alignable adjustment cases are selected as the initial transformation. The KD-tree is used to find the corresponding points of the source point cloud and the target point cloud. The rigid body transformation matrix is ​​calculated by SVD decomposition and the position of the iterative trajectory point cloud is updated.

[0213] In step S301-5, the number of iterations is set to 2. After 2 iterations, the Euclidean distance value of the origin of the coordinates in the source point cloud and the target point cloud in the alignment is recorded for 24 alignment adjustment cases.

[0214] In steps S301-6, N is 3, that is, the rotation and translation homogeneous matrices of the three cases with the smallest fitness_score values ​​are selected and output.

[0215] In the coarse registration step S301, three good initial positions were obtained that can be used for the subsequent precise registration step. After coarse registration, the trajectory has been transformed to a position that is relatively close to the target. Next, fine registration is performed in step S302.

[0216] Specifically, step S302 includes:

[0217] Step S302-1: Using the rotation and translation homogeneous matrices of the N cases output in step S301, perform initial transformations on the source point cloud to obtain N initial source point clouds.

[0218] The initial source point cloud obtained from the initial transformation is trace init =R*trace+t;

[0219] Trace is the source point cloud, and R and t are the rigid body rotation and translation matrices in the homogeneous rotation and translation matrices for one of the cases.

[0220] Step S302-2: Using N initial source point clouds, perform a total of 7 registrations between the source point cloud and the target point cloud, including point-to-point, point-to-surface, surface-to-surface, point-to-point followed by surface-to-surface, point-to-surface followed by surface-to-surface, surface-to-surface followed by point-to-point, and surface-to-surface followed by point-to-surface, to obtain 7*N rigid body rotation matrices and Euclidean distances.

[0221] With N=3, a total of 21 rigid body rotation matrices and Euclidean distances were obtained.

[0222] Step S302-3: Select the rigid body rotation matrix with the best robustness as the fine registration result output.

[0223] In this embodiment, the average value of the Euclidean distance is used to evaluate robustness; the shorter the distance, the better the robustness.

[0224] Specifically, in step S302-2, the point-to-point registration process between the source point cloud and the target point cloud, using each initial source point cloud, includes:

[0225] 1) Calculate the pairwise Euclidean distance between points in the initial source point cloud and points in the target point cloud;

[0226] 2) Use KD-tree to find the point pairs with the shortest Euclidean distance between points in the initial source point cloud and points in the target point cloud, and form a point pair set;

[0227] 3) For all points in the initial source point cloud, calculate their covariance matrix to describe the transformation of the local structure around each point;

[0228] The covariance matrix describes the changes in the local structure around the point and reflects the uncertainty of the local region of the point cloud. It can be used as a weight matrix to reduce the influence of bad corresponding points.

[0229] 4) Set the objective function E(R,t) to minimize the weighted sum of squared distances of all point pairs in the point pair set; the distance is the distance from the initial source point cloud point to the target point cloud point after rotation and translation in the point pair.

[0230]

[0231] Among them, w ij s is an element in the inverse of the covariance matrix; i Let point S i The coordinates, t i Let point T j The coordinates are given by R, where R is the rigid body rotation matrix and t is the translation matrix.

[0232] 5) Use gradient descent to solve the objective function to obtain the optimal rigid body rotation matrix R and translation matrix t;

[0233] 6) Update the initial source point cloud position using the calculated rigid body rotation matrix R and translation matrix t, and then repeat the above steps until the termination condition is met; record the final Euclidean distance value of point-to-point registration, fitness_score;

[0234] The termination condition is when the number of iterations reaches a preset maximum value or the transformation amount is less than a threshold.

[0235] Specifically, in step S302-2, the point-area registration process between the source point cloud and the target point cloud using each initial source point cloud includes:

[0236] 1) Before registration, estimate the normal vectors of the points in the target point cloud; by finding the neighboring points of each point and fitting a local plane, the normal vector of the current point is obtained; where, the normal vector of point T in the target point cloud is... j The normal vector is n j ;

[0237] 2) Use KD-tree to find the point pairs with the shortest Euclidean distance between points in the initial source point cloud and points in the target point cloud, and form a point pair set;

[0238] 3) Using the coordinates of the initial source point cloud and the target point cloud in the point pair, and the normal vector of the target point cloud, calculate the distance from the initial source point cloud to the plane containing the target point cloud in the point pair;

[0239] Distance from the initial source point cloud points to the plane containing the target point cloud points

[0240] 4) For all points in the initial source point cloud, calculate their covariance matrix to describe the transformation of the local structure around each point;

[0241] The covariance matrix describes the changes in the local structure around the point and reflects the uncertainty of the local region of the point cloud. It can be used as a weight matrix to reduce the influence of bad corresponding points.

[0242] 5) Set the objective function E(R,t) to minimize the weighted sum of squared distances between points and opposite points in all point pairs in the point pair set; the distance between points and opposite points is the distance from the initial source point cloud point in the point pair after rotation and translation to the plane where the target point cloud point is located.

[0243]

[0244] Among them, w ij s is an element in the inverse of the covariance matrix; i Let point S i The coordinates, t i Let point T j The coordinates are given by R, where R is the rigid body rotation matrix and t is the translation matrix; n j Let point T j The normal vector;

[0245] 6) Use gradient descent to solve the objective function to obtain the optimal rigid body rotation matrix R and translation matrix t;

[0246] 7) Update the initial source point cloud position using the calculated rigid body rotation matrix R and translation matrix t, and then repeat the above steps until the termination condition is met; record the final Euclidean distance value of the point and surface registration, fitness_score;

[0247] The termination condition is when the number of iterations reaches a preset maximum value or the transformation amount is less than a threshold.

[0248] Specifically, in step S302-2, the surface-to-surface registration process between the source point cloud and the target point cloud, using each initial source point cloud, includes:

[0249] 1) Before registration, normal vectors are estimated for points in the initial source and target point clouds; the normal vector of the current point is obtained by finding the neighboring points of each point and fitting a local plane; where, in the initial source point cloud data, point T in the target point cloud... j The normal vector is n j ;

[0250] 2) Establish point S in the initial source point cloud i With plane Π in the target point cloud j The correspondence between them; the plane Π j From point T j and normal vector n j Define; calculate point S i Projected onto plane Π j point p′ on i Point S i and p′ i Establish point pairs to form point pair sets;

[0251] Point Si Projected onto plane Π j point p′ on i satisfy:

[0252]

[0253] Point S i Projected onto plane Π j point p′ on i It can be calculated using the following formula:

[0254]

[0255] 3) For all points in the initial source point cloud data, calculate their covariance matrix to describe the transformation of the local structure around each point;

[0256] The covariance matrix describes the changes in the local structure around the point and reflects the uncertainty of the local region of the point cloud. It can be used as a weight matrix to reduce the influence of bad corresponding points.

[0257] 4) Define the objective function E(R,t) to minimize the weighted sum of squared distances between points and their opposites in all pairs of points in the point pair set; the distance between the points and their opposites is the distance between points S in the point pair after rotation and translation. i to p′ i Plane Π j The distance;

[0258]

[0259] Among them, w ij s is an element in the inverse of the covariance matrix; i Let point S i The coordinates, p′ i Let point S i Projected onto plane Π j The coordinates on the coordinate system are given by R, where R is the rigid body rotation matrix and t is the translation matrix; n j Let point T j The normal vector;

[0260] 5) Use gradient descent to solve the objective function to obtain the optimal rigid body rotation matrix R and translation matrix t;

[0261] 6) Update the initial source point cloud position using the calculated rigid body rotation matrix R and translation matrix t, and then repeat the above steps until the termination condition is met; record the final Euclidean distance value fitness_score between the surfaces;

[0262] The termination condition is when the number of iterations reaches a preset maximum value or the transformation amount is less than a threshold.

[0263] Specifically, in step S302-2, each initial source point cloud is used to perform point-to-point and then surface-to-surface registration between the source point cloud and the target point cloud. First, point-to-point registration is performed. The obtained rigid body rotation matrix R and translation matrix t are used to transform the initial source point cloud to obtain point cloud data, which is then used as the initial point cloud for the subsequent surface-to-surface registration. This results in the registered rigid body rotation matrix and the final Euclidean distance value fitness_score.

[0264] Specifically, in step S302-2, each initial source point cloud is used to perform point-to-surface and then surface-to-surface registration between the source point cloud and the target point cloud. First, point-to-surface registration is performed. The obtained rigid body rotation matrix R and translation matrix t are used to transform the initial source point cloud to obtain point cloud data, which is then used as the initial point cloud for the subsequent surface-to-surface registration. This results in the registered rigid body rotation matrix and the final Euclidean distance value fitness_score.

[0265] Specifically, in step S302-2, each initial source point cloud is used to perform face-to-face and then point-to-point registration between the source point cloud and the target point cloud. First, face-to-face registration is performed. The obtained rigid body rotation matrix R and translation matrix t are used to transform the initial source point cloud to obtain point cloud data, which is then used as the initial point cloud for subsequent point-to-point registration. This results in the registered rigid body rotation matrix and the final Euclidean distance value fitness_score.

[0266] Specifically, in step S302-2, each initial source point cloud is used to perform face-to-face and then point-to-face registration between the source point cloud and the target point cloud. First, face-to-face registration is performed. The obtained rigid body rotation matrix R and translation matrix t are used to transform the initial source point cloud to obtain point cloud data, which is then used as the initial point cloud for subsequent point-to-face registration. This results in the registered rigid body rotation matrix and the final Euclidean distance value fitness_score.

[0267] In summary, the embodiments of the present invention can achieve the following beneficial effects:

[0268] 1. By combining curvature analysis and clustering, an efficient, accurate, and reliable solution for extracting frontal point clouds is provided. Curvature analysis can capture subtle changes on the head surface, thereby enabling accurate identification and extraction of frontal feature point sets from the head point cloud. The automated calculation process reduces manual intervention and improves data processing efficiency. Furthermore, data verification ensures the accuracy of the segmentation operation.

[0269] 2. By performing equidistant interpolation on the trajectory point cloud, a uniform spatial distribution and enhanced density of the point cloud data were achieved. This improvement effectively solves the problem of uneven point spacing caused by inconsistent operating speeds of doctors in traditional hand-drawn trajectory point clouds, thereby improving the overall quality of the point cloud data. Uniform and dense point cloud data not only provides richer geometric information but also significantly improves the accuracy and stability of subsequent registration processes.

[0270] 3. Automatic point cloud registration enables the automatic processing of 3D point cloud data, reducing manual intervention and improving registration efficiency and accuracy. Furthermore, combining multiple registration strategies enhances the robustness of the registration results and expands the applicability of the registration method.

[0271] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. An automatic registration method for neurosurgical facial navigation point clouds, characterized in that, include: Step S1: Segment and extract the head point cloud to obtain the front head point cloud as the target point cloud; Step S2: The trajectory point cloud obtained by preprocessing the navigation trajectory of the facial navigation is used as the source point cloud; Step S3: Perform automatic registration from the source point cloud to the target point cloud; In the automatic registration, coarse registration is first used to obtain the rotation and translation homogeneous matrices for N cases in which the Euclidean distance between the origins of the two coordinate systems is minimized among various cases in which the source point cloud coordinate system and the target point cloud coordinate system are aligned. Then, through fine registration, the N homogeneous rotation and translation matrices of the coarse registration are used as the initial transformation matrices of the source point cloud. M types of registration are then performed between the source point cloud and the target point cloud, including point-to-point, point-to-plane, and plane-to-plane relationships, as well as pairwise combinations of the above three relationships. This yields M*N rigid body rotation matrices and Euclidean distances. The rigid body rotation matrix with the best robustness is then selected as the automatic registration result output.

2. The automatic registration method for neurosurgical facial navigation point clouds according to claim 1, characterized in that, Step S1 includes: Step S101: Solve for the curvature values ​​of the points in the head point cloud to obtain the curvature distribution data of the head point cloud; filter out the curvature maxima in the head point cloud to obtain the set of curvature maxima. Step S102: Perform clustering operation on the set of points with maximum curvature to obtain 6 point clusters including the ears, nose, eye sockets and lips; Step S103: Based on the centroid coordinates of the ear and nose point clusters, construct a section to divide the head point cloud, divide the head point cloud into two parts, and extract the point cloud data containing the point clusters to obtain the front head point cloud. The cut surface is perpendicular to the XOY plane of the head coordinate system, and intersects the centroid of the isosceles triangle formed by the projections of the centroids of the two ears and the nose on the XOY plane, and is parallel to the line connecting the centroids of the two ears.

3. The automatic registration method for neurosurgical facial navigation point clouds according to claim 2, characterized in that, Methods for obtaining curvature distribution data of head point clouds include: Step S101-1: Select the neighboring points around each point in the head point cloud; establish a set N of neighboring points. p ; Step S101-2: In the neighborhood point set N p The neighborhood points are used to fit a local plane; Step S101-3: Construct the covariance matrix C of the neighborhood points using the coordinates of the neighborhood points in the local plane; Step S101-4: Calculate the eigenvalues ​​and eigenvectors of the covariance matrix C; Step S101-5: Sort the eigenvalues ​​of the covariance matrix C of each neighboring point, extract the maximum and minimum values ​​of the eigenvalues, calculate the Gaussian curvature of each neighboring point, and obtain the curvature distribution data of the head point cloud.

4. The automatic registration method for neurosurgical facial navigation point clouds according to claim 1, characterized in that, In step S2, the navigation trajectory is the navigation trajectory obtained by the electromagnetic navigation registration device through electromagnetic positioning when the surgeon operates on the anterior part using surgical instruments including an electromagnetic navigation registration device; the preprocessing is equidistant interpolation processing.

5. The automatic registration method for neurosurgical facial navigation point clouds according to claim 4, characterized in that, The equidistant interpolation process in step S2 includes: Step S201: Perform interpolation initialization; Determine the number of interpolation points N and the interpolation point spacing intv; set the trajectory between each pair of adjacent points of the input navigation trajectory as an interpolation segment, divide the trajectory of M points into M-1 segments, set the trajectory start point as the segment interpolation start point of the first interpolation segment, and the trajectory points adjacent to the trajectory start point as the segment end point of the first interpolation segment; initialize the segment number i = 1; Step S202: Determine whether segment number i = M is true; if yes, then the equidistant interpolation has reached the end of the trajectory and the interpolation is stopped; if no, then calculate the segment direction vector X = {X(1), X(2), X(3)} of the current interpolation segment with segment number i; calculate the Euclidean distance inter1 between the segment interpolation start point and the segment end point of the current interpolation segment. Step S203: Determine whether inter1≥intv is true. If true, proceed to step S204; otherwise, proceed to step S207. Intv is the interpolation point spacing. Step S204: In the segment direction of the current interpolation segment, perform equidistant interpolation within the segment successively from the segment interpolation start point until the segment interpolation end point; Step S205: Determine whether the segment interpolation endpoint coincides with the segment endpoint of the current interpolation segment; if yes, determine the segment interpolation starting point of the next segment as the segment starting point; if no, search for a point in the next segment that is at a distance of intv-dist1 from the current segment endpoint as the segment interpolation starting point of the next segment; dist1 is the distance from the segment interpolation endpoint of the current segment to the current segment endpoint; the current segment endpoint is the segment starting point of the next segment. Step S206: After determining the starting point of the next segment interpolation, set i = i + 1 and return to step S202; Step S207: Do not interpolate the current segment; search for a point in the next segment that is at a distance of intv-dist2 from the end point of the current segment as the starting point for segment interpolation of the next segment; dist2 is the length of the current segment; the end point of the current segment is the starting point of the next segment. Step S208: After determining the starting point of the next segment interpolation, set i = i + 1 and return to step S202; Interpolated data is stored in a storage array of length N. The first position of the storage array stores the starting point of the input navigation trajectory. Starting from the second position, the interpolated points generated in step S204 are stored sequentially until the storage array is full, thus forming a source point cloud.

6. The automatic registration method for neurosurgical facial navigation point clouds according to claim 1, characterized in that, Step S3 includes: Step S301: Perform coarse registration; perform principal component analysis on the target point cloud and the source point cloud respectively to obtain the coordinate systems of the target point cloud and the source point cloud. Using the coordinate system of the target point cloud as the reference, adjust the pointing direction of each axis of the source point cloud coordinate system in turn to align it with the coordinate system of the target point cloud; calculate each alignment case and output the rotation and translation homogeneous matrix of the N cases with the smallest Euclidean distance between the origin of the source point cloud and the target point cloud during alignment. Step S302: Perform fine registration; use the N homogeneous rotation and translation matrices from the coarse registration as the initial transformation matrices of the source point cloud, and perform M types of registration between the source point cloud and the target point cloud, including point-to-point, point-to-surface, and surface-to-surface relationships, as well as pairwise combinations of the above three relationships, to obtain M*N rigid body rotation matrices and Euclidean distances. Select the rigid body rotation matrix with the best robustness as the fine registration result output.

7. The automatic registration method for neurosurgical facial navigation point clouds according to claim 6, characterized in that, The coarse registration process in step S301 includes: Step S301-1: Perform principal component analysis on the target point cloud and the source point cloud respectively to obtain the coordinate systems of the target point cloud and the source point cloud; Step S301-2: Using the coordinate system of the target point cloud as a reference, adjust the pointing directions of the +X, -X, +Y, -Y, +Z and -Z axes of the source point cloud coordinate system in sequence, and form various adjustment situations through permutation and combination; Step S301-3: Calculate the rotation and translation homogeneous matrices corresponding to various adjustment cases using the three-point registration method; and eliminate cases where the coordinate systems cannot be aligned after transformation, and select the adjustment cases that can be aligned. Step S301-4: Input the source point cloud and the target point cloud, select the rotation and translation homogeneous matrices corresponding to the alignment adjustment cases as the initial transformation, use KD-tree to find the corresponding points of the source point cloud and the target point cloud, calculate the rigid body transformation matrix by SVD decomposition, and update the position of the iterative trajectory point cloud. Step S301-5: After satisfying the iteration limit, record the Euclidean distance value (fitness_score) between the origin of the source point cloud and the target point cloud in various alignment adjustment situations. Step S301-6: Sort the Euclidean distance values ​​(fitness_score) between the origin and the source point cloud and the target point cloud, and output the rotation and translation homogeneous matrices for the N cases with the smallest fitness_score values.

8. The automatic registration method for neurosurgical facial navigation point clouds according to claim 7, characterized in that, Methods for determining principal component components include: 1) Point cloud data centralization; First, calculate the average value of the point cloud, and then subtract the average value from each point in the point cloud; 2) Construct the covariance matrix; use the centered dataset to calculate the covariance matrix C; 3) Calculate the eigenvalues ​​and eigenvectors; solve for the eigenvalues ​​and corresponding eigenvectors of the covariance matrix C to obtain the principal component vectors; 4) Sort the feature values; arrange the feature vectors in descending order of feature value; the direction corresponding to the largest feature value is the direction of the greatest data change; 5) Determine the coordinate system; select the eigenvectors corresponding to the first 3 largest eigenvalues ​​as the 3 principal component vectors of the point cloud.

9. The automatic registration method for neurosurgical facial navigation point clouds according to claim 8, characterized in that, The fine registration process in step S302 includes: Step S302-1: Using the rotation and translation homogeneous matrices of the N cases output in step S301, perform initial transformations on the source point cloud to obtain N initial source point clouds. Step S302-2: Using N initial source point clouds, perform a total of 7 registrations between the source point cloud and the target point cloud, including point-to-point, point-to-surface, surface-to-surface, point-to-point followed by surface-to-surface, point-to-surface followed by surface-to-surface, surface-to-surface followed by point-to-point, and surface-to-surface followed by point-to-surface, to obtain 7*N rigid body rotation matrices and Euclidean distances. Step S302-3: Select the rigid body rotation matrix with the best robustness as the fine registration result output.

10. The automatic registration method for neurosurgical facial navigation point clouds according to claim 9, characterized in that, In step S302-2, during the registration process between the source point cloud and the target point cloud, first between points and then between surfaces, First, point-to-point registration is performed. Using the obtained rigid body rotation matrix R and translation matrix t, the initial source point cloud is transformed to obtain point cloud data, which is then used as the initial point cloud for subsequent face-to-face registration. This process yields the registered rigid body rotation matrix and the final Euclidean distance value, fitness_score. In step S302-2, during the registration process between the source point cloud and the target point cloud, first between points and surfaces, then between surfaces, First, point-to-surface registration is performed. The obtained rigid body rotation matrix R and translation matrix t are used to transform the initial source point cloud to obtain point cloud data, which is then used as the initial point cloud for subsequent face-to-surface registration. This process yields the registered rigid body rotation matrix and the final Euclidean distance value, fitness_score. In step S302-2, during the registration process between the source point cloud and the target point cloud, first between surfaces and then between points,... First, perform face-to-face registration. Then, use the obtained rigid body rotation matrix R and translation matrix t to transform the initial source point cloud to obtain point cloud data, which is used as the initial point cloud for subsequent point-to-point registration. This process yields the registered rigid body rotation matrix and the final Euclidean distance value, fitness_score. In step S302-2, during the registration process between the source point cloud and the target point cloud, first between surfaces and then between points and surfaces, First, perform face-to-face registration. Then, use the obtained rigid body rotation matrix R and translation matrix t to transform the initial source point cloud to obtain point cloud data, which is used as the initial point cloud for subsequent point-to-face registration. This process yields the registered rigid body rotation matrix and the final Euclidean distance value, fitness_score.

Images