Point cloud registration-based auricle positioning method of otology surgical robot

Abstract

The invention discloses an auricle positioning method of an otology surgical robot based on point cloud registration, and relates to the field of surgical robots. The problem that auricle positioning is difficult to achieve rapidly and accurately in the existing otology operation process is solved. The method comprises the following steps of: 1, performing registration preparation work: planning a target point before an operation, acquiring a three-dimensional image of the auricle of a patient through CT scanning, marking an operation target point in the three-dimensional image of the auricle, and acquiring a three-dimensional point cloud, namely a source point cloud, of the auricle part through a structured light camera; 2, noise reduction and down-sampling are carried out on the point cloud data, and features are extracted; 3, performing coarse registration on the target point cloud and the source point cloud by using an improved RANSAC algorithm; 4, performing fine registration of the target point cloud and the source point cloud by using an LM improved NDT algorithm; and 5, obtaining a rigid body transformation matrix from the source point cloud to the target point cloud, and completing registration positioning. The method is used for point cloud-based registration of the otology surgical robot.

Description

Technical Field

[0001] The present invention relates to the technical field of otology robots, and in particular to an auricle positioning method of an otology surgery robot based on point cloud registration, which is used for positioning the auricle to ensure smooth operation of the otology surgery robot. Background Art

[0002] In recent decades, surgical robots have made great progress and are widely used in multiple clinical departments. Surgical robots have the characteristics of good stability, high precision, and intelligence. They can avoid hand tremors during surgery, reduce postoperative complications, and shorten operation time. Therefore, they have become a hot topic in the field of medical and engineering research. Their huge economic benefits and good performance have attracted more and more attention. Although few ear surgery robots have been approved for clinical use, many breakthroughs have been made.

[0003] For example: The announcement number is: CN116098762A, and the name is an inner ear injection-sampling actuator for an ear surgery robot. It is an actuator installed at the end of an inner ear-sampling surgical robot through the ear canal. It enters the middle ear through a 3mm hole pre-opened in the tympanic membrane through the external auditory canal (which can heal itself after surgery) and bypasses the inner ear structure to move to the round window to deliver the drug to the surgical target location and complete the inner ear tissue fluid biopsy.

[0004] This kind of high-precision surgical robot often needs to be implemented in conjunction with a high-precision auricle positioning method during actual surgery. In the prior art, there are few otology surgery navigation solutions. For example, the few methods for auricle positioning are announced with the number CN117576214A, and the patent name is an invention patent for an ear acupoint positioning method and system based on 3DMM. The ear acupoint positioning method includes multiple steps such as data labeling, registration, and coefficient optimization. First, obtain a real human ear model and a standard human ear model, and annotate the acupoints of the standard human ear model to generate label acupoints; then, through center alignment and size normalization, use ICP point cloud registration to achieve rigid deformation registration; then introduce the 3DMM model for non-rigid deformation registration, calculate the chamfer distance between the real human ear model and the standard human ear model, and use the gradient descent method to optimize the 3DMM coefficients, thereby achieving 3DMM registration; finally, based on the annotated acupoint labels, the nearest neighbor principle is used to map and propagate the acupoints to achieve acupoint positioning on a personalized real auricle model. However, this positioning method requires data annotation before surgery, and the ICP algorithm is time-consuming. It may also fail to converge or the local optimum cannot be iterated to the target point cloud. Even though it is a method for auricle positioning, the auricle positioning accuracy and time are difficult to meet the design requirements. Summary of the invention

[0005] The purpose of the present invention is to solve the problem that it is difficult to quickly and accurately locate the auricle during existing ear surgery, and to provide an auricle positioning method based on point cloud registration for an ear surgery robot.

[0006] The technical solution of the present invention is:

[0007] An auricle positioning method for an ear surgery robot based on point cloud registration comprises the following steps:

[0008] Step 1: Preparation for registration:

[0009] The target points are planned before surgery, and a three-dimensional image of the patient's auricle is obtained through CT scanning, and the surgical target points are marked in the three-dimensional image of the auricle;

[0010] During the operation, a structured light camera is used to obtain a three-dimensional point cloud of the patient's auricle during the operation, namely, the target point cloud;

[0011] Step 2: Preprocess the point cloud data:

