Point cloud registration method applied to neurosurgical navigation

Abstract

The invention provides a point cloud registration method applied to neurosurgical operation navigation, which comprises the following steps: firstly, carrying out initial registration on a preoperative point cloud and an intraoperative point cloud, and converting the preoperative point cloud into a reference point cloud; and then optimizing a rotation matrix and a translation vector required by optimal registration of the reference point cloud and the point cloud in the operation based on an objective function and an optimization algorithm, and after an optimal registration rotation matrix and an optimal registration translation vector are obtained, performing pose transformation on the reference point cloud, so that registration of the point cloud before the operation and the point cloud in the operation is completed. According to the scheme, an invasive physical mark does not need to be applied to a patient for registration; through a mode of initial registration and optimization registration, registration errors caused by weak anti-noise capability and easy falling into local optimum of an algorithm can be greatly eliminated, and the registration precision of the model is improved; meanwhile, the dependence on data in the whole registration process is greatly reduced, and the robustness is higher.

Description

Technical Field

[0001] This invention generally relates to the field of point cloud registration technology, and specifically to a point cloud registration method for neurosurgical navigation. Background Technology

[0002] Neurosurgical navigation systems are crucial for complex neurosurgical procedures such as skull base tumor resection and epileptic focus localization. By establishing a spatial mapping between high-resolution preoperative medical images and real-time intraoperative anatomical structures, they provide surgeons with three-dimensional anatomical localization guidance, directly impacting the accuracy of lesion resection, the protection of key neurovascular structures, and surgical safety. Spatial registration, a core component of the navigation system, requires sub-millimeter-level coordinate system alignment, placing extremely high demands on the robustness and anti-interference capabilities of the algorithm.

[0003] Currently, spatial registration methods in neurosurgery are mainly divided into two categories: one is labeled registration that relies on physical markers, and the other is labelless surface registration based on head surface point cloud data.

[0004] While marker-based registration is relatively mature in clinical applications, it has inherent drawbacks: bone-implanted markers (such as skull screws) require preoperative invasive procedures (drilling for implantation), increasing the risk of infection and patient trauma; manually affixed markers (such as optical reflectors and electromagnetic sensors) are prone to positional drift due to skin deformation and patient head movement, and marker recognition relies on manual matching, which is not only time-consuming but also prone to human error; headframe integrated marker solutions require intraoperative fixation of the stereotactic headframe, which reduces patient comfort and may cause headframe installation deviations to be directly transmitted to the registration results, leading to decreased accuracy and stability. These solutions generally face common problems such as invasiveness, reliance on manual operation, and susceptibility to accuracy drift, which contradicts the demands of minimally invasive and efficient modern neurosurgical procedures.

[0005] While label-free surface registration avoids the intrusiveness of physical labels, it is constrained by both algorithm and data quality. Traditional Iterative Closest Point (ICP) algorithms, which optimize by minimizing the Euclidean distance between corresponding points, are extremely sensitive to intraoperative noise caused by device measurement errors, hair occlusion, and outliers caused by device artifacts. They are also prone to getting trapped in local optima, leading to decreased registration accuracy. Some ICP improvements incorporating kernel correlation entropy still fail to balance noise suppression and local detail matching due to the limited adaptability of single kernel functions and insufficient optimization of multi-kernel combination strategies. Furthermore, the inherent precision of intraoperative point cloud acquisition equipment is insufficient, and thermal drift from prolonged operation further introduces data bias, amplifying algorithmic defects. While recent deep learning registration methods (such as PointNetLK and DCP) can automatically extract point cloud features, they heavily rely on large-scale labeled datasets. In the neurosurgery field, due to significant individual anatomical differences among patients and ethical and privacy restrictions on intraoperative data acquisition, the models' generalization ability is severely limited, making direct application in clinical scenarios difficult.

[0006] In summary, existing registration methods suffer from core bottlenecks: labeled methods are highly invasive, heavily reliant on manual operation, and exhibit poor accuracy stability; labeled methods are susceptible to algorithmic noise resistance, local optima traps, and equipment errors, making it difficult to meet the sub-millimeter registration requirements of neurosurgery. Furthermore, manual coarse-grained registration is time-consuming, and deep learning methods are heavily data-dependent. Therefore, developing a neurosurgical registration method that is non-invasive, highly robust, and highly accurate is a key direction for overcoming existing technological barriers. Summary of the Invention

[0007] In view of the above-mentioned defects or deficiencies in the prior art, it is desirable to provide a point cloud registration method for neurosurgical navigation.

