# Fast Cylindrical Fitting Method Using Point Cloud’s Normals Estimation

Authors: Yunlong Wu; Qing Zhang; Shuxuan Zhang
Journal: Mathematical Problems in Engineering
 (2018)

---

## 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.


---

## Body

## 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 Figure1.Figure 1
Neighborhood point cloud search schematic. For any data point in a cylinder, a certain number of neighborhood points are selected for normals estimation.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 C is(1)C=1k∑i=1kpi-p¯·pi-p¯T(2)C·vt→=λt·vt→,t∈0,1,2In formula (1), k is the number of nearest neighbor points of the selected point P. p¯ is the three-dimensional centroid of a nearest neighbor point of point P, λt is the t eigenvalue of the covariance matrix, and vt→ is the t eigenvector corresponding to the values of the first eigenvalue. When selecting the k 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 k.The eigenvectors and eigenvalues of the covariance matrix are analyzed, and the eigenvectorv0→, which corresponded to the minimum eigenvalue λ0 of the covariance matrix C 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(3)x-x0a=y-y0b=z-z01Thez axis direction of cylindrical axis is fixed to 1. At this time, the axis direction vector of the cylinder is w→=a0,b0,1; the vector is perpendicular to the point cloud surface normal ni→=xi,yi,zi. Set function:(4)f=∑i=1mxi·a0+yi·b0+zi2In this formula, m is the total number of points in the point cloud.The order functionf is derived from the parameters a0 and b0 and then sets zero. The initial values a0 and b0 of a and b are obtained as follows: (5)
a
0
b
0
=
-
∑
i
=
1
m
x
i
2
∑
i
=
1
m
x
i
·
y
i
∑
i
=
1
m
x
i
·
y
i
∑
i
=
1
m
y
i
2
-
1
·
∑
i
=
1
m
x
i
·
z
i
∑
i
=
1
m
y
i
·
z
iBased on the dual quaternion theory, Plücker coordinates are generally composed of six parametersL,M,N,L0,M0,N0, and the vector β=L,M,N and γ=L0,M0,N0 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 x1,y1,z1, x2,y2,z2, then there is(6)L,M,N=x2-x1,y2-y1,z2-z1L0,M0,N0=x1,y1,z1×x2,y2,z2=y1z2-y2z1,z1x2-z2x1,x1y2-x2y1In 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.(7)l=L+εL0i+M+εM0j+N+εN0kIn this formula,ε is a dual unit and satisfies ε2=0, ε≠0 and i, j, and k are the imaginary units and satisfy i2=j2=k2=-1.Plücker straight linel1 is transformed into Plücker line through space transformation. The spiral displacement equation of l1′ is(8)l1′=q^·l1·q^-1In this formula,q^ is a unit dual of four elements and q^-1 is a conjugate of q^.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 thez axis in the O-xyz coordinate system, the initial value of the axis direction vector is w→=a,b,1, and the coordinate conversion matrix is(9)R=1-u21+w-u·v1+w-u-u-u·w1+w1-v21+wuvwIn this formula, (10)u=a0a02+b02+1;v=b0a02+b02+1;w=1a02+b02+1.The cylindrical axis of the above coordinate transformation is approximately parallel to thez axis, and the x and y coordinates of all the surface points of the cylinder are equal to the coordinates of the projection to the plane z=0, so the x position and the y position after the coordinate transformation can be used for a round fitting. Assuming that the center coordinate is x0′,y0′ after the round fitting, the coordinates of a point on the axis can be(11)x0y0z0=R-1·x0′y0′z0′In order to reduce the influence of the fitting error of the cylindrical axis, the average value of thez coordinates of all points after the coordinate transformation is selected as the z coordinates of the axis; that is, in this paper, the z coordinate of the points after the coordinate transformation is selected as(12)zo′=∑i=1mzi′m