[0012] De-noise, downsample and extract features from the source point cloud formed by the target points obtained before the operation and the target point cloud obtained during the operation;

[0013] Step 3: Use the improved RANSAC algorithm to perform rough registration of the target point cloud and the source point cloud;

[0014] Step 4: Use the LM-improved NDT algorithm to perform precise registration of the target point cloud and the source point cloud;

[0015] Step 5: Obtain the rigid body transformation matrix from the source point cloud to the target point cloud to complete the registration and positioning.

[0016] Furthermore, the preprocessing of the point cloud data in step 2 includes the following steps:

[0017] Step 21: Remove interference from other environmental factors:

[0018] Use the straight-through filtering method to filter out points whose coordinate values ​​in the three dimensions of x, y and z are not within the given value range, perform preliminary interference removal on the target point cloud acquired by the structured light camera, and remove the point cloud of the non-auricular part of the target point cloud;

[0019] Step 22: Downsampling:

[0020] Use fast uniform sampling method to downsample the target point cloud and source point cloud at the same time, sampling a point every certain number of points to retain the shape information of the point cloud;

[0021] Step 2 and 3: Use statistical filtering to further reduce the noise of the target point cloud;

[0022] Step 24: Further downsampling the target point cloud and source point cloud:

[0023] The target point cloud and the source point cloud are further downsampled using an improved voxel filtering method, and the target point cloud and the source point cloud are divided into three-dimensional voxel grids, and the point closest to the center of gravity is selected as the representative point;

[0024] Step 25: Use the FPFH algorithm to extract feature points of the target point cloud and the source point cloud.

[0025] Furthermore, step 22 also includes extracting source point cloud information from the structured light camera scanning target point cloud and the CT scanning image, and downsampling the source point cloud and the target point cloud at the same time.

[0026] Furthermore, the process of further denoising the source point cloud in steps 2 and 3 is as follows:

[0027] S231: When a structured light camera collects point clouds, sparse outliers will be generated. Outliers are caused by noise and are characterized by sparse distribution in space. Therefore, if a point cloud at a certain point is less than a certain density, the point cloud is invalid.

[0028] S232: Calculate the average distance from each point to its nearest points, then the distances of all points in the point cloud should form a Gaussian distribution; given the mean and variance, remove points outside the range.

[0029] Furthermore, in step 25, extracting feature points of the target point cloud and the source point cloud using the FPFH algorithm includes the following steps:

[0030] S251: Select the projection coordinate system:

[0031] u=n i

[0032] v=(p j -p i )×u

[0033] ω=u×v

[0034] Among them, u, V and w represent three vectors that establish a rectangular coordinate system axis at the point to be solved;

[0035] S252: Use the following formula to find the relationship between the point to be found and the surrounding points:

[0036] α=v·n j

[0037]

[0038] θ=arc tan(w·n j , u·n j )

[0039] Among them, p i 、p j is a pair of points whose geometric relationship is to be determined, n i 、n j Estimate the normal for two points, with p i Establish the coordinate system o-uvw for the origin, (α, φ, θ) is the feature information of the interaction between the two points, and the result of this step is called the simplified feature histogram SPFH;

[0040] S253: Re-determine the k neighboring points of each point and use the neighboring feature histogram weights to calculate the final histogram SPFH (p k ):

[0041]

[0042] Furthermore, the rough registration of the target point cloud and the source point cloud using the improved RANSAC algorithm in step 3 includes the following steps:

[0043] Step 31: Find n ≥ 3 random sampling points in the source point cloud P and their corresponding points in the target point cloud Q through the nearest neighbor matching of the invariant feature descriptor histogram SPFH;

[0044] Step 32: Calculate the dissimilarity vector δ of the sampling polygon. If the dissimilarity vector δ is less than the set threshold, return to step 31;

[0045] Step 33: Use n sampling points to estimate the hypothesized transformation T;

[0046] Step 34: Apply the hypothesized transformation T to the source point cloud P;

[0047] Step 35: Use the Euclidean distance threshold to find the inner point by the nearest neighbor search between the transformed source point cloud P and the target point cloud Q. If the number of inner points is less than the ratio of the total number of points preset by the algorithm, return to step 31:

[0048]

