A method for 3D reconstruction of propeller blade surface based on point cloud data

Abstract

A method for three-dimensional reconstruction of propeller blade surfaces based on point cloud data. A truncated least squares cost function is used to redefine the relationship between the point cloud to be registered and the target point cloud; an adaptive voting algorithm is used to estimate the scaling ratio S through translation and rotation invariant measurement methods; a semidefinite relaxation algorithm is used to estimate the rotation matrix R through scaling ratio S and translation invariant measurement methods; an adaptive voting algorithm is used to estimate the translation matrix t through scaling ratio S, rotation matrix R and truncated least squares function; point cloud registration is achieved based on scaling ratio S, rotation matrix R and translation matrix t; the actual propeller blade is compared with the errors between the registered geometric parameters to determine whether they are all less than a set threshold; finally, the registered propeller blade point cloud data is subjected to bicubic NURBS surface fitting. The present invention effectively eliminates the influence of most false corresponding points in the registration process, and can significantly improve the efficiency and accuracy of reconstructing the propeller blade surface.

Description

Technical Field

[0001] The invention relates to the field of three-dimensional model reconstruction, and in particular to a method for three-dimensional reconstruction of a propeller blade surface based on point cloud data. Background Art

[0002] With the development of computer vision and the advancement of hardware technology, 3D reconstruction technology has been increasingly used in our daily lives. Registration algorithm is one of the important steps that affect the accuracy of 3D reconstruction. In the process of 3D reconstruction, due to the limitation of the position of the instrument, the scanner cannot fully collect the 3D information of the object surface at one time. Therefore, to achieve 3D reconstruction of the object, the key is to convert the point cloud data collected from different angles or different coordinate systems to a unified angle or reference system, and the key to solving this conversion relationship is to achieve point cloud registration.

[0003] The iterative closest point (ICP) algorithm is considered a milestone in point cloud registration. However, the ICP algorithm is limited by the large amount of computation, slow convergence speed, and the requirement that the two point clouds to be registered have a high degree of overlap, so the registration effect is often not very ideal. Clark et al. used Riemann optimization to align point clouds as continuous functions. Le et al. used semidefinite relaxation to generate assumptions for random methods. But all these methods cannot guarantee global optimality. Therefore, how to accurately, comprehensively, and highly accurately complete the registration of propeller blade point cloud data and complete surface reconstruction is an urgent problem to be solved. Summary of the invention

[0004] In view of the above-mentioned deficiencies in the prior art, an object of the present invention is to provide a method for three-dimensional reconstruction of a propeller blade surface based on point cloud data.

[0005] A method for three-dimensional reconstruction of a propeller blade surface based on point cloud data comprises the following steps:

[0006] Step 1: Use the truncated least squares cost function to redefine the relationship between the point cloud to be registered and the target point cloud;

[0007] Step 2: Use an adaptive voting algorithm to estimate the scaling factor S using translation and rotation invariant measurements.

[0008] Step 3: Use the semidefinite relaxation algorithm to estimate the rotation matrix R by scaling S and translation-invariant measurement method;

[0009] Step 4: Adopt an adaptive voting algorithm to estimate the translation matrix t by using the scaling factor S, the rotation matrix R and the truncated least squares function;

[0010] Step 5: According to the scaling ratio S, the rotation matrix R and the translation matrix t, the point cloud is registered to obtain the complete propeller blade point cloud data;

[0011] Step 6: Calculate the geometric parameters of the propeller blades, compare the actual propeller blades with the errors between the geometric parameters after registration, and determine whether they are all less than the set threshold. If so, proceed to the next step 7. Otherwise, return to step 2 and reset the noise boundary β i ;

[0012] Step 7: Reconstruct the bicubic NURBS surface of the registered propeller blade point cloud data.

[0013] As an improvement, the formula for redefining the relationship between the point cloud to be registered and the target point cloud in step 1 is as follows:

[0014] First set a boundary β for the noise i , so that ||ε i ||≤β i , where ε i is a noise point, and the following truncated least squares cost function is used to redefine the point cloud to be registered and target point cloud Relationship:

[0015]

[0016] In the formula, a i is the i-th point to be registered; b i is the i-th target point cloud; S, R and t are the scaling ratio, rotation matrix and translation matrix calculated in steps 2, 3 and 4 respectively; N is the total number of point cloud data to be registered; is a threshold value set according to the noise. Usually This function makes the point cloud set insensitive to a large number of false corresponding points (including noise points) during the registration process.

[0017] As an improvement, the steps of the adaptive voting algorithm in step 2 and step 4 are as follows:

[0018] Step 1: Given the s variable, define the consistent set For any s, there are at most 2K-1 different non-empty consistent sets;

[0019] Step 2: Given a symbol s k , that is, every s in the kth interval satisfies in, It refers to the ratio of the distance between two points in the target point cloud to the distance between two points in the point cloud to be registered;

[0020] Step 3: Because the consensus set does not change in a single interval, the center position of the interval is calculated;

[0021] Step 4: Calculate the cost of each consensus set and return the one with the minimum cost as the optimal set.

[0022] As an improvement, the step of estimating the scaling factor S in step 2 is as follows:

[0023] Step i: Compute the norm of the translation-invariant measure:

[0024] ||b ij ||=||SRa ij +ε ij ||

[0025] Since the measurement noise is bounded, we can get ||b ij ||=||SRa ij ||+ε ij ;

[0026] Step ii: Divide both sides by ||a ij ||, then we get

[0027] Step iii: The scaling ratio S can be calculated according to the following formula:

[0028]

[0029] As an improvement, the steps for estimating the rotation matrix R in step 3 are as follows:

[0030] Step 3-1: Define variables matrix

[0031] Step 3-2: Use unit quaternion to convert the rotation truncated least squares cost function into mixed integer programming and quaternion form, and finally obtain non-convex rank constraint programming. On this basis, cancel the rank-1 constraint and add redundant constraints, and then obtain the relaxed model of semidefinite programming with redundant constraints as follows:

[0032]

[0033]

[0034] In the formula, tr(QZ)=tr(Qxx T )=x T Qx; [Z] 00 Refers to the upper left 4×4 matrix in Z; [Z] kk refers to the kth diagonal matrix in Z;

[0035] Step 3-3: Using the translation invariant measurement principle, the i-th point a in the point cloud to be registered and the target point cloud i and b i and the jth point a j and bj Eliminate the corresponding two points a by the following formula i and b i Translation distance:

[0036] b j -b i =SR(a j -a i )+(ε j -ε i )

[0037] Step 3-4: The rotation matrix R can be calculated according to the following formula:

[0038]

[0039] In the formula, a k =a j -a i The physical meaning of b is the distance between two points in the point cloud to be registered. k =b j -b i The physical meaning represented is the distance between two points in the target point cloud, ε k =ε j -ε i The physical meaning represented is the distance between two noise points.

[0040] As an improvement, the translation matrix t in step 4 is obtained by bringing the scaling ratio S and the rotation matrix R obtained in steps 2 and 3 into the truncated least squares cost function described in claim 2, and combining it with the adaptive algorithm in claim 3. The specific formula is as follows:

[0041]

[0042] As an improvement, the bicubic NURBS surface fitting in step 7 includes:

[0043] Step 7-1: Determine the node vector;

[0044] Step 7-2: Solve the fitting boundary conditions and correct them;

[0045] Step 7-3: Inversely calculate the control points of the interpolation curve and obtain the parameter coordinates P(u,v) of each point on the surface to generate the propeller blade bicubic NURBS surface.

[0046] Beneficial effects:

[0047] Compared with the prior art, the present invention discloses a method for three-dimensional reconstruction of propeller blade surface based on point cloud data, which has the following advantages:

[0048] (1) The method for 3D reconstruction of propeller blade surface based on point cloud data of the present invention uses a truncated least squares cost function to redefine the relationship between the point cloud to be registered and the target point cloud, thereby solving the problem that there are too many false corresponding points in point cloud data and they are difficult to handle, making the registration process insensitive to a large number of false corresponding points (including noise points), thereby reducing the data preprocessing process and improving the efficiency of point cloud registration;

[0049] (2) The method for 3D reconstruction of propeller blade surfaces based on point cloud data of the present invention uses an adaptive voting algorithm to accurately solve the truncated complete least squares problem of scalar examples in polynomial time, and can efficiently estimate the scaling factor S and the translation matrix t; the accuracy of the estimation result far exceeds the singular value decomposition in the traditional ICP algorithm, thereby further improving the accuracy and precision of 3D reconstruction;

[0050] (3) The method for three-dimensional reconstruction of propeller blade surfaces based on point cloud data of the present invention uses a semidefinite relaxation algorithm to ensure global optimality when solving convex feasible sets in polynomial time. Semidefinite programming is a generalization of linear programming, and its duality theory is weaker than that of linear programming. Therefore, many nonlinear convex optimization problems can be reduced to semidefinite programming for solution. Its solution conditions are more general and the solution results are more accurate, thereby further reducing the complexity of three-dimensional reconstruction and improving the accuracy of three-dimensional reconstruction;

[0051] (4) The three-dimensional reconstruction method of the propeller blade surface based on point cloud data of the present invention specifically verifies the model after point cloud registration according to the geometric parameters of the propeller blade, and quantifies the actual error value, thereby improving the efficiency of surface reconstruction and shortening the research cycle. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 It is a schematic flow chart of the three-dimensional reconstruction process of the present invention;

[0053] Figure 2 It is a schematic diagram of a multi-view point cloud of a propeller;

[0054] Figure 3 It is a multi-view scanning propeller viewing angle pattern;

[0055] Figure 4 is the cross-sectional geometry of the propeller blade;

[0056] Figure 5 It is a complete propeller blade point cloud data image after registration, where (a) is a top view of the complete propeller blade point cloud, and (b) is a front view of the complete propeller blade point cloud;

[0057] Figure 6 This is the propeller blade image after bicubic NURBS surface fitting. DETAILED DESCRIPTION

[0058] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution in the embodiments of the present invention will be clearly described below in conjunction with the accompanying drawings and specific embodiments. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0059] Figure 1 The present invention is a schematic flow chart of a method for three-dimensional reconstruction of a propeller blade surface based on point cloud data according to an embodiment of the present invention, which includes any one or more of the following steps:

[0060] The implementation premise of the method of the present invention is that the on-site staff uses a three-dimensional laser scanner to scan around the propeller blade from multiple directions, and obtains a point cloud in each direction, and obtains a total of Figure 2 To ensure that there is an overlapping area between two adjacent point clouds, the initial scanning direction of the scanner and the moving direction (i.e. Figure 3 The viewing angle direction) is 135°.

[0061] Steps 1-6 will take point cloud sets A and B as examples to describe the specific implementation methods in detail:

[0062] Step 1: Use the truncated least squares cost function to redefine the relationship between the point cloud to be registered and the target point cloud;

[0063] Assume that the point cloud to be registered is Yes, the target point cloud is (m, n are the number of point clouds respectively), ε i It is noise:

[0064] Then the relationship between point clouds A and B is as follows:

[0065]

[0066] In order to make point clouds A and B insensitive to a large number of false corresponding points (including noise points) during the registration process, a boundary β is first set for the noise points. i , let ||ε i ||≤β i , and use the following function to redefine the relationship between the point cloud A to be registered and the target point cloud B:

[0067]

[0068] In the formula, a i is the i-th point to be registered; b iis the i-th target point cloud; S, R and t are the scaling ratio, rotation matrix and translation matrix calculated in steps 2, 3 and 4 respectively; N is the total number of point cloud data to be registered; is a threshold value set according to the noise. Usually This function makes the point cloud set insensitive to a large number of false corresponding points (including noise points) during the registration process.

