Global registration method based on spherical harmonic transform and iterative optimization

Abstract

The technical solution of this invention provides a global registration method based on spherical harmonic transformation and iterative optimization. Its key features include the following assumptions: First, a predefined minimum percentage of correctly matched pairs in the input point cloud, represented as a limit on the number of outliers in the point cloud; second, a distance threshold for determining correctly matched pairs is preset by the scene, represented as a limit distance for outliers in the point cloud. In the algorithm provided by this invention, the point cloud first undergoes coarse registration to obtain several search domains. Then, the search domains are sorted using evaluation criteria, and incorrect search domains are eliminated using a branch and bound method to obtain the final registration result. The algorithm disclosed in this invention can be applied to 3D reconstruction. Because the algorithm of this invention is faster than all current global registration algorithms at the same accuracy, it can achieve higher accuracy than all current algorithms in the same time frame for 3D reconstruction.

Description

Technical Field

[0001] This invention relates to a global registration method. Background Technology

[0002] Global registration involves finding the optimal coordinate system transformation to make two overlapping point cloud datasets in different coordinate systems coincide. Global registration of point clouds is widely used in fields lacking initial localization or requiring high algorithm stability, such as motion estimation, 3D reconstruction, 3D object recognition, 3D object detection, autonomous driving, and intelligent robotics.

[0003] Due to the non-convexity of the registration optimization problem, traditional algorithms are prone to getting trapped in local minima. Global registration methods, through theoretical proof, guarantee that the algorithm's result will converge to the global optimum. Existing global registration algorithms mainly consist of two parts: those based on feature point correspondence and those based on raw point cloud data. Existing global registration implementations either require a long time to narrow down the search range or cannot achieve sufficient accuracy. Therefore, they are difficult to apply in practical scenarios where faster speed and higher accuracy are required. A good global registration algorithm should meet the following four conditions:

[0004] 1) The algorithm should focus on the original point cloud rather than feature points to maintain high accuracy;

[0005] 2) The algorithm complexity should be lower than existing algorithms;

[0006] 3) The algorithm should ensure that the search domain is rapidly reduced to the global optimum during the initial computation;

[0007] 4) When the algorithm is used in applications requiring high precision, the convergence speed should be accelerated, which means reducing the number of subsequent iterations and removing redundant search domains. Summary of the Invention

[0008] The technical problem that this invention aims to solve is that existing global registration implementations either require a long time to narrow down the search range or cannot achieve sufficient accuracy.

[0009] To address the aforementioned technical problems, the present invention provides a global registration method based on spherical harmonic transformation and iterative optimization. The method is characterized by the following assumptions: First, a minimum percentage of correctly matched pairs in the input point cloud is predefined, representing a limit on the number of outliers in the point cloud; second, a distance threshold for determining correctly matched pairs is preset by the scene, representing a limit distance for outliers in the point cloud. The global registration extraction method then includes the following steps:

[0010] Step 1: Align the point cloud in the frequency domain using spherical harmonics based on the second point assumption, and determine several small point cloud search domains based on the peak value of the phase difference and the first point assumption.

[0011] Step 2: Assign a search order to the point cloud search domains determined in Step 1 by analyzing the hit rate and search time of the search domains, where:

[0012] The formula for estimating hit rate is shown below:

[0013]

[0014] In the formula: LB represents the lower bound, N in Represents the number of point cloud matching pairs, residual i Let represent the matching error of the i-th matching pair, radius represent the radius of the search domain, C is a constant, caf is the centering decay factor; r represents the radius of the search domain, f(·) represents the selectability function; E * This indicates the size of the current optimal optimization function;

[0015] The formula for estimating search time is as follows:

[0016] F t (r)=r o *(LB-E * )

[0017] In the formula: o is an optional parameter;

[0018] Step 3: Using the branch and bound method, search sequentially in several point cloud search domains according to the order determined in Step 2. During the search, if the optimal solution is not within the current search domain, the branch and bound method will exclude the current search domain. If the optimal solution is within the current search domain, the branch and bound method will find the optimal solution. The following check condition is introduced during the branch and bound process:

[0019] N in ≤H T *N