[0049] Where T represents the estimated coordinate transformation; (Tp-q) 2 Represents the distance from point p of point cloud p to the corresponding point q of point cloud Q;

[0050] Step 36: Re-estimate the rigid transformation matrix using the inlier correspondence

[0051] Step 37: Calculate the distance between corresponding points of the inner point If the distance reaches the current minimum, is the final rigid transformation matrix, otherwise the algorithm will run the specified number of iterations to obtain the final result of coarse registration.

[0052] Furthermore, the calculation method of the difference vector δ in step 32 is as follows:

[0053] S321: Check the ratio between the side lengths of the virtual polygon formed by n sampling points on the object and the scene model; S322: Denote the sampling point corresponding to the FPFH feature as p i ,q i , the side length of the object polygon is:

[0054] d p，i =||p i+1modn -p i ||

[0055] Use the same method to calculate the scene polygon d q，i , and then calculate the difference vector δ by the ratio between the side lengths of n polygons:

[0056]

[0057] Each of the terms added in the above formula represents the ratio of the difference in each dimension to the maximum value of these differences.

[0058] Furthermore, in step 4, the precise registration of the target point cloud and the source point cloud using the LM-improved NDT algorithm includes the following steps:

[0059] Step 41: Divide the point cloud space into several three-dimensional cubes of uniform size, and place the point clouds in the corresponding three-dimensional cubes;

[0060] Step 42: Approximate the point cloud distribution in each cube to obey the normal distribution, and calculate the probability density function PDF parameters in each cube, where PDF is the normal distribution probability density function of each voxel solved in the NDT algorithm, which is used as the matching criterion in NDT registration:

[0061]

[0062] Among them, μ represents the mean vector, ε is the covariance matrix, m is the number of points in the three-dimensional cube, and x k is a point in the cube;

[0063] Step 43: Transform each point of the source point cloud according to the estimated transformation matrix Transform the source point cloud and the target point cloud

[0064] In the same coordinate system;

[0065] Step 44: According to the specific location of the transformed point cloud to be registered in the cube of the target point cloud, calculate the probability distribution function representing the entire point cloud space to be registered:

[0066]

[0067] Step 4 and 5: Use the LM method to calculate the optimal value of all points. The objective function is:

[0068]

[0069] Compared with the prior art, the present invention has the following effects:

[0070] 1. Aiming at the need for positioning the auricle during otolaryngology surgery, the present invention proposes a point cloud-based registration scheme. The present invention uses straight-through filtering and statistical filtering to remove outliers, and uses fast uniform sampling and improved voxel filtering methods to downsample the original data; extracts features based on the Fast Point Feature Histograms (FPFH) algorithm, and uses the improved Random Sample Consensus (RANSAC) algorithm to achieve coarse registration of the point cloud; and uses the Normal Distributions Transform (NDT) algorithm improved based on the Levenberg-Marquardt (LM) method to complete the precise registration of the point cloud. After adopting the above registration method, the target point cloud and the source point cloud data obtained are more accurate.

[0071] 2. The present invention adopts the noise reduction and downsampling method of three-dimensional point cloud to reduce the calculation amount and interference of registration; uses the FPFH feature extraction method to extract the point cloud feature histogram and obtain shape information; successively uses the improved RANSAC algorithm and the LM accelerated NDT algorithm to achieve coarse registration and fine registration, and finally obtains the rigid transformation matrix from the source point cloud to the target point cloud. It is realized that the target points planned in the source point cloud before the operation can be transformed into the target point cloud obtained during the operation through the transformation matrix, and finally realizes the posture coordinates of the preoperative planning points in the intraoperative world coordinate system. The present invention uses feature-based coarse registration to ensure the completion of the registration task, and the LM improved NDT algorithm also improves the robustness; the human ear model of a specific patient is obtained by CT scanning before the operation, and there is no need to perform scaling operations on the standard human ear model, which not only simplifies the registration process, but also improves the registration accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 It is a flowchart of the auricle positioning method of the present invention;

[0073] Figure 2 It is a schematic diagram of the process of preprocessing point cloud data according to the present invention;

[0074] Figure 3 It is a schematic diagram of a process of further downsampling the target point cloud and the source point cloud using the improved voxel filtering method of the present invention;

[0075] Figure 4 It is a schematic diagram of the registration process of the present invention. DETAILED DESCRIPTION

