A method for matching a three-dimensional point cloud to a reconstruction model

Abstract

The application relates to a three-dimensional point cloud matching reconstruction model method, which comprises the following steps: S1, normal estimation on the three-dimensional point cloud; S2, marking of the normal abnormal value according to the overall normal distribution; S3, near-point iteration on the normal abnormal value and normal matching; and S4, iteration matching of the normal abnormal value to the normal point cloud position, so as to realize point cloud reconstruction. The three-dimensional point cloud matching reconstruction model method has the following advantages: through the matching reconstruction of the three-dimensional point cloud, the outlying part in the three-dimensional point cloud can be effectively restored to the original position, and the original three-dimensional shape can be recovered.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of data processing, in particular to a method for matching and reconstructing a model of three-dimensional point cloud. BACKGROUND

[0002] In the prior art, scanning devices such as sonar and laser radar have been gradually applied to underwater three-dimensional measurement of caissons, but there are the following problems in their application: due to the limited measurement conditions in the well, the connecting device will inevitably cause irregular point cloud of the well wall to be measured due to device shaking or installation deviation when measuring the caisson. SUMMARY

[0003] An object of the present application is to provide a method for matching and reconstructing a model of three-dimensional point cloud, which can match and reconstruct irregular point cloud.

[0004] The technical scheme adopted by the present application is as follows: a method for matching and reconstructing a model of three-dimensional point cloud, comprising the following steps: S1, normal estimation of three-dimensional point cloud; S2, marking normal abnormal values according to the overall normal distribution; S3, near-point iteration of normal abnormal values and normal matching; S4, iteration matching of points of normal abnormal values to normal point cloud positions to realize point cloud reconstruction.

[0005] In some embodiments of the present application, in step S1, a normal estimation function is used to calculate the normal vector for each data point of the three-dimensional point cloud, points within a radius A are found, a covariance matrix is constructed using the coordinates of the points within the radius A, eigenvalue decomposition is performed on the covariance matrix, and three eigenvalues and their corresponding eigenvectors are obtained. The eigenvector corresponding to the smallest eigenvalue is the normal vector Ni=(Nx i , Ny i , Nz i ).

[0006] Further, the radius A is 3 times the average distance of the point cloud.

[0007] Further, step S2 comprises: S21, slicing the point cloud horizontally according to the Z coordinate of the point cloud; S22, calculating the geometric centroid point of each horizontal slice point cloud; S23, constructing a horizontal slice plane coordinate system XOY with the slice geometric centroid point as the origin, the X axis direction of the horizontal slice plane coordinate system being the same as the X axis direction of the original point cloud coordinate, the Y axis direction of the horizontal slice plane coordinate system being the same as the Y axis direction of the original point cloud coordinate, and calculating the slice coordinate system coordinates (x i , y i ) of each point; S24, calculating the normal difference Δ NThe normal abnormal value is greater than the threshold value 0.3, and the point cloud of the normal abnormal value is a normal abnormal point cloud.

[0008] Further, in the step S3, the normal abnormal point cloud is taken as a source point cloud, and the normal point cloud is taken as a target point cloud, for each point in the normal abnormal point cloud, a neighborhood point of the point in the normal abnormal point cloud is found, and the neighborhood point is registered to a neighborhood point position of the normal point cloud by a generalized iterative closest point registration and Levenberg-Marquardt error optimization fusion algorithm.

[0009] Further, the step S3 includes: S31, initializing parameters, setting an initial transformation matrix T as a unit matrix, performing initial alignment of the source point cloud and the target point cloud, setting a damping factor lambda initial value = 10 6 , setting a maximum iteration number = 100 times, and setting a convergence threshold value = 10 -3 When the parameter update amount is less than the threshold value, it is considered that the algorithm has converged; S32, matching the closest neighbor points, for each point in the source point cloud, a closest point in the target point cloud is found as a corresponding point; S33, minimizing the error, the error is minimized by an optimization method to obtain a new rigid body transformation matrix;

[0010] S34, applying the calculated rigid body transformation matrix T to the source point cloud to obtain the normal abnormal point cloud after registration with the normal point cloud.

[0011] Further, the step S33 includes:

[0012] S331, calculating a residual vector r: calculating the residual between the source point cloud and the corresponding target point cloud, for each pair of corresponding points p i and q i , the residual r i can be defined as:

[0013] r i =T(p i )-q i =Rp i +t-q i (i=1,2,..,n)

[0014]

[0015] Wherein, T is a transformation matrix, which includes a rotation R and a translation matrix t, and the residual vector r = (r1 … r n ) T ;