[0008] This invention provides a point cloud registration method for neurosurgical navigation, comprising: S1: Acquire preoperative point cloud and intraoperative point cloud; the preoperative point cloud is a point cloud taken before neurosurgery; the intraoperative point cloud is a point cloud taken during neurosurgery. S2: Based on the preoperative point cloud and the intraoperative point cloud, an initial registration rotation matrix and an initial registration translation vector are obtained; the initial registration rotation matrix and the initial registration translation vector are used to perform pose transformation on the point cloud; S3: Use the initial registration rotation matrix and the initial registration translation vector to perform pose transformation on the preoperative point cloud to obtain a reference point cloud; S4: Based on the objective function and optimization algorithm, the registration rotation matrix and registration translation vector required to transform the pose of the reference point cloud to the intraoperative point cloud are optimized to obtain the optimized registration rotation matrix and optimized registration translation vector; the objective function is used to determine the closeness between the point cloud transformed by the registration rotation matrix and registration translation vector in the iterative process and the intraoperative point cloud; S5: Use the optimized registration rotation matrix and optimized registration translation vector to perform pose transformation on the reference point cloud to complete the registration.

[0009] According to the technical solution provided by the present invention, the acquisition of preoperative point cloud includes: Obtain preoperative CT images; Threshold segmentation is performed on the preoperative CT images to obtain preoperative skin surface images; Three-dimensional reconstruction of preoperative skin surface images was performed to obtain preoperative three-dimensional skin data; Feature points extracted from the preoperative three-dimensional skin data are used to obtain the preoperative point cloud.

[0010] According to the technical solution provided by the present invention, obtaining the initial registration rotation matrix and the initial registration translation vector based on the preoperative point cloud and the intraoperative point cloud includes: Obtain preoperative and intraoperative images; Feature recognition was performed on preoperative images to obtain multiple preoperative two-dimensional feature points; Feature recognition was performed on intraoperative images to obtain multiple intraoperative two-dimensional feature points; Based on the camera intrinsic parameter matrix, the multiple preoperative two-dimensional feature points are transformed into preoperative three-dimensional feature points; Based on the camera intrinsic parameter matrix, the multiple intraoperative two-dimensional feature points are transformed into intraoperative three-dimensional feature points; Based on the optimal rotation matrix algorithm, preoperative 3D feature points, and intraoperative 3D feature points, the initial registration rotation matrix and the initial registration translation vector are obtained.

[0011] According to the technical solution provided by the present invention, the objective function has kernel parameters; the kernel parameters include: the covariance matrix of the anisotropic Gaussian kernel, the weight of the kernel of the cluster, and the center of the kernel of the cluster; The optimization algorithms include: asymptotic nonconvex algorithm, exponential trigonometric optimization algorithm, and quasi-Newton algorithm; The asymptotic nonconvex algorithm is used to update the kernel parameters in the objective function; The exponential triangulation optimization algorithm is used to perform a global search on the current registration rotation matrix and the current registration translation vector, so that the preoperative point cloud transformed by the current registration rotation matrix and the current registration translation vector is closer to the intraoperative point cloud. The current registration rotation matrix is ​​initially the initial registration rotation matrix, and the subsequent ones are the iterative registration rotation matrices obtained from the previous iteration; the current registration translation vector is initially the initial registration translation vector, and the subsequent ones are the iterative registration translation vectors obtained from the previous iteration. The quasi-Newton algorithm is used to perform local refinement based on the initial iterative registration rotation matrix and the initial iterative registration translation vector, so that the preoperative point cloud transformed by the current registration rotation matrix and the current registration translation vector is closer to the intraoperative point cloud.

[0012] According to the technical solution provided by the present invention, the optimization of the registration rotation matrix and registration translation vector required to transform the pose of the reference point cloud to the intraoperative point cloud based on the objective function and optimization algorithm, to obtain the optimized registration rotation matrix and optimized registration translation vector, includes: S4-1: Update the kernel parameters in the objective function using the asymptotic nonconvex algorithm based on the reference point cloud; S4-2: Use the exponential trigonometric optimization algorithm to perform a global search on the current registration rotation matrix and the current registration translation vector to obtain the preliminary iterative registration rotation matrix and the preliminary iterative registration translation vector; S4-3: Based on the preliminary iterative registration rotation matrix and the preliminary iterative registration translation vector, a quasi-Newton algorithm is used for local refinement to obtain the iterative registration rotation matrix and the iterative registration translation vector obtained in this iteration; S4-4: The objective function value is calculated based on the updated kernel parameters, the iterative registration rotation matrix, the iterative registration translation vector, and the objective function; S4-5: If the termination condition is met, the iterative registration rotation matrix corresponding to the historical maximum objective function value is used as the optimized registration rotation matrix, and the iterative registration translation vector corresponding to the historical maximum objective function value is used as the optimized registration translation vector; otherwise, repeat steps S4-1 to S4-5.

[0013] According to the technical solution provided by the present invention, the step of calculating the objective function value based on the updated kernel parameters, the iterative registration rotation matrix, the iterative registration translation vector, and the objective function includes: S4-4-1: The reference point cloud is pose transformed using the iterative registration rotation matrix and the iterative registration translation vector to obtain the iterative point cloud; S4-4-2: Based on the iterative point cloud and the intraoperative point cloud, calculate the residual between each point in the iterative point cloud and each point in the intraoperative point cloud; S4-4-3: Substitute the residual and the updated kernel parameters into the objective function to calculate the objective function value.