[0076] Specific implementation method 1: Combination Figures 1 to 4 To illustrate this embodiment, the embodiment includes the following steps:

[0077] Step 1: Preparation for registration:

[0078] The target points are planned before surgery, and a three-dimensional image of the patient's auricle is obtained through CT scanning, and the surgical target points are marked in the three-dimensional image of the auricle;

[0079] During the operation, a structured light camera is used to obtain a three-dimensional point cloud of the patient's auricle during the operation, namely, the target point cloud;

[0080] Step 2: Preprocess the point cloud data:

[0081] De-noise, downsample and extract features from the source point cloud formed by the target points obtained before the operation and the target point cloud obtained during the operation;

[0082] Step 3: Use the improved RANSAC algorithm to perform rough registration of the target point cloud and the source point cloud;

[0083] Step 4: Use the LM-improved NDT algorithm to perform precise registration of the target point cloud and the source point cloud;

[0084] Step 5: Obtain the rigid body transformation matrix from the source point cloud to the target point cloud to complete the registration and positioning.

[0085] Specific implementation method 2: Combination Figures 1 to 3 To illustrate this embodiment, the step 2 of this embodiment preprocesses the point cloud data including the following steps:

[0086] Step 21: Remove interference from other environmental factors:

[0087] Use the straight-through filtering method to filter out points whose coordinate values ​​in the three dimensions of x, y and z are not within the given value range, perform preliminary interference removal on the target point cloud acquired by the structured light camera, and remove the point cloud of the non-auricular part of the target point cloud;

[0088] Step 22: Downsampling:

[0089] Use fast uniform sampling method to downsample the target point cloud and source point cloud at the same time, sampling a point every certain number of points to retain the shape information of the point cloud;

[0090] Step 2 and 3: Use statistical filtering to further reduce the noise of the target point cloud;

[0091] Step 24: Further downsampling the target point cloud and source point cloud:

[0092] The target point cloud and the source point cloud are further downsampled using an improved voxel filtering method, and the target point cloud and the source point cloud are divided into three-dimensional voxel grids, and the point closest to the center of gravity is selected as the representative point;

[0093] Step 25: Use the FPFH algorithm to extract feature points of the target point cloud and the source point cloud.

[0094] With this arrangement, the feature histogram is extracted, which contains the 33-dimensional feature information of the points to be extracted, and provides preparation for the next step of feature-based rough registration. The other components and connection relationships are the same as those in the first specific implementation.

[0095] Specific implementation method three: Combination Figures 1 to 3 To illustrate this embodiment, step 22 of this embodiment also includes extracting source point cloud information from the structured light camera scanning target point cloud and the CT scanning image, and downsampling the source point cloud and the target point cloud at the same time.

[0096] In this way, RANSAC performs feature-based matching, and can achieve registration when the poses of two point clouds are greatly different. The improved RANSAC algorithm further simplifies the algorithm and improves the efficiency of the algorithm. Other components and connection relationships are the same as those of the first or second specific implementation.

[0097] Specific implementation method four: Combination Figure 3 To illustrate this implementation, the process of further denoising the source point cloud in steps 2 and 3 of this implementation is as follows:

[0098] S231: When a structured light camera collects point clouds, sparse outliers will be generated. Outliers are caused by noise and are characterized by sparse distribution in space. Therefore, if a point cloud at a certain point is less than a certain density, the point cloud is invalid.

[0099] S232: Calculate the average distance from each point to its nearest points, then the distances of all points in the point cloud should form a Gaussian distribution; given the mean and variance, remove points outside the range.

[0100] With such a configuration, the NDT algorithm is a probability-based registration algorithm that can achieve registration in a relatively short time. The LM acceleration is used to dynamically adjust the step size during the iteration process to avoid non-convergence or local optimality, thereby improving robustness and accuracy. The other components and connection relationships are the same as any one of the specific implementations one to three.

[0101] Specific implementation method five: Combination Figure 3 To illustrate this embodiment, in step 25 of this embodiment, using the FPFH algorithm to extract feature points of the target point cloud and the source point cloud includes the following steps:

[0102] S251: Select the projection coordinate system:

[0103] u=n i

[0104] v=(p j -p i )×u

[0105] ω=u×v