[0016] S332, constructing an error function: the error function E is the sum of squares of the residual vector, that is:

[0017]

[0018] S333, calculating the Jacobian matrix J: calculating the partial derivative of the residual with respect to the transformation parameters, constructing the Jacobian matrix J, each element of the Jacobian matrix J represents the rate of change of the residual with respect to each parameter,

[0019]

[0020] S334, updating the transformation matrix: using the Levenberg-Marquardt method to iteratively optimize the error function, finding the transformation matrix T that minimizes the error function, in each iteration, calculating the Jacobian matrix J and the residual vector r, then calculating the Gauss-Newton step, updating the transformation matrix T;

[0021] S335, adjusting the damping factor: in each iteration, adjusting the damping factor lambda value according to the change of the objective function value, if the objective function value decreases in the iteration, then: lambda new = lambda / 10, making the algorithm closer to the Gauss-Newton method, if the objective function value increases, then: lambda new = lambda*10, to reduce the step size and avoid excessive update;

[0022] S336, checking the convergence condition: if the parameter update amount is less than the preset threshold, or the number of iterations reaches the maximum value, the algorithm ends, that is, the target point cloud position is found.

[0023] Further, in the step S334, the Gauss-Newton step is calculated: Δx GN = (J T J) -1 J T r, the damping step is calculated: Δx LM = (J T J+λI) -1 J T , the transformation matrix is updated: T new = T+Δx LM , wherein J is the Jacobian matrix, I is the unit matrix, lambda is the adjustment factor, Δx is the parameter update amount, and r is the residual vector calculated in the previous step.

[0024] In some embodiments of the present application, step S0 is further included, which is located before step S1, and the step S0 is denoising the three-dimensional point cloud.

[0025] In some embodiments of the present application, step S5 is further included, which is located after step S4, and the step S5 is inputting the reconstructed three-dimensional point cloud into matlab to obtain three-dimensional perspective view and geometric dimension information.

[0026] The method for matching and reconstructing a model of a three-dimensional point cloud has the following advantages: through matching and reconstruction of the three-dimensional point cloud, the outlying part in the three-dimensional point cloud can be effectively restored to the original position, and the original three-dimensional shape can be recovered. BRIEF DESCRIPTION OF DRAWINGS

[0027] Figure 1 is a flowchart of embodiment 1 of the present application;

[0028] Figure 2 is a flowchart of step S3 in embodiment 1 of the present application;

[0029] Figure 3 is a schematic diagram of the original caisson point cloud in embodiment 1 of the present application;

[0030] Figure 4 is a schematic diagram of the normal estimation result in embodiment 1 of the present application;

[0031] Figure 5 is a schematic diagram of the matched and reconstructed caisson point cloud in embodiment 1 of the present application;

[0032] Figure 6 is a structural schematic diagram of the multi-beam sounding system in embodiment 1 of the present application;

[0033] Figure 7 is a multi-beam sounding principle diagram in embodiment 1 of the present application;

[0034] Figure 8 is a beam forming schematic diagram of the multi-beam system in embodiment 1 of the present application.

[0035] In the figure: 1, multi-beam sonar system; 2, rotating electromechanical system; 3, mounting bracket system; 4, cable; 5, sonar main machine; 6, computer. DETAILED DESCRIPTION

[0036] To further illustrate the technical means and effects adopted by the present application to achieve the predetermined invention purpose, the specific embodiments, structures, features and effects according to the present application are described in detail as follows in combination with the drawings and preferred embodiments.

[0037] Embodiment 1:

[0038] The method for matching and reconstructing a model of a three-dimensional point cloud provided in this embodiment, as shown in Figure 1 , Figure 2 includes the following steps:

[0039] S0, denoising the three-dimensional point cloud;

[0040] Specifically, the acquired multi-beam data is subjected to attitude correction and draft correction by using multi-beam post-processing software (CARIS HIPS AND SIPS), and underwater three-dimensional point cloud data of the caisson is obtained, and noise points are removed by using point cloud processing software CloudCompare.

[0041] S1, normal estimation is performed on the three-dimensional point cloud;

[0042] Specifically, for each data point of the three-dimensional point cloud, a normal estimation function is used to calculate a normal vector, points within a radius A are found, a covariance matrix is constructed using the coordinates of the points within the radius A, eigenvalue decomposition is performed on the covariance matrix, three eigenvalues and their corresponding eigenvectors are obtained, and the eigenvector corresponding to the smallest eigenvalue is the normal vector Ni=(Nx i , Ny i , Nz i ) of the point. The radius A is 3 times the average distance of the point cloud.

[0043] S2, mark the normal abnormal value according to the overall normal distribution;