[0014] According to the technical solution provided by the present invention, the termination condition is: The norm of the difference between the updated kernel parameters and the kernel parameters before the update is less than a first set threshold. Or the norm of the difference between the objective function value calculated in this iteration and the objective function value calculated in the previous iteration is less than the second set threshold; Or it may reach the maximum number of iterations.

[0015] The beneficial effects of this invention are as follows: To address the problems of large invasiveness, heavy reliance on manual operation, and poor accuracy and stability in existing registration methods, this invention designs a point cloud registration method for neurosurgical surgical navigation. First, initial registration is performed between preoperative and intraoperative point clouds, transforming the preoperative point cloud into a reference point cloud. Then, based on an objective function and optimization algorithm, the rotation matrix and translation vector required for optimized registration of the reference and intraoperative point clouds are optimized. After obtaining the optimized registration rotation matrix and optimized registration translation vector, the pose of the reference point cloud is transformed, completing the registration between the preoperative and intraoperative point clouds. According to this scheme, only preoperative and intraoperative point cloud images need to be captured; subsequent registration can be completed solely through algorithmic calculations. No invasive physical markings are required on the patient during the process. Moreover, the combination of initial and optimized registration significantly reduces errors caused by weak algorithm noise resistance and susceptibility to local optima, significantly improving registration accuracy. Simultaneously, the dependence on data is greatly reduced throughout the registration process, resulting in higher robustness. Attached Figure Description

[0016] Other features, objects, and advantages of the invention will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings: Figure 1 This is a flowchart illustrating a point cloud registration method applied to neurosurgical navigation. Figure 2 Point clouds are used to complete the registration process using a point cloud registration method applied to neurosurgical navigation. Detailed Implementation

[0017] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and not intended to limit it. Furthermore, it should be noted that, for ease of description, only the parts relevant to the invention are shown in the accompanying drawings.

[0018] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other. The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0019] refer to Figure 2 Blue represents intraoperative point cloud data, and red represents preoperative point cloud data after transformation according to the present invention.

[0020] It is evident that the areas captured by the intraoperative point cloud and preoperative point cloud images are the human face.

[0021] In the blue intraoperative point cloud, there are more point clouds on the forehead and fewer on the chin. In the red preoperative point cloud, there are fewer point clouds on the forehead and more on the chin.

[0022] refer to Figure 2 The blue area represents the intraoperative point cloud, and the red area represents the preoperative point cloud after transformation according to the present invention. It can be seen that the intraoperative and preoperative point clouds both depict the human face.

[0023] Among them, the blue intraoperative point cloud has a high point cloud density in the facial area, and a low point cloud density in the outlier points caused by sensor noise and intraoperative artifacts. In addition, the intraoperative point cloud is obtained by a structured light camera, and the point cloud shape is not exactly the same as the preoperative point cloud. The red preoperative point cloud has a relatively uniform point cloud density.

[0024] Furthermore, if the classic method of calculating Euclidean distance is used, the registration process will treat all points as equally important. Due to the uneven distribution of point cloud density, low-density outliers will cause significant interference and make the system overly sensitive to noise. The Euclidean distance of local outliers (such as intraoperative artifacts or noise points) will be amplified, affecting the overall registration accuracy.

[0025] However, actual registration requires more accurate point cloud registration of the actual facial region, eliminating the influence of noise and outliers. Therefore, calculating the Euclidean distance to minimize the sum of Euclidean distances will result in a large registration error. Only by eliminating the influence of noise in the point cloud can a higher registration accuracy be obtained.

[0026] To address the aforementioned technical issues, this invention eliminates the influence of noise in point cloud images on registration by employing an initial registration followed by optimized registration method to calculate the rotation matrix and translation vector required for the registration process. During this process, an objective function and optimization algorithm are used to optimize the rotation matrix and translation vector, ensuring that the final rotation matrix and translation vector reduce the impact of noise points on the registration process, resulting in a more consistent overall approximation between the two point clouds after registration.

[0027] refer to Figure 1 This invention provides a point cloud registration method for neurosurgical navigation, comprising: S1: Acquire preoperative point cloud and intraoperative point cloud; the preoperative point cloud is a point cloud taken before neurosurgery; the intraoperative point cloud is a point cloud taken during neurosurgery. Furthermore, the acquisition of preoperative point clouds includes: Obtain preoperative CT images; Threshold segmentation is performed on the preoperative CT images to obtain preoperative skin surface images; Three-dimensional reconstruction of preoperative skin surface images was performed to obtain preoperative three-dimensional skin data; Feature points extracted from the preoperative 3D skin data are used to obtain preoperative point clouds and preoperative images.

[0028] During the surgery, intraoperative point clouds and intraoperative images were obtained using a depth camera.

[0029] Specifically, threshold segmentation divides image pixels into different categories (such as foreground and background) by setting one or more grayscale thresholds. If a pixel's grayscale value meets the threshold condition (such as being greater than or less than the threshold), it is classified into one category; otherwise, it is classified into another. This allows for the extraction of preoperative skin surface images of the desired imaging area.

[0030] 3D reconstruction refers to the process of recovering the three-dimensional geometric structure and appearance information of an object or scene from two-dimensional images, point clouds, depth maps or other sensor data. Using 3D reconstruction methods, preoperative skin surface images can be restored to preoperative three-dimensional skin data.