## 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 C is(1)C=1k∑i=1kpi-p¯·pi-p¯T(2)C·vt→=λt·vt→,t∈0,1,2In formula (1), k is the number of nearest neighbor points of the selected point P. p¯ is the three-dimensional centroid of a nearest neighbor point of point P, λt is the t eigenvalue of the covariance matrix, and vt→ is the t eigenvector corresponding to the values of the first eigenvalue. When selecting the k 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 k.The eigenvectors and eigenvalues of the covariance matrix are analyzed, and the eigenvectorv0→, which corresponded to the minimum eigenvalue λ0 of the covariance matrix C 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(3)x-x0a=y-y0b=z-z01Thez axis direction of cylindrical axis is fixed to 1. At this time, the axis direction vector of the cylinder is w→=a0,b0,1; the vector is perpendicular to the point cloud surface normal ni→=xi,yi,zi. Set function:(4)f=∑i=1mxi·a0+yi·b0+zi2In this formula, m is the total number of points in the point cloud.The order functionf is derived from the parameters a0 and b0 and then sets zero. The initial values a0 and b0 of a and b are obtained as follows: (5)
a
0
b
0
=
-
∑
i
=
1
m
x
i
2
∑
i
=
1
m
x
i
·
y
i
∑
i
=
1
m
x
i
·
y
i
∑
i
=
1
m
y
i
2
-
1
·
∑
i
=
1
m
x
i
·
z
i
∑
i
=
1
m
y
i
·
z
iBased on the dual quaternion theory, Plücker coordinates are generally composed of six parametersL,M,N,L0,M0,N0, and the vector β=L,M,N and γ=L0,M0,N0 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 x1,y1,z1, x2,y2,z2, then there is(6)L,M,N=x2-x1,y2-y1,z2-z1L0,M0,N0=x1,y1,z1×x2,y2,z2=y1z2-y2z1,z1x2-z2x1,x1y2-x2y1In 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.(7)l=L+εL0i+M+εM0j+N+εN0kIn this formula,ε is a dual unit and satisfies ε2=0, ε≠0 and i, j, and k are the imaginary units and satisfy i2=j2=k2=-1.Plücker straight linel1 is transformed into Plücker line through space transformation. The spiral displacement equation of l1′ is(8)l1′=q^·l1·q^-1In this formula,q^ is a unit dual of four elements and q^-1 is a conjugate of q^.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 thez axis in the O-xyz coordinate system, the initial value of the axis direction vector is w→=a,b,1, and the coordinate conversion matrix is(9)R=1-u21+w-u·v1+w-u-u-u·w1+w1-v21+wuvwIn this formula, (10)u=a0a02+b02+1;v=b0a02+b02+1;w=1a02+b02+1.The cylindrical axis of the above coordinate transformation is approximately parallel to thez axis, and the x and y coordinates of all the surface points of the cylinder are equal to the coordinates of the projection to the plane z=0, so the x position and the y position after the coordinate transformation can be used for a round fitting. Assuming that the center coordinate is x0′,y0′ after the round fitting, the coordinates of a point on the axis can be(11)x0y0z0=R-1·x0′y0′z0′In order to reduce the influence of the fitting error of the cylindrical axis, the average value of thez coordinates of all points after the coordinate transformation is selected as the z coordinates of the axis; that is, in this paper, the z coordinate of the points after the coordinate transformation is selected as(12)zo′=∑i=1mzi′m


## 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 z coordinate value of a point on the fixed cylinder axis is z0, and the optimized objective parameter is X=a,b,x0,y0, and the optimized objective function is(13)ϕX=∑i=1mφi2X=∑i=1mLi2+Mi2+Ni2a2+b2+1-r02In this formula,(14)Li=ayi-y0-bxi-x0;Mi=azi-z0-xi-x0;Ni=bzi-z0-yi-y0.In formula (13), the value of r0 is recalculated after each iteration, and the formula is(15)r0=∑i=1mdimIn this formula, di is the distance between the point i and the axis [21].


### 3.2. Iterative Incremental Calculation

The incremental formula of Levenberg-Marquardt algorithm is [22](16)dp=JkTJk+λkI-1JkTεkIn formula (16), λk is the damping factor. The initial damping factor λ0 is 0.01, I is the m-order unit matrix, and Jk is the Jacobi matrix, and the formula is(17)
J
k
=
∂
φ
1
X
∂
a
∂
φ
1
X
∂
b
∂
φ
1
X
∂
x
0
∂
φ
1
X
∂
y
0
∂
φ
2
X
∂
a
∂
φ
2
X
∂
b
∂
φ
2
X
∂
x
0
∂
φ
2
X
∂
y
0
⋮
∂
φ
m
X
∂
a
∂
φ
m
X
∂
b
∂
φ
m
X
∂
x
0
∂
φ
m
X
∂
y
0
m
×
4εk is the computation error after each iteration, and its value is(18)εk=-φ1kX,-φ2kX,⋯,-φmkXTThe termination condition is JkTεk≤ε, and in this paper, the threshold ε is 10-5.


## 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 z coordinate value of a point on the fixed cylinder axis is z0, and the optimized objective parameter is X=a,b,x0,y0, and the optimized objective function is(13)ϕX=∑i=1mφi2X=∑i=1mLi2+Mi2+Ni2a2+b2+1-r02In this formula,(14)Li=ayi-y0-bxi-x0;Mi=azi-z0-xi-x0;Ni=bzi-z0-yi-y0.In formula (13), the value of r0 is recalculated after each iteration, and the formula is(15)r0=∑i=1mdimIn this formula, di is the distance between the point i and the axis [21].


## 3.2. Iterative Incremental Calculation