[0044] Specifically, S21, slice the point cloud horizontally according to the Z coordinate of the point cloud; S22, calculate the geometric centroid point of each horizontal slice point cloud; S23, construct a horizontal slice plane coordinate system XOY with the slice geometric centroid point as the origin, the X axis direction of the horizontal slice plane coordinate system is the same as the X axis direction of the original point cloud coordinate, the Y axis direction of the horizontal slice plane coordinate system is the same as the Y axis direction of the original point cloud coordinate, and the slice coordinate system coordinates (x i , y i ) of each point are calculated; S24, calculate the normal difference Δ N of the horizontal slice, and the points with the normal difference exceeding a threshold value of 0.3 are normal abnormal values, and the point cloud of the normal abnormal values is a normal abnormal point cloud.

[0045] S3, near point iteration is performed on the normal abnormal value, and normal matching is performed;

[0046] Specifically, the normal abnormal point cloud is taken as a source point cloud, and the normal point cloud is taken as a target point cloud, for each point in the normal abnormal point cloud, its neighborhood points in the normal abnormal point cloud are found, and the neighborhood points are registered to the position of the neighborhood points of the normal point cloud, and the steps of the fusion algorithm of the generalized iterative closest point (GICP) registration and Levenberg-Marquardt error optimization are as follows.

[0047] S31, initialize parameters

[0048] Initial transformation matrix: the initial transformation matrix T is set to the unit matrix, and the source point cloud and the target point cloud are initially aligned.

[0049] Set damping factor; damping factor lambda initial value = 10 6 To ensure that the algorithm is closer to gradient descent at the beginning, avoid divergence due to too large step size at the initial iteration.

[0050] Maximum iteration number: set the maximum iteration number = 100 times.

[0051] Convergence threshold: set the convergence threshold = 10 -3 When the parameter update is less than the threshold, it is considered that the algorithm has converged.

[0052] S32, match the nearest neighbor point

[0053] For each point in the source point cloud, find the nearest point in the target point cloud as the corresponding point.

[0054] S33, minimize error

[0055] Minimize the error by optimization method, get new rigid transformation matrix, including rotation matrix and translation matrix.

[0056] Specifically,

[0057] S331, calculate residual vector r

[0058] Calculate the residual between the source point cloud and the corresponding target point cloud, each pair of corresponding points p i and q i , the residual r i can be defined as:

[0059] r i = T(p i )-q i = Rp i +t-q i (i = 1, 2,..,n)

[0060]

[0061] Where T is the transformation matrix, which includes rotation R and translation matrix t.

[0062] Residual vector r = (r1 … r n ) T ;

[0063] S332, construct error function

[0064] Error function E is the sum of the square of residual vector, that is:

[0065]

[0066] S333, calculate Jacobian matrix J

[0067] The partial derivative of the residual with respect to the transformation parameters is calculated, and a Jacobian matrix J is constructed, each element of which represents the rate of change of the residual with respect to each parameter.

[0068]

[0069] S334, update the transformation matrix

[0070] The Levenberg-Marquardt method is used to iteratively optimize the error function, and the transformation matrix T that minimizes the error function is found. In each iteration, the Jacobian matrix J and the residual vector r are calculated, and then the Gauss-Newton step is calculated, and the transformation matrix T is updated.

[0071] (1) Calculate the Gauss-Newton step: Δx GN = (J T J) -1 J T r;

[0072] (2) Calculate the damping step: Δx LM = (J T J+λI) -1 J T ;

[0073] (3) Update the transformation matrix: T new = T+Δx LM ;

[0074] Where J is the Jacobian matrix, I is the identity matrix, lambda is the adjustment factor, Δx is the parameter update, and r is the residual vector calculated in the previous step.

[0075] S335, adjust the damping factor

[0076] In each iteration, the damping factor lambda value is adjusted according to the change of the objective function value, if the objective function value decreases in the iteration: lambda new = lambda / 10, so that the algorithm is closer to the Gauss-Newton method, if the objective function value increases: lambda new = lambda*10, to reduce the step size and avoid excessive updates;

[0077] S336, check the convergence condition

[0078] If the parameter update is less than the preset threshold, or the number of iterations reaches the maximum value, the algorithm ends, that is, the target point cloud position is found, which is the neighborhood point cloud of the normal point cloud.

[0079] S34, apply the transformation

[0080] Apply the rigid transformation matrix T calculated to the source point cloud to obtain the normal abnormal point cloud after registration with the normal point cloud.

[0081] S4, iteratively match the normal abnormal value points to the normal point cloud position to realize point cloud reconstruction of multiple scans.