[0031] Taking facial recognition as an example, in this embodiment, the process of extracting feature points is as follows: based on the Dlib library, facial recognition is performed on preoperative and intraoperative facial photos, and multiple two-dimensional feature points of the face are obtained (for example, 41 are extracted in this embodiment). Thus, the approximate features of the face can be extracted, and then the initial registration rotation matrix and the initial registration translation vector can be calculated.

[0032] S2: Based on the preoperative point cloud and the intraoperative point cloud, an initial registration rotation matrix and an initial registration translation vector are obtained. The initial registration rotation matrix and the initial registration translation vector are used to perform pose transformation on the point cloud, including: S2-1: Obtain preoperative and intraoperative images; Preoperative images were acquired following the procedure in step S1. Intraoperative images were captured using a depth camera during the surgical procedure.

[0033] S2-2: Perform feature recognition on preoperative images to obtain multiple preoperative two-dimensional feature points; Feature recognition was performed on intraoperative images to obtain multiple intraoperative two-dimensional feature points; In this embodiment, taking the face as an example, the feature points are selected from areas such as the eyes, nose, and lips. Following the same method used with the Dlib library as described above, preoperative two-dimensional feature points and intraoperative two-dimensional feature points are extracted through face recognition.

[0034] S2-3: Converting multiple preoperative two-dimensional feature points into preoperative three-dimensional feature points based on the camera intrinsic parameter matrix, including: Based on the camera intrinsic parameter matrix, the multiple intraoperative two-dimensional feature points are transformed into intraoperative three-dimensional feature points; Based on the principle of pinhole camera inverse projection, the spatial ray direction vector corresponding to the two-dimensional pixel in the camera coordinate system is calculated, and the unit direction vector is obtained after normalizing the vector.

[0035] Then, taking the camera's optical center as the starting point of the ray, candidate points on the ray are retrieved in the 3D point cloud. The distance from each 3D point in the point cloud to the ray is calculated. Candidate points whose distance is less than a set threshold and are located in front of the ray (with positive projection length) are filtered out. The nearest point is selected as the matching 3D point by sorting by distance.

[0036] Ultimately, the matched 3D point becomes the 3D coordinates of the 2D feature point, completing the 2D to 3D conversion.

[0037] S2-4: Based on the optimal rotation matrix algorithm, preoperative 3D feature points, and intraoperative 3D feature points, the initial registration rotation matrix and the initial registration translation vector are obtained.

[0038] Among them, the optimal rotation matrix algorithm is the Kabsch algorithm. The Kabsch algorithm is a method for calculating the optimal rigid body transformation (rotation + translation) between two sets of corresponding points to minimize the root mean square error (RMSE) after alignment. It is often used in fields such as molecular structure alignment, computer vision, and point cloud registration.

[0039] Its core steps include: By removing the centroid, constructing the covariance matrix, performing singular value decomposition, and handling reflections, the rotation matrix and translation vector can be calculated. Applied to this embodiment, the initial registration rotation matrix and initial registration translation vector can be obtained.

[0040] S3: The pose transformation of the preoperative point cloud is performed using the initial registration rotation matrix and the initial registration translation vector to obtain the reference point cloud; the specific calculation process is matrix calculation, which will not be elaborated here.

[0041] In this embodiment, the preoperative point cloud is first transformed to a pose similar to the intraoperative point cloud to avoid an excessively large search space, which would cause the optimization registration algorithm to run too slowly.

[0042] S4: Based on the objective function and optimization algorithm, the registration rotation matrix and registration translation vector required to transform the pose of the reference point cloud to the intraoperative point cloud are optimized to obtain the optimized registration rotation matrix and optimized registration translation vector; the objective function is used to determine the closeness between the point cloud transformed by the registration rotation matrix and registration translation vector in the iterative process and the intraoperative point cloud; The optimization algorithm in this embodiment effectively avoids the problem of traditional registration methods (such as ICP) easily getting trapped in local optima through the synergistic effect of a multi-level mechanism. First, GNC (Asymptotic Non-convex) optimization is adopted. By gradually reducing the control parameter μ, the objective function transitions from a convex problem (easy to converge) to a non-convex problem, thereby dynamically escaping local extrema during the optimization process.

[0043] Secondly, by combining the adaptive population update strategy of the ETO (Exponential Triangle Optimization) algorithm, when the standard deviation of the population fitness is too small, the positions of some particles will be randomly reset to break the local convergence state.

[0044] Secondly, the Exponential Triangulation Optimization (ETO) algorithm is used to perform a global search on the current registration rotation matrix and the current registration translation vector. The ETO algorithm is improved by incorporating an adaptive population update strategy and a memory mechanism. When the fitness standard deviation is too small, the adaptive population update strategy will randomly reset the positions of some particles to avoid getting trapped in local optima. The memory mechanism can avoid some of the randomness of the search and use historical best solutions to guide the search process, so as to reach the optimal solution faster and better.

