Abstract
Cylindrical fitting is an essential step in Large Process Pipeline’s measurement process, and precision of initial values of cylindrical fitting is a key element in getting a correct fitting result. In order to get well initial values, covariance matrixes of all points in cylinder’s three-dimensional laser scanning point cloud should be firstly established to estimate normals of all points, and then cylinder’s axis vector can be calculated by using least squares method. Secondly, remaining parameters’ initial values of the cylinder can be got by coordinate transformation. Finally, Levenberg-Marquardt algorithm is used in iterative optimization process to get fitting result by using the above values as initial values. Experiments demonstrate that this method can get precise initial values of cylindrical fitting and improve the accuracy and speed of cylindrical fitting.
1. Introduction
With the application of 3D laser scanner in the fields of reverse engineering, industrial manufacturing, cultural relic protection and medical visualization, people have been able to obtain all kinds of original point cloud data at a lower cost as a basis for further processing, such as point based rendering, point based shape modeling, and surface weight, construction, and so on. A necessary attribute in point based representation is normal vector information. In addition to the precise and high quality point based rendering methods mainly relying on the normal vector [1, 2], there are many surface reconstruction algorithms that also need to get accurate reconstruction results with the help of accurate normal vector [3–6].
Large-scale process pipeline is the foundation of pipeline transportation industry, and the fitting of cylindrical parameters is a very necessary link in its construction monitoring process. At present, as a new measuring method, three-dimensional laser scanning is used to obtain the whole digital model [7] by measuring the point cloud data of the object. It can provide accurate measurement data for the monitoring of the large process pipeline. Practice proves that the selection of initial parameters of cylindrical fitting will greatly affect the convergence, speed, and precision of fitting. At present, the method of 3D point cloud vector estimation mainly includes the neighborhood search and fuzzy estimation [8, 9], and literature [10] analyses the normal estimation of the sharp feature scattered points. Literature [11] compares the point cloud data in the size detection field of the industrial plant pipeline system with the 3D CAD data. Literature [12] improves the error equation of the indirect adjustment with constraints and selects any initial values to get the correct fitting results, but the speed of adjustment is slow. Literature [13] accurately calculates the initial values by projecting the cylindrical point to the plane and analyzing the degree of the projection point closed to the circle, but the quantity of initial values is easily influenced by the projection distance; at the same time, the calculation speed of the adjustment is also slow. In literature [14], the initial values of the axis direction vector are calculated by the geometric relation between the cylindrical surface point and the axis direction vector of the cylindrical axis. In literature [15], the initial values of the axis direction vector are calculated by Gauss diagram. The essence of the two is to fit the two surfaces of the cylindrical surface according to the surface of the cylinder. Good initial values can be obtained, but the computation is large.
In this regard, a cylindrical fitting method is proposed to estimate the initial values of cylindrical parameters by estimating normal points of cylindrical dots. For all data points in the cylindrical point cloud, a certain number of neighborhood points are selected to estimate the normal direction of the point, as shown in Figure 1.
After estimating the normal line of the scanned point, the initial values of the axis direction vector are calculated by the least square method according to the vertical relation between the normal direction of the point and the axis direction of the cylinder, and then the initial values of a point and the cylindrical radius on the axis of the cylinder are calculated by coordinate transformation and the round least square fitting. After obtaining the initial values of all cylindrical parameters, the number of parameters, according to the constraint relation of cylindrical parameters, is reduced to replace the adjustment with constraint conditions, and the Levenberg-Marquardt algorithm is used to optimize the parameters of cylinder iteratively.
2. Initial Values Calculation Based on Point Cloud Normals Estimation
2.1. Normals Estimation Based on Near Neighbor
The normals estimation method based on the nearest neighbor point directly infers the surface normals according to the coordinates of data points in the point cloud. That is, the covariance matrix of the point is established in three-dimensional point cloud, and the normal direction [16] of the point is estimated approximately by covariance analysis. For any point P in the point cloud, the corresponding covariance matrix is
In formula (1), is the number of nearest neighbor points of the selected point P. is the three-dimensional centroid of a nearest neighbor point of point P, is the eigenvalue of the covariance matrix, and is the eigenvector corresponding to the values of the first eigenvalue. When selecting the nearest neighbor points of point P, if all the three-dimensional coordinates are traversed, the amount of calculation is very large, so the k-d tree [17] is used to divide the point cloud in space to realize the fast matching of the nearest neighbor of the .
The eigenvectors and eigenvalues of the covariance matrix are analyzed, and the eigenvector , which corresponded to the minimum eigenvalue of the covariance matrix is solved, that is, the normal vector direction of the point P. The obtained normal vector has two meanings; that is, it cannot determine whether it is positive or negative of the point cloud estimation method vector, but because the vector initial values of the axis direction are solved by the vertical relation with the surface normal vector of the point cloud, it is estimated that the two senses of the point cloud normals vector do not affect the estimation of the initial values of the cylindrical axis.
2.2. The Estimation of the Initial Values of the Axis Based on the Least Square Method
After the estimation of the point cloud surface normals, the least square method is used to calculate the initial value of the axis direction vector according to the surface normal vector of the point cloud, and the axis of the cylinder is expressed as
The axis direction of cylindrical axis is fixed to 1. At this time, the axis direction vector of the cylinder is ; the vector is perpendicular to the point cloud surface normal . Set function:In this formula, is the total number of points in the point cloud.
The order function is derived from the parameters and and then sets zero. The initial values and of and are obtained as follows:
Based on the dual quaternion theory, Plücker coordinates are generally composed of six parameters , and the vector and are the direction vector of the line and the moment of the line relative to the origin, respectively [18]. The direction vector of a line describes the direction of a line. The space position of the line is described by the moment relative to the origin. If Plücker passes straight point , , then there is
In order to describe the coordinate transformation of a line in three-dimensional space, a Plücker line is generally expressed as a unit dual vector.
In this formula, is a dual unit and satisfies , and , , and are the imaginary units and satisfy .
Plücker straight line is transformed into Plücker line through space transformation. The spiral displacement equation of is
In this formula, is a unit dual of four elements and is a conjugate of .
Therefore, Plücker coordinates are mainly used to describe the spiral motion of space vector, which is more stable for continuous spiral motion. In this paper, the least square method is used to estimate the axis of the cylinder, and all the points obtained by scanning are calculated. The spatial expression of the axis of the cylinder is obtained. This method is more accurate and robust than Plücker coordinates.
2.3. Solution of Initial Values of Residual Parameters
In order to obtain the coordinates of a point on the axis and the initial values of the radius of the cylinder, the coordinate transformation of the cylindrical point cloud is carried out, the axis direction is parallel to the axis in the coordinate system, the initial value of the axis direction vector is , and the coordinate conversion matrix isIn this formula,
The cylindrical axis of the above coordinate transformation is approximately parallel to the axis, and the and coordinates of all the surface points of the cylinder are equal to the coordinates of the projection to the plane , so the position and the position after the coordinate transformation can be used for a round fitting. Assuming that the center coordinate is after the round fitting, the coordinates of a point on the axis can be
In order to reduce the influence of the fitting error of the cylindrical axis, the average value of the coordinates of all points after the coordinate transformation is selected as the coordinates of the axis; that is, in this paper, the coordinate of the points after the coordinate transformation is selected as
3. Optimization of Cylinder Parameters Based on Levenberg-Marquardt Algorithm
3.1. Optimization Model of Cylinder Fitting
This paper reduces the number of parameters to be optimized according to the constraint relation of cylindrical parameters in order to avoid the establishment of optimal functions with constraints, and the parameter optimization of cylindrical fitting is carried out by using the Levenberg-Marquardt optimization algorithm [19] which is based on the initial value of the cylindrical parameters obtained. This algorithm can achieve the advantages of Gauss-Newton algorithm and gradient descent algorithm by modifying the parameters during execution and improve the shortcomings of both (for example, the inverse matrix of Gauss-Newton algorithm does not exist or the initial value is too far from the local minimum) [20]. The advantage of Levenberg-Marquardt method is that it can be adjusted: if the descent is too fast, use a smaller to make it closer to Gauss-Newton method; if the descent is too slow, use a larger to make it closer to the gradient descent method. The equation of cylindrical axis is shown in formula (3) in the second section. The coordinate value of a point on the fixed cylinder axis is , and the optimized objective parameter is , and the optimized objective function is
In this formula,
In formula (13), the value of is recalculated after each iteration, and the formula isIn this formula, is the distance between the point and the axis [21].
3.2. Iterative Incremental Calculation
The incremental formula of Levenberg-Marquardt algorithm is [22]
In formula (16), is the damping factor. The initial damping factor is 0.01, is the -order unit matrix, and is the Jacobi matrix, and the formula is is the computation error after each iteration, and its value isThe termination condition is , and in this paper, the threshold is .
4. Experimental Data Verification
4.1. Experimental Data Processing
The pipeline is scanned by LEICA HDS7000 3D laser scanner to obtain the point cloud data of the pipeline and establish the feature 3D model with the large pipeline in the pipeline storage and processing site as the experimental object. The model before 3D transformation is shown in Figure 2. Figure 3 is the model after 3D transformation.
Due to the small deviation between the position of the model and the actual position of the point cloud, the size of the 3D model can be properly amplified to ensure that the model contain all the pipeline point cloud data. At the same time, considering the distance between the pipeline and the adjacent pipelines, the size of the 3D feature model should not exceed 1.2 times. When the size of the pipeline is increased by 1.1 to 1.2 times, the generated feature model is fitted to the location of the target model which is interested in the measured data, and the data points outside the feature model are deleted at one time; thus the target model is obtained. Figure 4 shows the pipeline point cloud obtained from the extraction process of the target pipelines from the 3D scan data.
In order to build a complete data model, the point cloud data, obtained from all the viewpoints, must be integrated into a unified coordinate system and spliced into a whole. According to the point cloud stitching algorithm, the point cloud data can be quickly and accurately stitching. Figure 5 is a complete pipeline point cloud model after splicing.
Figure 6 is the result of flatness evaluation at any position on the pipeline model. The blue line is the optional feature line obtained by discretization and then fitting; the black mark position is beyond the standard and to be adjusted. Because it is necessary to rotate many times in order to ensure the accuracy of rotation in the process of model rotation, we set 2 mm as the standard. For the flatness evaluation results, if the difference is less than 2 mm, then it meets the standard; if it is more than 2 mm, it does not meet the standard. Meanwhile it will mark this area through black marking points.
At the same time, the MATLAB GUI interface is designed, and the fitting evaluation process and results are displayed intuitively in the interface. The interface display dialog box is shown in Figure 7. The MATLAB fitting and interface program is transformed into executable file. Through the secondary development of cloud compare, the group pair evaluation function is integrated into the point cloud data processing platform. The whole pipeline measurement data processing process is formed into a system.
Due to the problem of pipeline storage, a complete cylinder model cannot be observed, so some effective points of the pipeline are calculated. First, the pipeline point cloud is denoised and the outliers are removed. At the same time, in order to reduce the computation of cylindrical fitting, it is carried out after the point cloud data is sampled sparsely according to absolute spacing [23].
4.2. Cylinder Fitting Result
In order to ensure the accuracy of the normal estimation, the value of the number of neighborhood points in formula (1) should be guaranteed enough for the normals estimation of the cylindrical point cloud obtained by the scanning. After the initial calculation, when the value of in this paper is 20, the initial values of the cylinder fitting can reach a high precision. So =20 is selected as the number of nearest neighbor points of the point cloud normals estimation.
In order to verify the speediness and effectiveness of the cylindrical fitting method based on the point cloud linear estimation, the method of obtaining the initial values of the cylindrical parameters optimization in literature [13] is compared with the method of obtaining the initial parameters of the column in this paper, and the Levenberg-Marquardt algorithm is used in the two methods. Cylindrical fitting optimization is carried out for the extracted two-point cloud data. The termination threshold of iteration is set to according to the above section. The designed standard radius of the pipeline is 0.607m, and the allowable manufacturing error is (±2.5mm); that is, the actual diameter of the scanned pipeline is in the range of 0.6045-0.6095m. This value is used as the reference value of the fitting result. The fitting results are shown in Tables 1 and 2. Table 3 shows the maximum positive and negative distances from the point cloud to the fitting cylinder of the two groups of pipelines fitted by this method.
The values of latitude and longitude density of and of the hemispherical vector set in [13] method are 90 and 360, respectively. The first group of pipeline point cloud get the initial values, and then the iterative termination threshold is satisfied after 81 iterations, while the second group of pipeline point cloud has more iterations, so the and values of the second set of point cloud are changed to 240 and 960, and then, after 71 iterations, it satisfies the iteration termination threshold. However, after the initial values of the cylindrical parameters are obtained by normal estimation, the termination thresholds can be satisfied by only two iterations.
Through the fitting results, we can draw a conclusion that the initial values can be obtained by two methods, and the accurate cylindrical parameters can be obtained by using the Levenberg-Marquardt algorithm for iterative optimization. But compared with literature [13] method, the method of cylindrical point cloud estimation can be faster for getting the threshold values which greatly improves the speed of cylindrical fitting.
4.3. Analysis of the Efficiency of Cylinder Fitting
In order to verify the overall efficiency of the method of cylindrical fitting in this paper, MATLAB is used to time the initial values operation of Section 4.2 cylinder fitting. The time of obtaining initial values of the first group of data in the document [13] is 65.705s, while the time of obtaining initial values by this method is 13.640s, and the time of obtaining initial values by the second group of data in the document [13] is 63.637s, while the time of obtaining initial values by this method is 13.408s. It is obvious that the initial value of this method is faster than that of [13]. It can be found that for a large number of point cloud data, the method of this paper can obtain the first value of the cylinder fitting at a fast speed to the requirement of the threshold, and the number of iterations is less, thus reducing the time of the whole fitting of the cylinder, and the overall efficiency is high.
5. Conclusion
The key of the cylinder fitting is the initial values of the fitting. In this paper, the method of cylindrical point cloud method is used to obtain the initial values of the cylinder fitting. It can solve the very accurate cylindrical initial values by less calculation and less iterations. It can effectively improve the velocity of the cylinder fitting and greatly improve the overall calculation efficiency of the cylinder fitting. It is very suitable for engineering foundation. Cylindrical fitting of cylindrical 3D laser scanning data also plays an important role in the process of measuring related geometric parameters of large-scale process pipelines. It also has some reference significance for cylindrical fitting in other fields.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by Qingdao Institute of Marine Engineering, Tianjin University.
Supplementary Materials
The is the main program of the MATLAB, and all the calculations and drawings are done with it. The is the pipe mouth data, and the is the pipe wall data. The , , and are three measurements of point cloud data. The is the front-end point cloud data, the is the first set of pipeline point cloud data, and the is the second set of pipeline point cloud data. The is MATLAB code for measuring the radius and direction vector of the cylinder and the number of iterations. The is the MATLAB code for measuring the time used in this paper and the is the MATLAB code for measuring the time used in the [13]. (Supplementary Materials)