[0020] In the formula, H T This refers to the proportion of matching pairs in the point cloud in hypothesis 1, where N represents the size of the input point cloud.

[0021] Preferably, step 1 includes the following steps:

[0022] Step 101: Voxelize the two input point clouds and perform a spherical harmonic transform to obtain the frequency domain signals of the two input point clouds, where:

[0023] The voxelization formula used when voxelizing two input point clouds is shown below:

[0024]

[0025] In the formula: F([x, y, z]) T The obtained voxelized point cloud can be considered as a function; PC represents the input point cloud; max represents the maximum norm of the points in the input point cloud.

[0026] The spherical harmonic transformation is calculated using the following formula:

[0027]

[0028]

[0029]

[0030] In the formula: (x, y, z) represent the coordinates of a point, m represents the superscript of the spherical harmonic function, and l represents the subscript of the spherical harmonic function. Represents a spherical harmonic function. Let F denote the adjoint Legendre polynomial, d denote the differential, and F is F([x, y, z]). T ), Represents the basis coefficients of the spherical harmonic function;

[0031] The result obtained by calculating the spherical harmonic transformation is:

[0032] Step 102: Calculate the standard correlation function SC of the frequency domain signals of the two input point clouds. R (g), as shown in the following formula:

[0033] sC R (g)=∫ x,y,z MSA(g)TAdxdydz

[0034]

[0035]

[0036] In the formula: MSA(g) is... Let (α, β, γ) represent the spherical harmonic basis coefficient matrix of the model point cloud, (α, β, γ) represent the Euler angles corresponding to the rotation, and km and l represent subscripts and superscripts, which will be eliminated during summation. Denotes the Wigner function, d k m l Let represent the Wigner distribution function.

[0037] Step 103: Find all local maxima of the standard correlation function, and calculate the radius of the search domain based on the error function value of each local maximum. For a given local maximum, the radius of the corresponding search domain is... Where E is the size of the optimization function corresponding to the local maximum, and D is the restricted distance of the point cloud outlier in Hypothesis 2.

[0038] Preferably, the branch and bound method includes the following steps: For each search domain, calculate the error function value of the center of the current search domain, and estimate the lower bound of the function of the current search domain through the error function value: If the lower bound of the function is greater than the existing best result, then exclude the search domain; otherwise, recursively split the current search domain into smaller search domains and perform this process. When the size of the search domain is smaller than the target accuracy requirement, the center point of the search domain is the optimal solution.

[0039] In the algorithm provided by this invention, the point cloud is first coarsely registered to obtain several search domains. Then, the search domains are sorted using evaluation criteria, and incorrect search domains are eliminated using the branch and bound method to obtain the final registration result. The algorithm disclosed in this invention can be applied to 3D reconstruction. Since the algorithm of this invention is faster than all current global registration algorithms at the same accuracy, it can achieve higher accuracy than all current algorithms in the same amount of time for 3D reconstruction. Attached Figure Description

[0040] Figure 1 The global registration process of the present invention is illustrated. Detailed Implementation

[0041] The present invention will be further illustrated below with reference to specific embodiments. It should be understood that these embodiments are for illustrative purposes only and are not intended to limit the scope of the invention. Furthermore, it should be understood that after reading the teachings of this invention, those skilled in the art can make various alterations or modifications to the invention, and these equivalent forms also fall within the scope defined by the appended claims.

[0042] This invention analyzes the complexity of existing global registration algorithms and proposes two basic assumptions to accelerate the global registration solution: Assumption 1) a minimum percentage of correct matching pairs in the input point cloud is predefined, which is expressed as a limit on the number of outliers in the point cloud; Assumption 2) the distance threshold used to determine the correct matching pair is preset by the scene, which is expressed as a limit distance for outliers in the point cloud.

[0043] Based on the above two assumptions, this invention proposes a fast global registration algorithm with spherical harmonic transform and iterative optimization strategy, specifically including the following:

[0044] Step 1: This invention proposes a coarse-grained registration method that narrows the entire search domain to several smaller search domains. This invention aligns the point cloud in the frequency domain using spherical harmonics based on Hypothesis 2 (the limiting distance of point cloud outliers), and determines the smaller search domains based on the peak value of the phase difference and Hypothesis 1 (the limiting number of point cloud outliers). Through this method, this invention obtains several point cloud search domains from two input point clouds and two hypotheses, specifically including the following steps:

[0045] Step 101: Voxelize the two input point clouds and perform a spherical harmonic transform to obtain the frequency domain signals of the two input point clouds, where:

[0046] The voxelization formula used when voxelizing two input point clouds is shown below:

[0047]

[0048] In the formula: F([x, y, z]T) is the obtained voxelized point cloud, which can be regarded as a function; PC represents the input point cloud; max represents the maximum norm of the points in the input point cloud.

[0049] The spherical harmonic transformation is calculated using the following formula:

[0050]

[0051]

[0052]

[0053] In the formula: (x, y, z) represent the coordinates of a point, m represents the superscript of the spherical harmonic function, and l represents the subscript of the spherical harmonic function. Represents a spherical harmonic function. Let F denote the adjoint Legendre polynomial, d denote the differential, and F is F([x, y, z]). T ), Represents the basis coefficients of the spherical harmonic function;

[0054] The result obtained by calculating the spherical harmonic transformation is:

[0055] Step 102: Calculate the standard correlation function SC of the frequency domain signals of the two input point clouds. R (g), as shown in the following formula:

[0056] SC R (g)=∫ x,y,z MSA(g)TAdxdydz

[0057]

[0058]

[0059] In the formula: MSA(g) is... Let (α, β, γ) represent the spherical harmonic basis coefficient matrix of the model point cloud, (α, β, γ) represent the Euler angles corresponding to the rotation, and km and l represent subscripts and superscripts, which will be eliminated during summation. Denotes the Wigner function, d k m l Let represent the Wigner distribution function.

[0060] Step 103: Find all local maxima of the standard correlation function, and calculate the radius of the search domain based on the error function value of each local maximum. For a given local maximum, the radius of the corresponding search domain is... Where E is the size of the optimization function corresponding to the local maximum, and D is the restricted distance of the point cloud outlier in Hypothesis 2.

[0061] Step 2: Based on the new evaluation criteria, assign a search order to the several point cloud search domains determined in Step 1, thereby finding the optimal solution more quickly. This invention calculates this evaluation criterion by analyzing the hit rate and search time of the search domains, wherein:

[0062] The formula for estimating hit rate is shown below:

[0063]

[0064] In the formula: LB represents the lower bound, N in Represents the number of point cloud matching pairs, residual i Let represent the matching error of the i-th matching pair, radius represent the radius of the search domain, C is a constant (2√3 in this embodiment), caf is the centering decay factor; r represents the radius of the search domain, f(·) represents the selectable function (in this embodiment, f(·) is a linear function); E * This indicates the size of the current optimal optimization function.

[0065] The formula for estimating search time is as follows:

[0066] F t (r)=r o *(LB-E * )

[0067] In the formula: o is an optional parameter.

[0068] This invention uses hit rate estimation / search time estimation as the evaluation criteria to arrange the order of searches in the search domain.

[0069] Step 3: Using the branch and bound method, search sequentially through several point cloud search domains in the order determined in Step 2. During the search, if the optimal solution is not within the current search domain, the branch and bound method will exclude the current search domain. If the optimal solution is within the current search domain, the branch and bound method will find the optimal solution.

[0070] The speed of the branch and bound method is highly dependent on the accuracy of the bounds. Compared to other existing methods, this invention proposes a more accurate residual estimation method, making the lower bound of the branch and bound method more accurate. Furthermore, based on the two assumptions proposed in this invention, a new lower bound is introduced to accelerate the operation of the branch and bound method.

[0071] The branch and bound method involves the following steps: For each search domain, calculate the error function value at the center of the current search domain, and estimate the lower bound of the function for the current search domain using this error function value. If the lower bound is greater than the existing best result, then exclude the search domain; otherwise, recursively split the current search domain into smaller search domains. When the size of the search domain is smaller than the target accuracy requirement, the center point of the search domain is the optimal solution. This recursive process guarantees that if the function has one and only one optimal solution, then the branch and bound method will eventually converge to a single smallest search domain, which is the optimal solution.