[0045] Meanwhile, the enhanced multi-kernel correlation entropy (EMKC) models the error distribution through multi-scale kernel functions, avoiding the limitations of a single kernel function. It uses adaptive calculation of point cloud density weights to provide more robust initial parameters for the algorithm.

[0046] Furthermore, the objective function includes kernel parameters; the kernel parameters include the covariance matrix of the anisotropic Gaussian kernel. The weights of the cluster kernels The center of the cluster nucleus .

[0047] The optimization algorithms include: asymptotic nonconvex algorithm, exponential trigonometric optimization algorithm, and quasi-Newton algorithm; The asymptotic nonconvex algorithm is used to update the kernel parameters in the objective function; The exponential triangulation optimization algorithm is used to perform a global search on the current registration rotation matrix and the current registration translation vector, so that the preoperative point cloud transformed by the current registration rotation matrix and the current registration translation vector is closer to the intraoperative point cloud. The current registration rotation matrix is ​​initially the initial registration rotation matrix, and the subsequent ones are the iterative registration rotation matrices obtained from the previous iteration; the current registration translation vector is initially the initial registration translation vector, and the subsequent ones are the iterative registration translation vectors obtained from the previous iteration. The quasi-Newton algorithm is used to perform local refinement based on the initial iterative registration rotation matrix and the initial iterative registration translation vector, so that the preoperative point cloud transformed by the current registration rotation matrix and the current registration translation vector is closer to the intraoperative point cloud.

[0048] Specifically, S4 includes: S4-1: Update the kernel parameters in the objective function using the asymptotic nonconvex algorithm based on the reference point cloud; Specifically, the Generative Nonconvexity (GNC) algorithm is a special strategy for solving highly nonconvex optimization problems. Its algorithmic process includes: constructing a surrogate family and performing asymptotic solutions, ultimately outputting the computational results.

[0049] The core idea is to first replace the original non-convex objective with a convex (or relatively smooth) surrogate function to obtain a global solution; then use this solution as an initial value to gradually restore the non-convexity of the original function until it approximates the true optimal solution. Since each step starts with a high-quality solution from the previous step, the probability of falling into shallow local minima can be significantly reduced.

[0050] Specifically, the steps for updating kernel parameters using an asymptotic nonconvex algorithm include: First, calculate the current transformation parameters. Below, refer to the points in the point cloud. The residuals of the nearest neighbor points in the intraoperative point cloud :

[0051] in, Represents the rotation matrix for iterative registration. This represents the i-th point in the reference point cloud. Represents the iterative registration translation vector. This represents the corresponding nearest neighbor.

[0052] Then, the residual weights are calculated using the GNC algorithm for optimization. :

[0053] in The mean standard deviation of the covariance matrix is ​​used to filter out larger residuals. μ This represents the control parameters for the filtered residual. K-Means clustering yields m clusters, and the mean residual of each cluster is used as the kernel center. The covariance matrix of each cluster is calculated using PCA decomposition. :

[0054] The eigenvector matrix, Let eigenvalues ​​be a diagonal matrix, and take... As anisotropic bandwidth This represents the transpose of a matrix; finally, it involves solving a least-squares problem with regularization terms.

[0055] in, This represents the weight of the kernel of the j-th cluster. The variable that minimizes the objective function. G For Gram matrices, h For the kernel response vector, The regularization coefficient is used to update the weights of the kernel of the j-th cluster. This is to enhance the ability of multi-kernel models to fit multimodal residuals.

[0056] S4-2: Use the Exponential Triangulation Optimization (ETO) algorithm to perform a global search on the current registration rotation matrix and the current registration translation vector to obtain the preliminary iterative registration rotation matrix and the preliminary iterative registration translation vector; Specifically, the exponential trigonometric optimization algorithm is a metaheuristic optimization algorithm. Its core idea is to combine the two very different mathematical characteristics of exponential functions ("rapid decay") and trigonometric functions ("periodic oscillation") to form a search framework of "coarse jumping first, then fine oscillation", thereby achieving a gradual balance between exploration (global) and development (local).

[0057] It comprises an exploration phase and a development phase. The exploration phase, as the early iteration of the algorithm, uses an exponential function to generate large step jumps, rapidly covering the solution space and reducing the risk of getting trapped in local minima. The development phase, the later iteration, uses trigonometric functions to generate small-amplitude periodic perturbations, performing a refined search within the discovered high-quality regions to improve convergence accuracy. The switching threshold and transformation coefficient decrease with the iteration exponent, achieving a "gradual transition."

[0058] In some implementations, the step of performing a global search for the current registration rotation matrix and the current registration translation vector using the Exponential Triangulation Optimization (ETO) algorithm specifically includes: First, define the global objective function:

[0059] Where N is the number of points in the preoperative point cloud. For Centered on Anisotropic Gaussian kernels with bandwidth; then individual values ​​are updated using adaptive rules. The j-th dimension parameter ( Corresponding to Euler angles of rotation, (corresponding translation components), combined with a memory mechanism:

[0060] in, This represents the value of the i-th individual (a search unit in the population) in the j-th dimension during the t-th iteration. This represents the updated value of the i-th individual in the j-th dimension at the (t+1)-th iteration. , As a learning factor, , The numbers are random and all belong to the interval [0,1]. For the individual's historical best, To guide the search direction based on the global historical best; Simultaneously, an adaptive population strategy is employed, which is used when the population fitness variance is too small (e.g., less than 10). -3 To maintain diversity, some individuals are reset during iteration, and a globally approximate optimal solution is obtained through iterative optimization. .

[0061] In this embodiment, a global search can be performed on the registration rotation matrix and the registration translation matrix, enabling the initial iteration process of the optimization algorithm in this embodiment to achieve high accuracy.

[0062] S4-3: Based on the preliminary iterative registration rotation matrix and the preliminary iterative registration translation vector, a quasi-Newton algorithm is used for local refinement to obtain the iterative registration rotation matrix and the iterative registration translation vector obtained in this iteration; In this embodiment, the quasi-Newton algorithm uses the BFGS algorithm. The BFGS algorithm is one of the most commonly used and effective quasi-Newton methods for solving unconstrained nonlinear optimization problems. It avoids calculating and finding the Hessian matrix every time by iteratively constructing an inverse approximation of the Hessian matrix, thus combining the fast convergence of Newton's method with the low computational cost of the gradient method.

[0063] In this embodiment, a more refined local search can be performed on the registration rotation matrix and registration translation matrix, enabling the optimization algorithm in this embodiment to obtain a better rotation matrix and translation vector.

[0064] Specifically, local refinement is performed using a quasi-Newton algorithm, including: by Using the initial values, the BFGS algorithm is used to transform the parameters. Perform local fine-tuning to improve registration accuracy: using the global objective function To optimize the objective, the parameters are iteratively updated using an approximate Hessian matrix and gradient information:

[0065] Where, θ k Let θ be the transformation parameter for the k-th iteration. k+1 Let be the transformation parameters for the (k+1)th iteration. Let be the approximate Hessian matrix for the k-th iteration. This corresponds to the gradient.

[0066] The iteration terminates when the termination condition is met, yielding the optimal transformation parameters. .

[0067] Based on optimal transformation parameters Calculate the optimal registration rotation matrix and translation vector.

[0068] Finally, the initial registration rotation matrix and initial registration translation vector are combined with the optimized registration rotation matrix and optimized registration translation vector to obtain the final transformed rotation matrix and translation vector, thus achieving high-precision registration between the source point cloud and the target point cloud.

[0069] S4-4: The objective function value is calculated based on the updated kernel parameters, the iterative registration rotation matrix, the iterative registration translation vector, and the objective function, including: S4-4-1: The reference point cloud is pose-transformed using the iterative registration rotation matrix and the iterative registration translation vector to obtain the iterative point cloud; this process is also a matrix transformation.

[0070] S4-4-2: Based on the iterative point cloud and the intraoperative point cloud, calculate the residual between each point in the iterative point cloud and each point in the intraoperative point cloud; Specifically, after obtaining the iterative point cloud, the nearest neighbor of each point in the iterative point cloud is found in the intraoperative point cloud. The distance between a point in the iterative point cloud and its corresponding nearest neighbor after optimized registration transformation is defined as the residual. The expression is as follows:

[0071] in, Represents the rotation matrix for iterative registration. Represents points in an iterative point cloud. Represents the iterative registration translation vector. This represents the corresponding nearest neighbor.

[0072] The residuals directly reflect the registration deviation between corresponding points in two point clouds, providing a basic deviation metric for the subsequent calculation of the objective function value.

[0073] For each point in the iterative point cloud, K nearest neighbors are selected to construct a local neighborhood. PCA decomposition is performed on the local neighborhood of each point to obtain the feature value of each point in the three-dimensional dimension. Based on the feature value in the three-dimensional dimension, the point-level anisotropic bandwidth is defined. The specific value of the point-level anisotropic bandwidth is adjusted by pruning operation.

[0074] The expression for point-level anisotropic bandwidth is:

[0075] in, The bandwidth represents the point-level anisotropy, and d represents the dimension. Indicates the scale factor. Represents the feature values ​​in three dimensions. This represents the minimum value of the point-level anisotropic bandwidth. This represents the maximum value of the point-level anisotropic bandwidth. This represents the clipping function. Time to take , Time to take .

[0076] and Based on the point cloud resolution setting, in this embodiment... It is 0.5 times the average spacing of the point cloud. It is twice the average spacing of the point cloud.

[0077] The K-means algorithm was used to cluster all residuals, resulting in multiple clusters. Calculate the kernel-level anisotropic bandwidth for each cluster:

[0078] in, Indicates the core-level anisotropic bandwidth. This indicates the number of points contained in the j-th cluster. It represents the sum of features of all points in the j-th cluster along the d-th dimension.

[0079] After obtaining the kernel-level anisotropic bandwidth for constructing clusters, construct the covariance matrix of the anisotropic Gaussian kernel corresponding to each cluster. :

[0080] in, Construct a diagonal matrix, where the diagonal elements are respectively , , , which is the square of the kernel-level bandwidth, with off-diagonal elements being 0.

[0081] Considering that the residuals of actual point cloud registration usually exhibit a multimodal distribution, a multimodal kernel model is constructed by weighting multiple anisotropic Gaussian kernels on the point cloud. The specific steps are as follows: Define the j-th anisotropic Gaussian kernel function :

[0082] in, π Pi Covariance matrix The determinant, for The inverse matrix of , exp(·) is an exponential function with the natural exponent as the base, where the exponential term is the negative half of the square of the Mahalanobis distance, used to describe the probability distribution of the residuals in the multidimensional space. This is the transpose of the matrix. The center of the kernel of the j-th cluster (take the mean of all residuals within the j-th cluster, i.e.) ).