[0082] S5, input the reconstructed three-dimensional point cloud into matlab to obtain three-dimensional perspective view and geometric dimension information.

[0083] In the embodiment, the normal estimation function adopts the normal estimation function (estimate_normals) in the feature extraction module of the Open3D open source library, and the website of the Open3D open source library is https: / / www.open3d.org / . The algorithm is prior art, and will not be described herein.

[0084] The application can perform matching reconstruction operation on the three-dimensional caisson point cloud collected by devices such as laser radar and sonar, including normal estimation algorithm, point cloud registration algorithm fusing generalized iterative closest point registration and Levenberg-Marquardt least square error optimization, etc. The normal of each well wall and the bottom surface can be calculated through normal estimation calculation for subsequent matching process, and the matching reconstruction of point cloud can be performed according to the calculated normal through iterative closest point calculation. The application can perform matching reconstruction on the three-dimensional point cloud of the caisson scanned multiple times, effectively solving the problems of unmatched well wall point cloud, irregular contour, etc. caused by multiple acquisition of well wall point cloud due to device jitter, rotation error, etc. during point cloud data acquisition.

[0085] In the embodiment, the three-dimensional caisson point cloud data is used, and through matching reconstruction of the three-dimensional point cloud, the outlying part in the three-dimensional caisson point cloud can be effectively restored to the original position, and the original three-dimensional shape of the caisson can be restored, as shown in Figure 3 、 Figure 4 、 Figure 5 .

[0086] The three-dimensional point cloud data can be point cloud data collected by scanning devices such as sonar and laser radar detection, and the three-dimensional caisson point cloud data collection in the embodiment adopts a multi-beam sounding system, as shown in Figure 6As shown, including a multi-beam sonar system 1, a rotating electromechanical system 2, a mounting bracket system 3 and corresponding software, the software includes: multi-beam data acquisition software, Viterbi angle acquisition software, point cloud calculation generation software, the multi-beam sonar system 1 includes NORBIT WBMS multi-beam sonar and its supporting control and acquisition software; the rotating electromechanical system 2 drives the connecting rod and the sonar transducer to realize 360° uniform rotation through the uniform rotation of the rotating motor, and then realizes the omnidirectional measurement of the caisson. When working, the transducer is placed in the appropriate water depth through the bracket, the cable 4 transmits the data collected by the sonar transducer to the sonar host 5 and the computer 6 for processing and analysis, and the multi-beam data acquisition software can display the data collected by the sonar and the angle data of the angle sensor module in real time.

[0087] The working principle of the multi-beam sounding system is: the core device of the multi-beam sounding system is a transceiver system and a transducer, and the transducer realizes the conversion between electrical signals and acoustic signals. As shown in Figure 7 The position of the transducer is known, the angle θ and the slant range r of the measured target and the transducer are measured, and the position (hi) of the target relative to the transducer can be calculated. At present, multi-beam is roughly divided into two categories, one is to measure the slant range r with the known transmission angle θ, and the other is to measure the slant range ri of the target to multiple transducers to solve the angle θ.

[0088] When the multi-beam works, the transducer transmits a fan-shaped acoustic wave in a set direction (usually vertically downward), the acoustic wave projects onto the water bottom to form a strip-shaped beam footprint, as shown in Figure 8 The acoustic wave is reflected from the water bottom and received by the transducer, and the relative three-dimensional coordinates of each beam footprint can be calculated through the time of the acoustic wave from the transducer to the water bottom, the speed of the acoustic wave and the angle of the acoustic wave transmission. The three-dimensional water depth data of the planar region is usually obtained through the movement of the transducer. In the multi-beam sounding system, as shown in Figure 7 The transmitting transducer and the receiving transducer both have directivity angles and form orthogonal fans, and the size of the beam footprint is determined by the product of the directivity angles of the two.

[0089] Unlike the traditional shipborne multi-beam underway measurement, the system rigidly connects the angle sensor and the sonar transducer on a vertical shaft through the bracket and the connecting rod, and realizes the three-dimensional scanning of the underwater terrain through synchronous rotation measurement by the rotating electromechanical system.

[0090] The above merely describes the preferred embodiments of the present application, and is not intended to limit the present application in any form. Although the present application has been disclosed with the preferred embodiments as above, it is not intended to limit the present application. Any person skilled in the art can make some changes or modifications to the above disclosed technical content to obtain equivalent embodiments with equivalent changes, as long as the changes or modifications do not deviate from the technical solution of the present application. Any modification, change, equivalent change and modification of the above embodiments made according to the technical essence of the present application still belong to the scope of the technical solution of the present application.