[0072] The calculation method for the lower bound has been explained above and will not be repeated here. Furthermore, this invention introduces a check condition during the branch-bound process (corresponding to the previous assumptions):

[0073] N in ≤H T *N

[0074] In the formula, H T This refers to the proportion of matching point cloud pairs in the assumption, where N represents the size of the input point cloud. This additional check condition accelerates the search process of the branch and bound method.

[0075] Experimental results on the Stanford dataset and the RGB-D object dataset show that the proposed algorithm is on average 53.6% and 31.3% faster than the state-of-the-art non-correspondence-based and correspondence-based implementations, respectively.

Claims

1. A method for global registration based on spherical harmonic transform and iterative optimization, characterized in that, Assumptions: First, a minimum percentage of correctly matched pairs in the input point cloud is predefined, representing a limit on the number of outliers in the point cloud; second, the distance threshold used to determine correctly matched pairs is preset by the scene, representing a limit distance for outliers in the point cloud. Then, the global registration method includes the following steps: Step 1: Align the point cloud in the frequency domain using spherical harmonics based on the second point assumption, and determine several small point cloud search domains based on the peak value of the phase difference and the first point assumption, including the following steps: Step 101: Voxelize the two input point clouds and perform a spherical harmonic transform to obtain the frequency domain signals of the two input point clouds, where: The voxelization formula used when voxelizing two input point clouds is shown below: wherein: is the resulting voxelized point cloud; , is the input point cloud; denotes the maximum norm of a point in the input point cloud; The spherical harmonic transformation is calculated using the following formula: In the formula: Represents the coordinates of a point. The superscript indicates a spherical harmonic function. Indicates the subscript of a spherical harmonic function. Represents a spherical harmonic function. This represents the accompanying Legendre polynomial. To represent the differential, That is , Represents the basis coefficients of the spherical harmonic function; By computing the spherical harmonic transform in this way, the result is ; Step 102, compute the standard cross-correlation function of the frequency domain signals of the two input point clouds As shown in the following formula: In the formula: That is , This represents the spherical harmonic basis coefficient matrix of the model point cloud. This represents the Euler angle corresponding to the rotation. Represents the Wigner function. Represent the Wigner distribution function; Step 103: Find all local maxima of the standard correlation function, and calculate the radius of the search domain based on the error function value of each local maximum. For a given local maximum, the radius of the corresponding search domain is... ,in, The size of the optimization function corresponding to this local maximum. The distance constraint for point cloud outliers in Hypothesis 2; Step 2: Assign a search order to the point cloud search domains determined in Step 1 by analyzing the hit rate and search time of the search domains, where: The formula for estimating hit rate is shown below: In the formula: Indicates the lower bound. , This indicates the number of point cloud matching pairs. Indicates the first The matching error of each matching pair Indicates the radius of the search domain. It is a constant. It is a centralized decay factor; Indicates the radius of the search domain. Indicates an optional function; This indicates the size of the current optimal optimization function; The formula for estimating search time is as follows: In the formulae: are optional parameters; Step 3: Using the branch and bound method, search sequentially in several point cloud search domains according to the order determined in Step 2. During the search, if the optimal solution is not within the current search domain, the branch and bound method will exclude the current search domain. If the optimal solution is within the current search domain, the branch and bound method will find the optimal solution. The following check condition is introduced during the branch and bound process: In the formula, This refers to the proportion of point cloud matching pairs in the first hypothesis. This indicates the size of the input point cloud.

2. A global registration method based on spherical harmonic transform and iterative optimization as claimed in claim 1, wherein, The branch and bound method includes the following steps: For each search domain, calculate the error function value of the center of the current search domain, and estimate the lower bound of the function of the current search domain through the error function value: If the lower bound of the function is greater than the existing best result, then exclude the search domain; otherwise, recursively split the current search domain into smaller search domains and perform this process. When the size of the search domain is smaller than the target accuracy requirement, the center point of the search domain is the optimal solution.

Images