[0083] Then, a multimodal model is constructed. :

[0084]

[0085] Where m is the total number of anisotropic Gaussian kernels, The weight of the kernel of the j-th cluster is obtained by minimizing the weighted kernel matching loss.

[0086] After constructing the multimodal kernel model of the point cloud, local-global density weights are calculated to enhance the robustness of the final objective function. The steps are as follows: First, calculate the i-th point in the iterative point cloud. Local density :

[0087] Where K0 represents The number of K-nearest neighbors, Indicates the i-th , Indicates the i-th The k-th nearest neighbor, Indicates an indicator function, express The global standard deviation is denoted by s, where s is the density threshold adjustment factor.

[0088] Calculate points in iterative point cloud global density :

[0089] in, The norm of a vector and Indicates the i-th and j-th , meaning any two points in a point cloud.

[0090] The ratio of local density to global density is converted into density weights using a Gaussian function. :

[0091] in, This is the width coefficient of the Gaussian function.

[0092] Finally, integrate the multimodal Gaussian kernel model. With density weight We construct an objective function to optimize registration by maximizing the expected kernel correlation entropy (EMKC). :

[0093] Substitution and Then, we obtained:

[0094] S4-4-3: Substitute the residual and the updated kernel parameters into the objective function to calculate the objective function value.

[0095] At this point, the density weights, which are kernel parameters, Anisotropic Gaussian kernel function and residual All of these are known; by substituting them into the objective function, the value of the objective function can be calculated.

[0096] S4-5: If the termination condition is met, the iterative registration rotation matrix corresponding to the historical maximum objective function value is used as the optimized registration rotation matrix, and the iterative registration translation vector corresponding to the historical maximum objective function value is used as the optimized registration translation vector; otherwise, repeat steps S4-1 to S4-5.

[0097] Specifically, the termination condition is: the norm of the difference between the updated kernel parameters and the kernel parameters before the update is less than a first set threshold. Or the norm of the difference between the objective function value calculated in this iteration and the objective function value calculated in the previous iteration is less than the second set threshold; Or it may reach the maximum number of iterations.

[0098] In this embodiment, the first set threshold is 10. -3 The second threshold is set at 10. -8 The maximum number of iterations is set to 200.

[0099] As mentioned earlier, the control parameter μ in the GNC algorithm decays to less than 10. -3 If the current kernel parameters are stable enough, then it can be considered that the current kernel parameters are stable enough, and further iteration will not bring about significant changes.

[0100] Specifically, the rotation matrix and translation vector corresponding to the historical maximum objective function value are used as the optimized registration rotation matrix and optimized registration translation vector. This ensures that the final optimized registration rotation matrix and optimized registration translation vector of the optimization algorithm can bring the reference point cloud closer to the intraoperative point cloud after the second pose transformation, ultimately achieving a more accurate registration between the preoperative point cloud and the intraoperative point cloud.

[0101] S5: Use the optimized registration rotation matrix and optimized registration translation vector to perform pose transformation on the reference point cloud to complete the registration.

[0102] The above description is merely a preferred embodiment of the present invention and an explanation of the technical principles employed. Those skilled in the art should understand that the scope of the invention is not limited to the specific combination of the above-described technical features, but also includes other technical solutions formed by arbitrary combinations of the above-described technical features or their equivalents without departing from the inventive concept. For example, technical solutions formed by substituting the above-described features with (but not limited to) technical features with similar functions disclosed in this invention.

Claims

1. A point cloud registration method for neurosurgical navigation, characterized in that, include: S1: Acquire preoperative point cloud and intraoperative point cloud; the preoperative point cloud is a point cloud taken before neurosurgery; the intraoperative point cloud is a point cloud taken during neurosurgery. S2: Based on the preoperative point cloud and the intraoperative point cloud, an initial registration rotation matrix and an initial registration translation vector are obtained; the initial registration rotation matrix and the initial registration translation vector are used to perform pose transformation on the point cloud; S3: Use the initial registration rotation matrix and the initial registration translation vector to perform pose transformation on the preoperative point cloud to obtain a reference point cloud; S4: Based on the objective function and optimization algorithm, the registration rotation matrix and registration translation vector required to transform the pose of the reference point cloud to the intraoperative point cloud are optimized to obtain the optimized registration rotation matrix and optimized registration translation vector; the objective function is used to determine the closeness between the point cloud transformed by the registration rotation matrix and registration translation vector in the iterative process and the intraoperative point cloud; S5: Use the optimized registration rotation matrix and optimized registration translation vector to perform pose transformation on the reference point cloud to complete the registration.