The incremental formula of Levenberg-Marquardt algorithm is [22](16)dp=JkTJk+λkI-1JkTεkIn formula (16), λk is the damping factor. The initial damping factor λ0 is 0.01, I is the m-order unit matrix, and Jk is the Jacobi matrix, and the formula is(17)
J
k
=
∂
φ
1
X
∂
a
∂
φ
1
X
∂
b
∂
φ
1
X
∂
x
0
∂
φ
1
X
∂
y
0
∂
φ
2
X
∂
a
∂
φ
2
X
∂
b
∂
φ
2
X
∂
x
0
∂
φ
2
X
∂
y
0
⋮
∂
φ
m
X
∂
a
∂
φ
m
X
∂
b
∂
φ
m
X
∂
x
0
∂
φ
m
X
∂
y
0
m
×
4εk is the computation error after each iteration, and its value is(18)εk=-φ1kX,-φ2kX,⋯,-φmkXTThe termination condition is JkTεk≤ε, and in this paper, the threshold ε is 10-5.


## 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 Figure2. Figure 3 is the model after 3D transformation.Figure 2
Model before 3D transformation.Figure 3
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. Figure4 shows the pipeline point cloud obtained from the extraction process of the target pipelines from the 3D scan data.Figure 4
Fast acquisition of pipelines. This group of point cloud data is scanned by 3D laser scanner.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. Figure5 is a complete pipeline point cloud model after splicing.Figure 5
Complete point cloud of pipeline. This point cloud data is obtained by using the fast extraction technology of the above model.Figure6 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.Figure 6
Evaluation of flatness. 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.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 Figure7. 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.Figure 7
Evaluation interface of flatness. Through the evaluation interface of flatness, the flatness evaluation result of pipeline model can be obtained directly.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 pointsk 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 k in this paper is 20, the initial values of the cylinder fitting can reach a high precision. So k =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 10-5 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.Table 1
Data fitting results of the first group of pipeline point cloud; unit: m.
Method of obtaining initial value Unit vector of cylindrical axis direction A point on the axis of a cylinder Radius (r) Number of iterations u v w This paper method 0.3639 0.3414 0.4681 (-3.0371,-2.0910,0.5211) 0.6074 2  Literature [13] method 0.3639 0.3414 0.4681 (-3.0332,-2.0886,0.5251) 0.6068 81Table 2
Data fitting results of the second group of pipeline point cloud; unit: m.
Method of obtaining initial value Unit vector of cylindrical axis direction A point on the axis of a cylinder Radius (r) Number of iterations u v w  This paper method 0.3748 0.3307 0.4877 (-5.3132,-3.1070,-0.4741) 0.6078 2  Literature [13] method 0.3748 0.3307 0.4878 (-5.3121,-3.1072,-0.4729) 0.6070 71Table 3
The distance statistics of the point cloud to the fitting cylinder; unit: m.
Fitting point cloud Maximum value of positive distance Maximum value of negative distance Average distance The first set of point cloud 0.0046 -0.0018 4.0643 × 1 0 - 4 The second set of point cloud 0.0019 -0.0020 4.2539 × 1 0 - 4The values of latitude and longitude density ofm1 and m2 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 m1 and m2 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 Section4.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.


## 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 Figure2. Figure 3 is the model after 3D transformation.Figure 2
Model before 3D transformation.Figure 3
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. Figure4 shows the pipeline point cloud obtained from the extraction process of the target pipelines from the 3D scan data.Figure 4
Fast acquisition of pipelines. This group of point cloud data is scanned by 3D laser scanner.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. Figure5 is a complete pipeline point cloud model after splicing.Figure 5
Complete point cloud of pipeline. This point cloud data is obtained by using the fast extraction technology of the above model.Figure6 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.Figure 6
Evaluation of flatness. 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.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 Figure7. 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.Figure 7
Evaluation interface of flatness. Through the evaluation interface of flatness, the flatness evaluation result of pipeline model can be obtained directly.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 pointsk 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 k in this paper is 20, the initial values of the cylinder fitting can reach a high precision. So k =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 10-5 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.Table 1
Data fitting results of the first group of pipeline point cloud; unit: m.
Method of obtaining initial value Unit vector of cylindrical axis direction A point on the axis of a cylinder Radius (r) Number of iterations u v w This paper method 0.3639 0.3414 0.4681 (-3.0371,-2.0910,0.5211) 0.6074 2  Literature [13] method 0.3639 0.3414 0.4681 (-3.0332,-2.0886,0.5251) 0.6068 81Table 2
Data fitting results of the second group of pipeline point cloud; unit: m.
Method of obtaining initial value Unit vector of cylindrical axis direction A point on the axis of a cylinder Radius (r) Number of iterations u v w  This paper method 0.3748 0.3307 0.4877 (-5.3132,-3.1070,-0.4741) 0.6078 2  Literature [13] method 0.3748 0.3307 0.4878 (-5.3121,-3.1072,-0.4729) 0.6070 71Table 3
The distance statistics of the point cloud to the fitting cylinder; unit: m.
Fitting point cloud Maximum value of positive distance Maximum value of negative distance Average distance The first set of point cloud 0.0046 -0.0018 4.0643 × 1 0 - 4 The second set of point cloud 0.0019 -0.0020 4.2539 × 1 0 - 4The values of latitude and longitude density ofm1 and m2 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 m1 and m2 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 Section4.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.



---

*Source: 8904653-2018-12-12.xml*