[0069] Step 2: Use an adaptive voting algorithm to estimate the scaling factor S using translation and rotation invariant measurements.

[0070] The principle of the adaptive voting algorithm is as follows:

[0071] Step 1: Given the s variable, define the consistent set For any s, there are at most 2K-1 different non-empty consistent sets;

[0072] Step 2: Given a symbol s k , that is, every s in the kth interval satisfies in, It refers to the ratio of the distance between two points in the target point cloud to the distance between two points in the point cloud to be registered;

[0073] Step 3: Because the consensus set does not change in a single interval, the center position of the interval is calculated;

[0074] Step 4: Calculate the cost of each consensus set and return the one with the minimum cost as the optimal set.

[0075] The principle of the translation and rotation invariant measurement method is as follows:

[0076] Step i: Compute the norm of the translation-invariant measure:

[0077] ||b ij ||=||SRa ij +ε ij ||

[0078] Since the measurement noise is bounded, we can get ||b ij ||=||SRa ij ||+ε ij ;

[0079] Step ii: Divide both sides by ||a ij ||, then we get

[0080] Step iii: The scaling ratio S can be calculated according to the following formula:

[0081]

[0082] Step 3: Use the semidefinite relaxation algorithm to estimate the rotation matrix R by scaling S and translation-invariant measurement method;

[0083] Step 3-1: Define variables matrix

[0084] Step 3-2: Use unit quaternion to convert the rotation truncated least squares cost function into mixed integer programming and quaternion form, and finally obtain non-convex rank constraint programming. On this basis, cancel the rank-1 constraint and add redundant constraints, and then obtain the relaxed model of semidefinite programming with redundant constraints as follows:

[0085]

[0086]

[0087] In the formula, tr(QZ)=tr(Qxx T )=x T Qx; [Z] 00 Refers to the upper left 4×4 matrix in Z; [Z] kk refers to the kth diagonal matrix in Z;

[0088] Step 3-3: Using the translation invariant measurement principle, the i-th point a in the point cloud to be registered and the target point cloud i and b i and the jth point a j and b j Eliminate the corresponding two points a by the following formula i and b i Translation distance:

[0089] b j -b i =SR(a j -a i )+(ε j -ε i )

[0090] Step 3-4: The rotation matrix R can be calculated according to the following formula:

[0091]

[0092] In the formula, a k =a j -a i The physical meaning of b is the distance between two points in the point cloud to be registered. k =b j -b i The physical meaning represented is the distance between two points in the target point cloud, ε k =εj -ε i The physical meaning represented is the distance between two noise points.

[0093] Step 4: Adopt an adaptive voting algorithm to estimate the translation matrix t by using the scaling factor S, the rotation matrix R and the truncated least squares function;

[0094] The translation matrix t is obtained by substituting the scaling ratio S and the rotation matrix R obtained in steps 2 and 3 into the truncated least squares cost function described in step 1 and combining it with the adaptive algorithm in step 2. The specific formula is as follows:

[0095]

[0096] Step 5: According to the scaling ratio S, the rotation matrix R and the translation matrix t, the point cloud is registered to obtain the complete propeller blade point cloud data;

[0097] Step 6: According to Figure 4 The geometric parameters of the propeller blades can be calculated through the geometric relationship, and the errors between the actual propeller blades and the registered geometric parameters can be compared to determine whether they are all less than the set threshold. If so, proceed to the next step 7. Otherwise, return to step 2 and reset the noise boundary β i .

[0098] Specifically, the propeller geometric parameter calculation formula and setting threshold are as follows:

[0099]

[0100] Assumptions Figure 4 The coordinates of point E are (x1, y1), and the coordinates of point F are (x n ,y n )

[0101] Then in the radius: x'=x n -x1,y'=y n -y1

[0102] Disk ratio: g(ξ i )=y xi,k -y xi,j

[0103] Pitch Angle:

[0104] Step 7: Register the propeller blade point cloud data (see Figure 5 ) performs bicubic NURBS surface fitting, including the following steps:

[0105] NURBS surface is based on the tensor product form of NURBS curves, so the bicubic NURBS surface can be expressed as:

[0106]

[0107] Step 7-1: Determine the node vector;

[0108] When a row (or column) of sample points is used as the type value point and a bicubic NURBS curve is used for fitting, the repetition degree at both ends is set to 4. Even if the first and last points of the curve pass through the type value point, the middle point P i There is only one parameter u i+3 Correspondingly, the node vector u can be obtained by using the cumulative chord length method:

[0109]

[0110] Where: n is the number of type value points; P i is the type value point; S is the total length of the chord,

[0111] Step 7-2: Solve the fitting boundary conditions and correct them;

[0112] Common boundary conditions include free endpoint conditions, endpoint tangent vector conditions, parabola conditions, etc. The endpoint tangent vector conditions are used here. Taking the node vector u as the parameter, the high-order curve passing through multiple shape value points at the near end is expressed as a parametric formula:

[0113] R(u)=[X(u),Y(u),Z(u)]

[0114] Then the tangent vector equation is:

[0115] T(u)=[X'(u),Y'(u),Z'(u)]

[0116] Where: X, Y, Z are the three coordinates of the shape point, and u is the corresponding node vector.

[0117] Then the endpoint tangent vector of the fitting curve can be approximated by the endpoint tangent vector of the higher-order curve. Considering that the endpoint tangent vector of the higher-order curve may not be on the tangent plane passing through the corresponding endpoint, its projection on the tangent plane passing through the endpoint can be taken as the approximate tangent vector as a correction.

[0118] Step 7-3: Inversely calculate the control points of the interpolation curve and obtain the parameter coordinates P(u,v) of each point on the surface;

[0119] In order to make the interpolation calculation more concise and efficient, we introduce the matrix representation of the i-th segment of the NURBS curve:

[0120]

[0121] Where: T = [1, t1, t2, t3]; M is a 4×4 coefficient matrix;

[0122] D=[ω i d i ,ω i+1 d i+1 ,ω i+2 d i+2 ,ω i+3 d i+3 ] T W=[ω i ,ω i+1 ,ω i+2 ,ω i+3 ] T

[0123] According to the characteristic that the NURBS curve segment is connected at the shape value point, the above formula can be further simplified to obtain:

[0124] [m 11 ω i ,m 12 ω i+1 ,m 12 ω i+2 ][d i ,d i+1 ,d i+2 ] T =[m 11 ω,m 12 ω i+1 ,m 12 ω i+2 ]R i

[0125] From the above properties, we can get n+1 equations, but there are n+3 unknowns. Adding the endpoint tangent vector conditions obtained in step 2, we can solve all the controlling vertices of the curve.

[0126] According to the distribution characteristics of the sampling points, the control vertex Q of the curve family along the U direction can be obtained by repeatedly calculating according to the above steps. ij , and weight factors, and then Q ij The control vertex V of the NURBS surface can be obtained by calculating along the V direction according to the above method. ij , and weight factors, and then substitute them into P(u,v) to obtain the parameter coordinates of each point on the surface, thereby obtaining the propeller blade bicubic NURBS surface. The result is as follows Figure 6 shown.

[0127] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the technical principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A method for three-dimensional reconstruction of propeller blade surface based on point cloud data, characterized in that: The following steps are involved: Step 1: Use the truncated least squares cost function to redefine the relationship between the point cloud to be registered and the target point cloud; Step 2: Use an adaptive voting algorithm to estimate the scaling factor S using translation and rotation invariant measurements. Step 3: Use the semidefinite relaxation algorithm to estimate the rotation matrix R by scaling S and translation-invariant measurement method; Step 4: Adopt an adaptive voting algorithm to estimate the translation matrix t by using the scaling factor S, the rotation matrix R and the truncated least squares cost function; Step 5: According to the scaling ratio S, the rotation matrix R and the translation matrix t, the point cloud is registered to obtain the complete propeller blade point cloud data; Step 6: Calculate the geometric parameters of the propeller blades, compare the actual propeller blades with the errors between the geometric parameters after registration, and determine whether they are all less than the set threshold. If so, proceed to the next step 7. Otherwise, return to step 2 and reset the noise boundary β i ; Step 7: Reconstruct the bicubic NURBS surface of the registered propeller blade point cloud data; The formula for redefining the relationship between the point cloud to be registered and the target point cloud in step 1 is as follows: First set a boundary β for the noise i , so that ||ε i ||≤β i , where ε i is a noise point. Use the following function to redefine the point cloud to be registered and target point cloud Relationship: In the formula, a i is the i-th point to be registered; b i is the i-th target point; S, R and t are the scaling ratio, rotation matrix and translation matrix calculated in steps 2, 3 and 4 respectively; N is the total number of point cloud data to be registered; is a threshold value set according to the noise. Usually This function makes the point cloud set insensitive to a large number of false corresponding points (including noise points) during the registration process.

2. The method for three-dimensional reconstruction of propeller blade surface based on point cloud data according to claim 1, characterized in that: The steps of the adaptive voting algorithm in step 2 and step 4 are as follows: Step 1: Given the s variable, define the consistent set For any s, there are at most 2K-1 different non-empty consistent sets; Step 2: Given a symbol s k , that is, every s in the kth interval satisfies in, It refers to the ratio of the distance between two points in the target point cloud to the distance between two points in the point cloud to be registered; Step 3: Because the consensus set does not change in a single interval, the center position of the interval is calculated; Step 4: Calculate the cost of each consensus set and return the one with the minimum cost as the optimal set.

3. The method for three-dimensional reconstruction of propeller blade surface based on point cloud data according to claim 1, characterized in that: The steps to estimate the scaling factor S in step 2 are as follows: Step i: Compute the norm of the translation-invariant measure: ||b ij ||(||SRa ij +ε ij || Since the measurement noise is bounded, we can get ||b ij ||=||SRa ij ||+ε ij ; Step ii: Divide both sides by ||a ij ||, then we get Step iii: The scaling ratio S can be calculated according to the following formula:

4. The method for three-dimensional reconstruction of propeller blade surface based on point cloud data according to claim 1, characterized in that: The steps to estimate the rotation matrix R in step 3 are as follows: Step 3-1: Define variables matrix Step 3-2: Use unit quaternion to transform the rotation truncated least squares cost function into mixed integer programming and quaternion form, and finally obtain non-convex rank constraint programming, cancel rank-1 constraint, add redundant constraint, and obtain the relaxed model of semidefinite programming with redundant constraint as follows: In the formula, tr(QZ)=tr(Qxx T )=x T Qx; [Z] 00 Refers to the upper left 4×4 matrix in Z; [Z] kk refers to the kth diagonal matrix in Z; Step 3-3: Using the translation invariant measurement principle, the i-th point a in the point cloud to be registered and the target point cloud i and b i and the jth point a j and b j Eliminate the corresponding two points a by the following formula i and b i Translation distance: b j -b i =SR(a j -a i )+(ε j -ε i ) Step 3-4: The rotation matrix R can be calculated according to the following formula:

5. The method for three-dimensional reconstruction of propeller blade surface based on point cloud data according to claim 1, characterized in that: The translation matrix t in step 4 is obtained by bringing the scaling ratio S and the rotation matrix R obtained in steps 2 and 3 into the truncated least squares cost function described in claim 2, and combining it with the adaptive algorithm in claim 3. The specific formula is as follows:

6. The method for three-dimensional reconstruction of propeller blade surface based on point cloud data according to claim 1, characterized in that: The bicubic NURBS surface fitting in step 7 includes: Step 7-1: Determine the node vector; Step 7-2: Solve the fitting boundary conditions and correct them; Step 7-3: Inversely calculate the control points of the interpolation curve and obtain the parameter coordinates P(u,v) of each point on the surface to generate the propeller blade bicubic NURBS surface.

Images