[0106] Among them, u, V and w represent three vectors that establish a rectangular coordinate system axis at the point to be solved;

[0107] S252: Use the following formula to find the relationship between the point to be found and the surrounding points:

[0108] α=v·n j

[0109]

[0110] θ=arctan(w·n j , u·n j )

[0111] Among them, p i 、p j is a pair of points whose geometric relationship is to be determined, n i 、n j Estimate the normal for two points, with p i Establish the coordinate system o-uvw for the origin, (α, φ, θ) is the feature information of the interaction between the two points, and the result of this step is called the simplified feature histogram SPFH;

[0112] S253: Re-determine the k neighboring points of each point and use the neighboring feature histogram weights to calculate the final histogram SPFH (p k ):

[0113]

[0114] The other components and connection relationships are the same as any one of the specific implementation modes 1 to 4.

[0115] Specific implementation method six: Combination Figure 4 To illustrate this embodiment, the step 3 of this embodiment uses the improved RANSAC algorithm to perform rough registration of the target point cloud and the source point cloud, including the following steps:

[0116] Step 31: Find n ≥ 3 random sampling points in the source point cloud P and their corresponding points in the target point cloud Q through the nearest neighbor matching of the invariant feature descriptor histogram SPFH;

[0117] Step 32: Calculate the dissimilarity vector δ of the sampling polygon. If the dissimilarity vector δ is less than the set threshold, return to step 31;

[0118] Step 33: Use n sampling points to estimate the hypothesized transformation T;

[0119] Step 34: Apply the hypothesized transformation T to the source point cloud P;

[0120] Step 35: Use the Euclidean distance threshold to find the inner point by the nearest neighbor search between the transformed source point cloud P and the target point cloud Q. If the number of inner points is less than the ratio of the total number of points preset by the algorithm, return to step 31:

[0121]

[0122] Where T represents the estimated coordinate transformation; (Tp-q) 2 Represents the distance from point p of point cloud p to the corresponding point q of point cloud Q;

[0123] Step 36: Re-estimate the rigid transformation matrix using the inlier correspondence

[0124] Step 37: Calculate the distance between corresponding points of the inner point If the distance reaches the current minimum, is the final rigid transformation matrix, otherwise the algorithm will run the specified number of iterations to obtain the final result of coarse registration.

[0125] The other components and connection relationships are the same as any one of the specific implementation modes 1 to 5.

[0126] Specific implementation method seven: Combination Figure 4 To illustrate this embodiment, the calculation method of the difference vector δ in step 32 of this embodiment is as follows:

[0127] S321: Check the ratio between the side lengths of the virtual polygon formed by n sampling points on the object and the scene model; S322: Denote the sampling point corresponding to the FPFH feature as p i ,q i , the side length of the object polygon is:

[0128] d p，i =||p i+1modn -p i ||

[0129] Use the same method to calculate the scene polygon d q，i , and then calculate the difference vector δ by the ratio between the side lengths of n polygons:

[0130]

[0131] Each of the terms added in the above formula represents the ratio of the difference in each dimension to the maximum value of these differences.

[0132] The other components and connection relationships are the same as any one of the specific implementation modes 1 to 6.

[0133] Specific implementation method eight: Combination Figure 4 To illustrate this embodiment, in step 4 of this embodiment, using the LM-improved NDT algorithm to perform precise registration of the target point cloud and the source point cloud includes the following steps:

[0134] Step 41: Divide the point cloud space into several three-dimensional cubes of uniform size, and place the point clouds in the corresponding three-dimensional cubes;

[0135] Step 42: Approximate the point cloud distribution in each cube to obey the normal distribution, and calculate the probability density function PDF parameters in each cube, where PDF is the normal distribution probability density function of each voxel solved in the NDT algorithm, which is used as the matching criterion in NDT registration:

[0136]

[0137] Among them, μ represents the mean vector, ε is the covariance matrix, m is the number of points in the three-dimensional cube, and x k is a point in the cube;

[0138] Step 43: Transform each point of the source point cloud according to the estimated transformation matrix Transform the source point cloud and the target point cloud

[0139] In the same coordinate system;

[0140] Step 44: According to the specific location of the transformed point cloud to be registered in the cube of the target point cloud, calculate the probability distribution function representing the entire point cloud space to be registered:

[0141]

[0142] Step 4 and 5: Use the LM method to calculate the optimal value of all points. The objective function is:

[0143]

[0144] The other components and connection relationships are the same as any one of the specific implementation modes 1 to 7.

[0145] Combination Figures 1 to 4 Describe the embodiments of the present invention:

[0146] Step 1: Plan the target points before surgery and obtain a 3D image of the patient's auricle through CT scanning. Mark the target points in the image and obtain a 3D point cloud of the patient's auricle during surgery through a structured light camera.

[0147] Step 2: Preprocess the source point cloud obtained before surgery and the target point cloud obtained during surgery;

[0148] Step 21: Use the straight-through filtering method to filter out points whose coordinate values ​​in the three dimensions of x, y and z are not within the given value range, and filter the target point cloud obtained by the structured light camera (such as Figure 2 (a) Perform preliminary interference removal to remove most of the non-auricular part of the source point cloud, such as Figure 2 (b)

[0149] Step 22: Use the fast uniform sampling method to downsample. It is suitable for scenes with a large number of point clouds. Sampling a point every certain number of points can retain the shape information of the point cloud. Extract point cloud information from the structured light camera scan point cloud and CT scan image ( Figure 2 (c)), downsampling the source point cloud and the target point cloud simultaneously, such as Figure 2 As shown in (d) and (e), the amount of subsequent calculations is reduced and the efficiency is improved;

[0150] Step 2 and 3: Use statistical filtering to further reduce the noise of the target point cloud, such as Figure 3 (a). When a structured light camera collects point clouds, sparse outliers will be generated. Outliers are often caused by noise and are characterized by being sparsely distributed in space. It can be defined that if a point cloud is less than a certain density, the point cloud is invalid. Calculate the average distance from each point to its nearest points. The distances of all points in the point cloud should form a Gaussian distribution. Given the mean and variance, points outside the range can be removed;

[0151] Step 24: Use the improved voxel filtering method to further downsample the two point clouds and divide the point clouds into three-dimensional voxel grids, such as Figure 3 (b)(d), select the point closest to the center of gravity as the representative point, such as Figure 3 (c)(e), improve the accuracy of point cloud representation after downsampling.

[0152] Step 25: Use FPFH algorithm to extract two point cloud feature points:

[0153] Select a projected coordinate system:

[0154] u=n i

[0155] v=(p j -p i )×u

[0156] ω=u×v

[0157] Use the following formula to find the relationship between the point to be found and the surrounding points:

[0158] α=v·n j

[0159]

[0160] θ=arctan(w·n j , u·n j )

[0161] Among them, p i 、p j is a pair of points whose geometric relationship is to be determined, n i 、n j Estimate the normal for two points, with p i Establish the coordinate system o-uvw for the origin, (α, φ, θ) is the feature information of the interaction between the two points, and the result of this step is called the simplified feature histogram SPFH;

[0162] Step 26: Next, re-determine the k neighboring points of each point, and use the neighboring SPFH and weight them to calculate the final histogram:

[0163]

[0164] Step 3: Use the improved RANSAC algorithm to perform rough registration of the two point clouds;

[0165] Step 31: Find n ≥ 3 random sampling points in the source point cloud P and their corresponding points in the target point cloud Q through the nearest neighbor matching of the invariant feature descriptor FPFH;

[0166] Step 3: Calculate the difference vector δ of the sampling polygon. If it is less than the set threshold, return to step 31;

[0167] in:

[0168] The ratio between the side lengths of the virtual polygon formed by n sampling points on the object and the scene model is checked.

[0169] The sampling point corresponding to the feature is denoted as p i ,q i , the side length of the object polygon is:

[0170] d p，i =||p i+1modn -p i ||

[0171] The same method can be used to calculate the scene polygon dqi, and then the difference vector δ is calculated by the ratio between the edge lengths of n polygons:

[0172]

[0173] Step 33: Use n sampling points to estimate the assumed rigid transformation matrix T;

[0174] Step 34: Apply the assumed rigid transformation matrix to the source point cloud P;

[0175] Step 35: Use the Euclidean distance threshold to find the inner point by the nearest neighbor search between the transformed source point cloud P and the target point cloud Q. If the number of inner points is less than the ratio of the total number of points preset by the algorithm, return to step 31:

[0176]

[0177] Step 36: Re-estimate the hypothesized transformation using the inlier correspondence;

[0178] Step 37: Calculate the distance ε(T) between the corresponding points of the inner point. If the distance reaches the current minimum value, T is the final transformation. Otherwise, the algorithm will run the specified number of iterations. The final result of the rough registration is as follows: Figure 4 (b) as shown.

[0179] Step 4: Use the LM-improved NDT algorithm to perform precise registration of the two point clouds;

[0180] Step 41: Divide the point cloud space into several three-dimensional cubes of uniform size, and place the point clouds in the corresponding three-dimensional cubes;

[0181] Step 42: The point cloud distribution in each cube is approximated to obey the normal distribution, so the PDF (probability density function) parameters in each cube can be calculated:

[0182]

[0183] Where μ represents the mean vector, ε is the covariance matrix, m is the number of points in the three-dimensional cube, and x k is a point in the cube.

[0184] Step 43: Transform each point of the source point cloud according to the estimated transformation matrix Transform them so that they are in the same coordinate system;

[0185] Step 44: According to the specific location of the transformed point cloud to be registered in the cube of the target point cloud, calculate the probability distribution function representing the entire point cloud space to be registered:

[0186]

[0187] Step 4 and 5: Use the Levenberg-Marquardt (LM) method to calculate the optimal value of all points. The objective function is:

[0188]

[0189] The final result of fine alignment is as follows Figure 4 (c) as shown.

[0190] Step 5: Get the rigid body transformation matrix from the source point cloud to the target point cloud.

[0191] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some of the technical features may be replaced by equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. An auricle positioning method for an ear surgery robot based on point cloud registration, characterized in that: It includes the following steps: Step 1: Preparation for registration: The target points are planned before surgery, and a three-dimensional image of the patient's auricle is obtained through CT scanning, and the surgical target points are marked in the three-dimensional image of the auricle; During the operation, a structured light camera is used to obtain a three-dimensional point cloud of the patient's auricle during the operation, namely, the target point cloud; Step 2: Preprocess the point cloud data: De-noise, downsample and extract features from the source point cloud formed by the target points obtained before the operation and the target point cloud obtained during the operation; Step 3: Use the improved RANSAC algorithm to perform rough registration of the target point cloud and the source point cloud; Step 4: Use the LM-improved NDT algorithm to perform precise registration of the target point cloud and the source point cloud; Step 5: Obtain the rigid body transformation matrix from the source point cloud to the target point cloud to complete the registration and positioning.

2. The auricle positioning method based on point cloud registration of an ear surgery robot according to claim 1, characterized in that: The preprocessing of point cloud data in step 2 includes the following steps: Step 21: Remove interference from other environmental factors: Use the straight-through filtering method to filter out points whose coordinate values ​​in the three dimensions of x, y and z are not within the given value range, perform preliminary interference removal on the target point cloud acquired by the structured light camera, and remove the point cloud of the non-auricular part of the target point cloud; Step 22: Downsampling: Use fast uniform sampling method to downsample the target point cloud and source point cloud at the same time, sampling a point every certain number of points to retain the shape information of the point cloud; Step 2 and 3: Use statistical filtering to further reduce the noise of the target point cloud; Step 24: Further downsampling the target point cloud and source point cloud: The target point cloud and the source point cloud are further downsampled using an improved voxel filtering method, and the target point cloud and the source point cloud are divided into three-dimensional voxel grids, and the point closest to the center of gravity is selected as the representative point; Step 25: Use the FPFH algorithm to extract feature points of the target point cloud and the source point cloud.

3. The auricle positioning method based on point cloud registration of an ear surgery robot according to claim 2, characterized in that: Step 22 also includes extracting source point cloud information from the structured light camera scanning target point cloud and the CT scanning image, and downsampling the source point cloud and the target point cloud at the same time.

4. The auricle positioning method based on point cloud registration of an ear surgery robot according to claim 3, characterized in that: The process of further denoising the source point cloud in steps 2 and 3 is as follows: S231: When a structured light camera collects point clouds, sparse outliers will be generated. Outliers are caused by noise and are characterized by sparse distribution in space. Therefore, if a point cloud at a certain point is less than a certain density, the point cloud is invalid. S232: Calculate the average distance from each point to its nearest points, then the distances of all points in the point cloud should form a Gaussian distribution; given the mean and variance, remove points outside the range.

5. The auricle positioning method based on point cloud registration of an ear surgery robot according to claim 4, characterized in that: In step 25, using the FPFH algorithm to extract feature points of the target point cloud and the source point cloud includes the following steps: S251: Select the projection coordinate system: u=n i v=(p j -p i )×u ω=u×v Among them, u, v and w represent three vectors that establish a rectangular coordinate system axis at the point to be solved; S252: Use the following formula to find the relationship between the point to be found and the surrounding points: α=v·n j θ=arctan(w·n j ,u·n j ) Among them, p i 、p j is a pair of points whose geometric relationship is to be determined, n i 、n j Estimate the normal for two points, with p i Establish the coordinate system o-uvw for the origin, (α, φ, θ) is the feature information of the interaction between the two points, and the result of this step is called the simplified feature histogram SPFH; S253: Re-determine the k neighboring points of each point and use the neighboring feature histogram weights to calculate the final histogram SPFH (p k ):

6. The auricle positioning method based on point cloud registration of an ear surgery robot according to claim 1 or 5, characterized in that: In step 3, the rough registration of the target point cloud and the source point cloud using the improved RANSAC algorithm includes the following steps: Step 31: Find n ≥ 3 random sampling points in the source point cloud P and their corresponding points in the target point cloud Q through the nearest neighbor matching of the invariant feature descriptor histogram SPFH; Step 32: Calculate the dissimilarity vector δ of the sampling polygon. If the dissimilarity vector δ is less than the set threshold, return to step 31; Step 33: Use n sampling points to estimate the hypothesized transformation Step 34: Transform the hypothesis Applied to the source point cloud P; Step 35: Use the Euclidean distance threshold to find the inner point by the nearest neighbor search between the transformed source point cloud P and the target point cloud Q. If the number of inner points is less than the ratio of the total number of points preset by the algorithm, return to step 31: Where T represents the estimated coordinate transformation; (Tp-q) 2 Represents the distance from point p of point cloud p to the corresponding point q of point cloud Q; Step 36: Re-estimate the rigid transformation matrix using the inlier correspondence Step 37: Calculate the distance between corresponding points of the inner point If the distance reaches the current minimum, is the final rigid transformation matrix, otherwise the algorithm will run the specified number of iterations to obtain the final result of coarse registration.

7. The auricle positioning method based on point cloud registration of an ear surgery robot according to claim 6, characterized in that: The calculation method of the dissimilarity vector δ in step 3.2 is as follows: S321: Checking the ratio between the side lengths of the virtual polygon formed by n sampling points on the object and the scene model; S322: The sampling point corresponding to the FPFH feature is represented as p i ,q i , the side length of the object polygon is: d p，i =||p i+1modn -p i || Use the same method to calculate the scene polygon d q,i , and then calculate the difference vector δ by the ratio between the side lengths of n polygons: Each of the terms added in the above formula represents the ratio of the difference in each dimension to the maximum value of these differences.

8. The auricle positioning method based on point cloud registration of an ear surgery robot according to claim 1 or 7, characterized in that: In step 4, the precise registration of the target point cloud and the source point cloud using the LM-improved NDT algorithm includes the following steps: Step 41: Divide the point cloud space into several three-dimensional cubes of uniform size, and place the point clouds in the corresponding three-dimensional cubes; Step 42: Approximate the point cloud distribution in each cube to obey the normal distribution, and calculate the probability density function PDF parameters in each cube, where PDF is the normal distribution probability density function of each voxel solved in the NDT algorithm, which is used as the matching criterion in NDT registration: Among them, μ represents the mean vector, ε is the covariance matrix, m is the number of points in the three-dimensional cube, and x k is a point in the cube; Step 43: Transform each point of the source point cloud according to the estimated transformation matrix Transform so that the source point cloud and the target point cloud are in the same coordinate system; Step 44: According to the specific location of the transformed point cloud to be registered in the cube of the target point cloud, calculate the probability distribution function representing the entire point cloud space to be registered: Step 4 and 5: Use the LM method to calculate the optimal value of all points. The objective function is:

Images