2. The point cloud registration method for neurosurgical navigation according to claim 1, characterized in that, The acquisition of preoperative point cloud includes: Obtain preoperative CT images; Threshold segmentation is performed on the preoperative CT images to obtain preoperative skin surface images; Three-dimensional reconstruction of preoperative skin surface images was performed to obtain preoperative three-dimensional skin data; Feature points extracted from the preoperative three-dimensional skin data are used to obtain the preoperative point cloud.

3. The point cloud registration method for neurosurgical navigation according to claim 1, characterized in that, The process of obtaining the initial registration rotation matrix and the initial registration translation vector based on the preoperative point cloud and the intraoperative point cloud includes: Obtain preoperative and intraoperative images; Feature recognition was performed on preoperative images to obtain multiple preoperative two-dimensional feature points; Feature recognition was performed on intraoperative images to obtain multiple intraoperative two-dimensional feature points; Based on the camera intrinsic parameter matrix, the multiple preoperative two-dimensional feature points are transformed into preoperative three-dimensional feature points; Based on the camera intrinsic parameter matrix, the multiple intraoperative two-dimensional feature points are transformed into intraoperative three-dimensional feature points; Based on the optimal rotation matrix algorithm, preoperative 3D feature points, and intraoperative 3D feature points, the initial registration rotation matrix and the initial registration translation vector are obtained.

4. The point cloud registration method for neurosurgical navigation according to claim 1, characterized in that, The objective function contains kernel parameters; Kernel parameters include: the covariance matrix of the anisotropic Gaussian kernel, the weights of the kernels of the clusters, and the center of the kernels of the clusters; The optimization algorithms include: asymptotic nonconvex algorithm, exponential trigonometric optimization algorithm, and quasi-Newton algorithm; The asymptotic nonconvex algorithm is used to update the kernel parameters in the objective function; The exponential triangulation optimization algorithm is used to perform a global search on the current registration rotation matrix and the current registration translation vector, so that the preoperative point cloud transformed by the current registration rotation matrix and the current registration translation vector is closer to the intraoperative point cloud. The current registration rotation matrix is ​​initially the initial registration rotation matrix, and the subsequent ones are the iterative registration rotation matrices obtained from the previous iteration; the current registration translation vector is initially the initial registration translation vector, and the subsequent ones are the iterative registration translation vectors obtained from the previous iteration. The quasi-Newton algorithm is used to perform local refinement based on the initial iterative registration rotation matrix and the initial iterative registration translation vector, so that the preoperative point cloud transformed by the current registration rotation matrix and the current registration translation vector is closer to the intraoperative point cloud.

5. A point cloud registration method for neurosurgical navigation according to claim 4, characterized in that, The process, based on an objective function and optimization algorithm, optimizes the registration rotation matrix and registration translation vector required to transform the pose of the reference point cloud to the intraoperative point cloud, resulting in optimized registration rotation matrix and optimized registration translation vector, including: S4-1: Update the kernel parameters in the objective function using the asymptotic nonconvex algorithm based on the reference point cloud; S4-2: Use the exponential trigonometric optimization algorithm to perform a global search on the current registration rotation matrix and the current registration translation vector to obtain the preliminary iterative registration rotation matrix and the preliminary iterative registration translation vector; S4-3: Based on the preliminary iterative registration rotation matrix and the preliminary iterative registration translation vector, a quasi-Newton algorithm is used for local refinement to obtain the iterative registration rotation matrix and the iterative registration translation vector obtained in this iteration; S4-4: The objective function value is calculated based on the updated kernel parameters, the iterative registration rotation matrix, the iterative registration translation vector, and the objective function; S4-5: If the termination condition is met, the iterative registration rotation matrix corresponding to the historical maximum objective function value is used as the optimized registration rotation matrix, and the iterative registration translation vector corresponding to the historical maximum objective function value is used as the optimized registration translation vector; otherwise, repeat steps S4-1 to S4-5.

6. The point cloud registration method for neurosurgical navigation according to claim 5, characterized in that, The calculation of the objective function value based on the updated kernel parameters, the iterative registration rotation matrix, the iterative registration translation vector, and the objective function includes: S4-4-1: The reference point cloud is pose transformed using the iterative registration rotation matrix and the iterative registration translation vector to obtain the iterative point cloud; S4-4-2: Based on the iterative point cloud and the intraoperative point cloud, calculate the residual between each point in the iterative point cloud and each point in the intraoperative point cloud; S4-4-3: Substitute the residual and the updated kernel parameters into the objective function to calculate the objective function value.

7. A point cloud registration method for neurosurgical navigation according to claim 5, characterized in that, The termination condition is: The norm of the difference between the updated kernel parameters and the kernel parameters before the update is less than a first set threshold. Or the norm of the difference between the objective function value calculated in this iteration and the objective function value calculated in the previous iteration is less than the second set threshold; Or it may reach the maximum number of iterations.

Images