Claims

1. A method of matching a three-dimensional point cloud to a reconstruction model, the method comprising: Comprise the following steps: S1, normal estimation is carried out to three-dimensional point cloud;S2, according to the overall normal distribution situation, mark normal abnormal value;S3, near point iteration is carried out to normal abnormal value, and normal matching is carried out;S4, the point iteration of normal abnormal value is matched to normal point cloud position, realizes point cloud reconstruction; In the step S1, a normal estimation function is used to calculate a normal vector for each data point of the three-dimensional point cloud, points within a radius A of the data point are found, a coordinate of the points within the radius A is used to construct a covariance matrix, eigenvalue decomposition is performed on the covariance matrix to obtain three eigenvalues and corresponding eigenvectors, and an eigenvector corresponding to a minimum eigenvalue is a normal vector N of the data point i = (Nx i , Ny i , Nz i ); and the radius A is 3 times of an average distance of the point cloud. The step S2 comprises: S21, slicing the point cloud horizontally according to the Z coordinate of the point cloud; S22, calculating the geometric centroid point of each horizontal slice point cloud; S23, taking the slice geometric centroid point as the origin, constructing a horizontal slice plane coordinate system XOY, the X axis direction of the horizontal slice plane coordinate system is the same as the X axis direction of the original point cloud coordinate, the Y axis direction of the horizontal slice plane coordinate system is the same as the Y axis direction of the original point cloud coordinate, and calculating the slice coordinate system coordinates (x i , y i ) of each point; S24, calculating the normal difference , in the horizontal slice, and if the normal difference exceeds the threshold value 0.3, the normal is an abnormal value, and the point cloud of the normal abnormal value is a normal abnormal point cloud; In the step S3, the normal abnormal point cloud is regarded as source point cloud, and the normal point cloud is regarded as target point cloud, for each point in the normal abnormal point cloud, find its neighborhood point in the normal abnormal point cloud, and register it to the neighborhood point position of the normal point cloud, then adopt generalized iterative closest point registration and Levenberg-Marquardt error optimization fusion algorithm; The step S3 comprises: S31, initializing parameters, setting an initial transformation matrix T as an identity matrix, initially aligning the source point cloud and the target point cloud, setting a damping factor lambda initial value = 10 6 , setting a maximum iteration number = 100 times, and setting a convergence threshold = 10 -3 ; when a parameter update amount is less than the threshold, it is considered that the algorithm has converged; S32, matching the nearest neighbor points, for each point in the source point cloud, finding the nearest point in the target point cloud as a corresponding point; S33, minimizing the error, minimizing the error through an optimization method to obtain a new rigid body transformation matrix; and S34, applying the calculated rigid body transformation matrix T to the source point cloud to obtain a normal abnormal point cloud after normal point cloud registration. The step S33 comprises: S331、compute a residual vector r: compute the residual between the source point cloud and the corresponding target point cloud, each pair of corresponding points p i and q i , the residual r i may be defined as: ； ； where T is a transformation matrix that includes a rotation R and a translation matrix t, and the residual vector ; S332, error function is constructed: error function E is the square sum of residual vector, that is: ； S333, calculate Jacobian matrix J: calculate the partial derivative of residual to transformation parameter, construct Jacobian matrix J, each element of Jacobian matrix J represents the change rate of residual to each parameter, ； S334, update transformation matrix: use Levenberg-Marquardt method to iteratively optimize error function, find transformation matrix T that minimizes error function, in each iteration, calculate Jacobian matrix J and residual vector r, then calculate Gauss-Newton step , update transformation matrix T; S335、Adjust the damping factor: in each iteration, adjust the damping factor lambda value according to the change of the objective function value, if the objective function value decreases in the iteration, then: Make the algorithm closer to the Gauss-Newton method, if the objective function value increases, then: Reduce the step size to avoid excessive update; S336, check convergence condition: if the parameter update amount is less than the preset threshold, or the iteration number reaches the maximum value, the algorithm ends, that is, the target point cloud position is found; In the step S334, a Gauss-Newton step length is calculated: , a damping step length is calculated: , a transformation matrix is updated: , wherein J is a Jacobian matrix, I is an identity matrix, λ is an adjustment factor, is a parameter update amount, and r is a residual vector calculated in the previous step. Further comprising step S0, step S0 is located before step S1, and the step S0 is denoising to three-dimensional point cloud.

2. The method of claim 1, wherein: Further comprising step S5, step S5 is located after step S4, and the step S5 is inputting the reconstructed three-dimensional point cloud into matlab, to obtain three-dimensional perspective view and geometric dimension information.

Images