CN115760954A - Method for rapidly calculating surface area of complex surface based on point cloud data - Google Patents

Abstract

The present invention provides a method for quickly calculating the surface area of complex surfaces based on point cloud data, comprising: step 1. dividing the original point cloud of the complex surface to be calculated into a series of small facets, and projecting the points to the fitting surface; step 2. The projected point cloud F _porj is subjected to a specially designed rigid body transformation to obtain the point cloud F _tran ; Step 3. Extract the edge point F _bou in F _tran ; Step 4. Calculate the center point p _cen of F _bou , for any point in F _bou A point p _bou , calculate the normal vector determined by p _cen and p _bou

Then calculate the included angle θ _i , sort and store in θ _sort , and reorder the sorted edge points after interpolation to obtain a new set of edge points; step 5. For any two phases in the set of edge points The area of the triangle or sector surrounded by the center point is calculated for the adjacent points, and the surface area of the complex surface is obtained by accumulating all the patches.

Description

Method for fast calculation of complex surface surface area based on point cloud data

technical field

The invention belongs to the technical field of three-dimensional laser scanning, in particular to a method for quickly calculating the surface area of complex surfaces based on point cloud data.

Background technique

In practical work, it is often necessary to calculate the surface area of an object, and surface area calculation plays an important role in fields such as medical treatment, land survey, and agriculture. The surface area of a regular object can be realized by simple mathematical calculations, while the calculation of the surface area of an irregular object is more complicated, and the results obtained by different methods are not the same. In the past, the calculation of the surface area was mainly carried out by various measuring instruments such as GPS, total station, laser tracker, and laser rangefinder. According to the obtained discrete points and side lengths, surface fitting methods, curve fitting methods, and irregular Network methods and other methods to calculate the area. However, the accuracy of single-point measurement is generally limited by the accuracy of laser ranging and dial angle measurement, and the influence of ranging and angle measurement errors is difficult to weaken. These methods suffer from low efficiency, poor accuracy, and high computational cost, especially for inaccessible, large-scale, rough, and curved target objects.

Contents of the invention

The present invention is carried out to solve the problems of low efficiency, poor precision and high cost in the calculation of the surface area of large-scale complex surfaces. , large scale, rough and curved target objects.

In order to achieve the above object, the present invention adopts the following scheme:

The present invention provides a kind of method based on point cloud data fast calculation complex surface surface area, it is characterized in that, comprises the following steps: Step 1. adopts adaptive super voxel segmentation algorithm to divide the original point cloud P _orig of the complex surface to be calculated Form a series of small facets F _fac , fit each facet F _fac , and project the points on the facet to the corresponding accurately fitted plane;

Step 2. Let F _porj be the point cloud after F _fac projection, and perform a specially designed rigid body transformation on F _porj : take the center point on F _porj as the origin of the local coordinate system, and take the center point and the first point on F _porj The normal vector is the X-axis of the local coordinate system, the normal vector of F _porj is the Z-axis of the local coordinate system, and the cross product of the X-axis of the local coordinate system and the Z-axis of the local coordinate system is used as the Y-axis of the local coordinate system, and the rigid body transformation is performed to make F _{The porj} normal vector is parallel to the Z axis of the global coordinate system;

Step 3. Let the point cloud obtained by F _porj undergo rigid body transformation be F _tran . After rigid body transformation, the plane with any attitude angle in space will be parallel to the XY plane, and the z coordinate value of each point is the same, and then extract F The edge point F _bou in _tran ;

为起始方向，对于F_bou中的任 意一点p_bou，计算出由p_cen和p_bou确定的法向量

再计算出

与

之间的夹 角θ_i，将计算得到的角度θ_i按预定顺序(例如，顺时针或逆时针)排序并存储在 θ_sort中，对应点的顺序也随之调整并存储在P_sort中，针对已排好序的边缘点P_sort， 因有些平面的边缘点比较稀疏，为获得更加可靠的面积计算结果，若两个相邻边 缘点的距离大于一定阈值(如0.01m)，则对其进行插值加密，对插值后的点重新 排序，得到新的边缘点集合P_interp；Step 4. Calculate the center point p _cen of F _bou , let

is the starting direction, for any point p _bou in F _bou , calculate the normal vector determined by p _cen and p _bou

Then calculate

and

The angle θ _i between , the calculated angle θ _i is sorted in a predetermined order (for example, clockwise or counterclockwise) and stored in θ _sort , and the order of the corresponding points is also adjusted accordingly and stored in P _sort , For the sorted edge points P _sort , because the edge points of some planes are relatively sparse, in order to obtain more reliable area calculation results, if the distance between two adjacent edge points is greater than a certain threshold (such as 0.01m), the Perform interpolation encryption, reorder the interpolated points, and obtain a new set of edge points P _interp ;

Step 5. For any two adjacent points p ⁱ and p ^j in P _interp , the triangle or sector area S _area formed by p ⁱ , p ^j and p _cen is calculated, and all triangles or The sector areas are summed to obtain the surface area S _all of the complex surface.

This method is mainly to calculate the surface area on the unorganized point cloud, it has nothing to do with the source of the point cloud data, as long as the data required in our steps can be obtained, for example, this method is applicable to the point cloud obtained by oblique photography , or artificially synthesized simulation point clouds, are more suitable for LiDAR technology to obtain accurate three-dimensional information on the surface of objects without contact. For example, for some inaccessible, large-scale, rough, and curved target objects, it is difficult to obtain the accurate surface area by common measurement methods. First, use a 3D laser scanner to obtain the surface point cloud of the object, and then use this method to calculate the surface area .

Further, the method for quickly calculating the surface area of complex surfaces based on point cloud data provided by the present invention can also have the following characteristics: in the supervoxel segmentation process of step 1, the number of adjacent points (ε ₎ of each point is set to 20, If the proportion of planar structure features in P _orig reaches more than 50%, set the distance scale ε _d to 0.4-0.6m, otherwise set the distance scale ε _d to 0.3-0.2m.

Further, the method for quickly calculating the surface area of complex surfaces based on point cloud data provided by the present invention may also have the following features: In step 1, for each generated face F _fac , a robust plane based on PCA (Principal Component Analysis) is used The fitting method fits it, and the plane equation of the point in F _fac is set as:

Ax+By+Cz+D=0 (1)

In the formula, A, B and C are three normal vectors of the plane and satisfy A ² +B ² +C ² =1;

The distance d _i from point p _i on F _fac to the plane is:

d _i ＝|Ax _i +By _i +Cz _i +D| (2)

In the formula, x _i , y _i and z _i are the coordinate values of point p _i ;

应满足如下条件：objective function

The following conditions should be met:

In the formula, N is the number of points in F _fac ;

最小值，首先组成以下函数：Using Lagrange Multiplier Method to Solve the Objective Function

minimum, first compose the following function:

In the formula, λ is an undetermined parameter;

Calculate the partial derivative of formula (4) with respect to D and set the result of the derivative to 0 to obtain the following formula:

Then there are:

In the formula,

Formula (2) can be rewritten as the following form:

Using equation (4) to calculate the partial derivatives of A, B and C respectively, the following relationship can be obtained:

In the formula,

make

分别为M_3×3的特征向量，λ₁≥λ₂≥ λ₃分别为M_3×3的特征值，M_3×3分解为：Since M _3×3 is always positive semi-definite, let

are the eigenvectors of M _3×3 , and λ ₁ ≥ λ ₂ ≥ λ ₃ are the eigenvalues of M _3×3 respectively, and M _3×3 is decomposed into:

和F_fac的质 心点确定平面的距离t_i，则标准偏差σ为：In order to obtain more accurate fitting parameters, calculate each point p _i in F _fac to be determined by

and the centroid point of F _fac determine the distance t _i of the plane, then the standard deviation σ is:

和质心点； 设置迭代次数为ε₁，经过ε₁次迭代求出平面法向量A、B和C的精确值，然后利 用式(6)求出D的值；If t _i >2σ, remove the corresponding point from F _fac and use the remaining points to recalculate

and the centroid point; set the number of iterations to ε ₁ , and obtain the exact values of the plane normal vectors A, B and C after ε ₁ iterations, and then use formula (6) to obtain the value of D;

Project the points on F _fac onto the fitting plane, and the projection equation is:

In the formula, x _i , y _i and z _i are coordinate values before projection, and x _j , y _j and z _j are coordinate values after projection.

Furthermore, the method for quickly calculating the surface area of complex surfaces based on point cloud data provided by the present invention may also have the following features: In step 1, given the original point cloud data P _orig , use kd-tree to manage P _orig and establish its correspondence The index of P _orig is divided into a series of small facets using the adaptive supervoxel segmentation algorithm, which involves two parameters that need to be adjusted: the number of neighboring points ε _k for each point and the distance scale ε _d , and ε _d Can be regarded as the average side length of the patch, ε _k is usually set to 20. ε _d cannot be determined by a fixed value, but depends on the specific situation. If the original point cloud data mainly consists of planar structures, _εd should be set slightly larger (e.g. 0.5m), and if the point cloud contains rich surface features, _εd should be set slightly smaller (e.g. 0.25m).

Further, the method for quickly calculating the surface area of complex surfaces based on point cloud data provided by the present invention may also have the following characteristics: in step 2, the central point of F _porj is:

In the formula, x _porj ⁱ , y _porj ⁱ and z _porj ⁱ are the coordinate values of any point in F _porj , and N is the number of points in F _porj ;

Next, the local coordinate system on F _porj is defined as follows:

和

分别表示局部坐标系的三个坐标轴的方向，x₁，y₁和 z₁表示F_porj中第一个点的坐标值；In the formula,

and

represent the directions of the three coordinate axes of the local coordinate system, x ₁ , y ₁ and z ₁ represent the coordinate value of the first point in F _porj ;

The target coordinate axis system is defined as follows:

In the formula,

Let M _3×3 be the rotation matrix, then M _3×3 is calculated as follows:

Let p _i (x _i , y _i , zi ₎ be a point in F _porj , p _j (x _j , y _j , z _j ) be obtained after performing rigid body transformation on p _i (x _i , y _i , _{zi )} , then p _j (x _j ,y _j ,z _j ) is calculated as follows:

Further, the method for quickly calculating the surface area of complex surfaces based on point cloud data provided by the present invention may also have the following characteristics: in step 3, assume that the point cloud obtained after F _porj undergoes rigid body transformation is F _tran , after rigid body transformation, A plane with an arbitrary attitude angle in space will be parallel to the XY plane, which means that the z coordinate value of each point is the same, thus reducing the three-dimensional data to two-dimensional data. Next, use the two-dimensional α-shape algorithm to extract the edge points in F _tran (only x and y coordinate values participate in the operation). Generally speaking, the α-shape algorithm is a method of extracting edge points in a set of discrete points. Its principle can be regarded as a circle with a radius of ε ₂ rolling around the outside of F _tran , and the circle passes through any two points in F _tran , if there are no other points in the circle, these two points are regarded as edge points, and the edge points are obtained iteratively in this way. Let the obtained edge point be F _bou . It should be pointed out that since the average point density between different point cloud patches may not be the same, and sometimes the difference will be very large, so in order to better extract edge points, when setting the parameter ε ₂ (ε ₂ means The radius value of the ring when the α-shape algorithm is used to extract the edge points in F _tran ) is usually related to the average point spacing d _aver of F _tran , for example, set to 10 times the average point spacing of F _porj . For any point p _i in F _tran , perform k-nearest neighbor search to find the corresponding n neighboring points (expressed as p _i ¹ , p _i ² ..., p _i ⁿ ), then d _aver can be calculated as follows:

In the formula, N is the number of points of F _tran , 1≤i≤N. n can be set to 2.

Further, the method for quickly calculating the surface area of complex surfaces based on point cloud data provided by the present invention may also have the following characteristics: in step 4, the key to area calculation is to divide the polygon composed of F _bou into a series of small triangles or small fans , the total area can be obtained by counting the sum of the areas of these small triangles or small sectors; therefore, before calculating the area, it is necessary to regularize the edge points; first calculate the center point p _cen of F _bou :

In the formula, x _bou ⁱ and y _bou ⁱ are the coordinates of a point in F _bou respectively, and M is the number of points in F _bou ;

为起始方向，对于F_bou中的任意一点p_bou，由p_cen和p_bou确定的法 向量

计算如下：make

is the starting direction, for any point p _bou in F _bou , the normal vector determined by p _cen and p _bou

Calculated as follows:

In the formula, 1≤i≤M;

与

之间的夹角θ_i为：

and

The angle θ _i between is:

In the formula,

The calculated angle θ _i is sorted clockwise or counterclockwise and stored in θ _sort , 1≤i≤M, and the order of corresponding points is also adjusted accordingly and stored in P _sort .

Furthermore, the method for quickly calculating the surface area of complex surfaces based on point cloud data provided by the present invention may also have the following features: In step 4, for the sorted edge points P _sort , if the distance between two adjacent edge points If it is greater than a certain threshold (such as 0.01m), it needs to be interpolated and encrypted, and an improved 3rd degree B-spline curve is used for interpolation calculation, and the 3rd degree B-spline curve passes through 4 given control points P ₁ ~ _P4 ;

In the formula, P _i is the characteristic point of the control curve, P ₁ ~ P ₄ respectively correspond to the four points adjacent to the current edge point P ₀ , P ₀ = P _sort ; Fi(t) represents the third degree B-spline curve Basis function; t is an independent variable in F _i (t), t∈[0,1], that is, t takes any real number between 0 and 1;

For a sparse edge point, use the improved 3-time B-spline curve interpolation to encrypt it to ensure that the distance between every two interpolated points is not greater than the threshold, and reorder the interpolated points to get a new edge Point set P _interp .

Further, the method for quickly calculating the surface area of complex surfaces based on point cloud data provided by the present invention may also have the following features: In step 5, let p _interp ⁱ and p _interp ^j be two adjacent points in P _interp , if If more than 50% of the target surface to be calculated is composed of planar structures (the proportion of planar structure features is more than 50%), triangles are used for calculation. The triangular area S _area formed by p _interp ⁱ , p _interp ^j and p _cen is:

In the formula,

In the formula, (x _interp ⁱ , y _interp ⁱ ) and (x _interp ^j , y _interp ^j ) are the point coordinates of p _interp ⁱ and p _interp ^j respectively.

If more than 50% of the target surface to be calculated is composed of curved surfaces (surface structural features account for more than 50%), then a sector is used for calculation, and the sector _area S area formed by p _interp ⁱ , p _interp ^j and p _cen is calculated as follows:

in,

The surface area S _all of the final complex surface is calculated as follows:

为当前第m个面片F_fac中第n个三角形或扇形面积。In the formula, S is the total number of facets F _fac , R is the number of points in P _interp of the current mth facet F _fac ,

is the area of the nth triangle or sector in the current mth facet F _fac .

Function and Effect of Invention

The method for quickly calculating the surface area of complex surfaces based on point cloud data proposed by the present invention proposes for the first time to calculate the surface area directly on the unorganized point cloud without the need for 3D reconstruction, and can be applied to a variety of complex surfaces. The invention uses point cloud data to automatically and quickly calculate the surface area of complex surfaces, aiming at the problems of low efficiency, poor accuracy and high calculation cost in most existing methods, especially when calculating the surface area of inaccessible, large-scale, rough and curved targets , the advantage of this method is particularly obvious. It not only greatly improves the work efficiency, but also has high reliability by comparing with the measurement data collected manually by reliable equipment.

Description of drawings

Fig. 1 is the flow chart of the method for calculating complex surface surface area quickly based on point cloud data that the embodiment of the present invention relates to;

Fig. 2 is a schematic diagram of super-voxel segmentation involved in an embodiment of the present invention, wherein (a) is the original point cloud data, and (b) is the result of super-voxel segmentation;

Fig. 3 is a comparison diagram before and after projection related to an embodiment of the present invention, wherein (a) is before projection, and (b) is after projection;

4 is a schematic diagram of rigid body transformation involved in an embodiment of the present invention;

FIG. 5 is a schematic diagram of a point cloud patch after rigid body transformation according to an embodiment of the present invention;

Fig. 6 is a schematic diagram of edge point extraction related to an embodiment of the present invention, wherein (a) is before extraction, and (b) is after extraction;

Fig. 7 is a schematic diagram of the plane projection involved in the embodiment of the present invention, wherein (a) is before projection, and (b) is after projection;

Fig. 8 is the improved 3 B-spline curve interpolation that the embodiment of the present invention relates to the schematic diagram of edge point encryption effect, wherein (a) is before encryption, (b) is after encryption;

Fig. 9 is a schematic diagram of the area calculation results involved in the embodiment of the present invention;

Fig. 10 is the simulation point cloud data related to the embodiment of the present invention, wherein (a) is a cube, (b) is a sphere, and (c) is a ring;

Fig. 11 is the calculation result of the method proposed by the present invention related to the embodiment of the present invention under different sampling rates, wherein (a) is calculation time, and (b) is relative accuracy;

Fig. 12 is the calculation result of the greedy triangle algorithm involved in the embodiment of the present invention at different sampling rates, wherein (a) is the calculation time, and (b) is the relative accuracy;

Fig. 13 is the calculation result of the Poisson surface reconstruction algorithm involved in the embodiment of the present invention at different sampling rates, where (a) is the calculation time, and (b) is the relative accuracy;

Figure 14 is the calculation results of the three methods involved in the embodiment of the present invention under different Gaussian noises, where (a) is the method proposed by the present invention, (b) is the greedy triangle algorithm, and (c) is the Poisson surface reconstruction algorithm ;

Fig. 15 is the actual experiment 1 that the embodiment of the present invention involves, wherein (a) is " 3S " sculpture, (b) is target area 1, (c) is the stadium of Hubei Institute of Science and Technology, (d) is target area 2;

Fig. 16 is the real experiment 2 involved in the embodiment of the present invention, wherein (a1) is the original geometric model of Horse, (b1) is the point cloud data of Horse, (c1) is the supervoxel segmentation result of Horse; (a2) is The original geometric model of Skeleton Hand, (b2) is the point cloud data of Skeleton Hand, (c2) is the supervoxel segmentation result of Skeleton Hand.

Detailed ways

The specific embodiment of the method for calculating complex surface surface area based on point cloud data that the present invention relates to is described in detail below in conjunction with accompanying drawing.

<Example>

As shown in Figure 1, the method for calculating complex surface surface area quickly based on point cloud data provided by the present embodiment comprises the following steps:

Step 1. In order to demonstrate the calculation process of the present invention's method conveniently, the cube simulation point cloud as shown in Figure 2 is selected as an example here, and a simulation point cloud P _orig is obtained by sampling six faces of a cube whose side length is 1m , each face contains 26,896 points, for a total of 161,376 points. Supervoxel segmentation involves two parameters: the number of nearest neighbors ε _k for each point and the distance scale ε _d , _which can be regarded as the average edge length of a patch. In our experiments, ε _k is set to 20. _εd cannot be determined by a fixed value, but depends on the specific situation. If the original point cloud data contains rich (more than 50%) planar features, _εd should be set slightly larger (for example, 0.5m), if the point cloud If the data contains abundant (more than 50%) surface features, then _εd should be set slightly smaller (eg 0.25 m). The proposed method involves two parameters: ε ₁ and ε ₂ , in this embodiment, ε ₁ is set to 200, and ε ₂ is set to 10 times the average point spacing of F _porj . Unless otherwise stated, the above parameter settings were used for all experiments involved in this example without adjustment. After super-voxel segmentation, a total of 28 facets F _fac are generated. For each generated surface F _fac , a robust plane fitting method based on PCA (Principal Component Analysis) is used for fitting, the number of iterations is set to 200, and all points on F _fac are projected onto the fitting surface. As shown in Figure 3, after projection, the small facet with noise becomes very flat, which is very beneficial to the extraction of edge points and area calculation.

Step 2. As shown in Figure 4, perform rigid body transformation on the projected patch F _porj so that its normal vector is parallel to the Z axis (that is, the plane after the rigid body transformation is parallel to the XY plane), and then the point after the rigid body transformation is obtained Cloud noodles F _bou .

Step 3. As shown in Figure 5, after the rigid body transformation, the plane with any attitude angle in space will be parallel to the X-Y plane, which means that the z coordinate values of each point are the same, thus reducing the three-dimensional data to two dimension data. Next, use the improved adaptive α-shape algorithm to extract the edge point F _bou in F _tran (only the x and y coordinate values participate in the operation), as shown in Figure 6, the edge point of the patch is successfully extracted.

Step 4. As shown in Figure 7, take the north direction as the positive direction, sort the edge points clockwise, calculate the Euclidean distance between every two points after sorting, if it is greater than 0.01m, use the improved 3 times B sample The curves are interpolated so that the distance between any two adjacent points is not greater than 0.01m, as shown in Figure 8, which is the interpolation of sparse edge points.

Step 5. As shown in FIG. 9 , calculate the area S _area of each patch, and summarize to obtain the total area S _all .

As shown in Figure 10, in order to further illustrate the reliability of this method, in addition to the cube, two other simulated point clouds were selected for testing. The radius of the sphere Sphere is 1m, the total number of points is 160,000, and ε _d is set to 0.25m. The inner diameter of Torus is 0.3m, the outer diameter is 1m, the total number of points is 160,000, and ε _d is set to 0.5m. Point clouds are generally very dense. In order to improve efficiency, it is usually necessary to downsample the point cloud before segmentation. In order to test the influence of different sampling rates on the calculation results, three simulation point clouds were down-sampled at five different density levels. Table 1 shows the statistical results of the surface area calculation method proposed by the present invention on three simulation point clouds, from which the following conclusions can be drawn: the proposed surface area calculation method is very efficient and the result is reliable, in planar structures and curved surfaces (such as spheres And the ring) can get the calculation results which are highly consistent with the theoretical value.

In order to test the influence of different sampling rates on the calculation results, three simulation point clouds were down-sampled at five different density levels. As shown in Fig. 11(a) and Fig. 11(b), proper sampling does not necessarily reduce the relative accuracy compared to performing surface area calculation directly on the original point cloud, and the computational efficiency is greatly improved with the increase of sampling rate, And can still maintain high reliability and precision.

Table 1

In order to further verify the performance of the method, the greedy triangle algorithm and the Poisson surface reconstruction algorithm are used to reconstruct the sampled point cloud (these two algorithms are prior art methods), and then the off-the-shelf software CloudCompare is used to calculate the value of each mesh grid. The surface area of the mesh. As shown in Fig. 12, when the surface is smooth and the point cloud is evenly distributed, the greedy triangle algorithm is not sensitive to the sampling rate, in this case, the area result calculated by the greedy triangle algorithm is highly reliable. The Poisson surface reconstruction method fits and approximates all points to generate a closed surface that is watertight and has good geometrical characteristics. Although Poisson surface reconstruction is a good overall reconstruction method, it is only suitable for smooth closed surfaces without sharp features. For the sake of fairness, only the reconstruction results of the first two scenarios are used for comparison experiments. As shown in Figure 13(a) and Figure 13(b), although the calculated surface area based on Poisson surface reconstruction has high accuracy for specific objects with smooth surfaces, this method is very time-consuming, and the method proposed in the present invention It is an order of magnitude (more than 10 times) faster than methods based on Poisson surface reconstruction. In fact, affected by various factors, the surface of point cloud is usually filled with a lot of random noise. In order to test the robustness of the proposed method, the area calculation is next performed from the noisy point cloud, as shown in Fig. 14(a), the proposed method can perform area well even in the presence of a large The calculation shows that the method has good noise immunity. Since the greedy triangle algorithm is very sensitive to noise, only a few small Gaussian noise values are set here. However, with the increase of noise, the calculation accuracy decreases rapidly (see Figure 14(b)). As shown in Figure 14(c), the Poisson surface reconstruction algorithm is not sensitive to Gaussian noise, but this algorithm is not suitable for all types of surfaces, to be precise, it is only suitable for some specific surfaces. The results compared with the prior art show that the method of the invention not only greatly improves the calculation speed, but also has better noise resistance and calculation accuracy, and is applicable to all types of surfaces.

As shown in Figure 15(a), in December 2020, we scanned a sculpture named "3S" with a Faro Focus ^s 150 laser scanner at Wuhan University, with a sampling rate of 1/20. The sculpture is composed of flat and curved parts , the middle part is a quadrangular prism. In order to facilitate comparison and avoid errors caused by point cloud registration, only the point cloud data collected from one scanning angle is used here, as shown in Figure 15(b), only the flat part in the middle is kept, and the number of points is 130,308. Next, the precise coordinates of the six corners of the target area are acquired using a LeicaAT960 laser tracker. The point measurement accuracy of the laser tracker is 15±6 microns, and the measured area is 6.9218m ² . Note that the parameter _εd is set to 0.5m, using triangles to calculate the area. As shown in Figure 15(c), a low-altitude UAV based on Changzhou Xinyi X650 flying platform and Regel Mini300 laser scanner was used to collect the 3D color point cloud of Hubei University of Science and Technology, and a relatively flat football field was selected as the target area ( See Figure 15(d)), the number of target area points is 2,974,086. Then use the Topcon total station to measure the coordinates of the four corners of the football field, and the measured area is 7154.722m ² . Note that because the point cloud of the target area is relatively sparse, the parameter _εd is set to 5m, and the sector is used to calculate the area. The area calculation results of the two target areas are listed in Table 2. It can be seen that the method proposed in the present invention has successfully calculated the area of the test point cloud data within an acceptable time, which is very close to the measured value, greatly improving work efficiency.

Table 2

As shown in Figure 16, in order to test the processing effect of the method proposed by the present invention on complex surfaces, two geometric models named Horse and Skeleton Hand are selected for detailed display. Both models were produced from Georgia Tech's large archive of geometric models, the first model was enlarged 10 times to approximate true scale, while the second model was left unchanged. The reference value is calculated by CloudCompare. The parameter ε _d is set to 0.2m, and the sector is used to calculate the area. As shown in Table 3, even for very complex curved surfaces, the method proposed in the present invention can still achieve good results.

table 3

The above embodiments are merely illustrations for the technical solution of the present invention. The method for quickly calculating the surface area of complex surfaces based on point cloud data involved in the present invention is not limited to the content described in the above embodiments, but is subject to the scope defined in the claims. Any modification or supplement or equivalent replacement made by those skilled in the art of the present invention on the basis of the embodiment is within the scope of protection required by the claims of the present invention.

Claims (7)

1. A method for rapidly calculating the surface area of a complex surface based on point cloud data is characterized by comprising the following steps:

step 1, adopting a self-adaptive hyper-voxel segmentation algorithm to carry out point cloud P on an original point of a complex surface to be calculated _orig Divided into a series of facets F _fac For each patch F _fac Fitting, and projecting the points on the surface patch to a corresponding accurately-fitted plane;

step 2, setting F _porj Is F _fac Projected point cloud, pair F _porj Performing specially designed rigid body transformation: with F _porj The central point on is the origin of the local coordinate system, and the central point and F are used _porj The normal vector of the first point is the X axis of the local coordinate system, and F is used _porj The normal vector of (A) is the Z axis of the local coordinate system, the cross product of the X axis of the local coordinate system and the Z axis of the local coordinate system is used as the Y axis of the local coordinate system, rigid body transformation is executed to enable F to be F _porj The normal vector is parallel to the Z axis of the whole coordinate system;

step 3, setting F _porj After rigid body transformation, the product is obtainedThe point cloud is F _tran After rigid body transformation, a plane with any attitude angle in space is parallel to the X-Y plane, the z coordinate value of each point is the same, and F is extracted _tran Edge point F in (1) _bou ；

Step 4, calculating F _bou Central point p of _cen Let us order

For the starting direction, for F _bou At any point p in _bou Calculate by p _cen And p _bou Determined normal vector

Then calculate out

And

angle theta therebetween _i Will calculate the obtained theta _i Sorting in a predetermined order and storing at theta _sort In the method, the sequence of corresponding points is adjusted and stored in P _sort In the method, the edge points P are arranged in order _sort If the distance between two adjacent edge points is greater than a certain threshold value, the two adjacent edge points are subjected to interpolation encryption, and points after interpolation are reordered to obtain a new edge point set P _interp ；

Step 5 for P _interp Any two adjacent points p in (2) ⁱ And p ^j All calculate by p ⁱ 、p ^j And p _cen Formed triangular or sector area S _area Adding all triangular or sector areas in all the panels to obtain the surface area S of the complex surface _all 。

2. The method for rapidly calculating the surface area of a complex surface based on point cloud data of claim 1, wherein:

wherein, in step 1In the process of superpixel segmentation, the number epsilon of adjacent points of each point _k Set to 20 if P _orig If the ratio of the structural features of the middle plane reaches more than 50%, the distance scale epsilon _d Set to 0.4-0.6 m, otherwise, the distance scale epsilon _d Set to 0.3 to 0.2m.

3. The method for rapidly calculating the surface area of a complex surface based on point cloud data of claim 1, wherein:

wherein, in step 1, for each patch F generated _fac Fitting the stable plane by adopting a stable plane fitting method based on PCA (principal component analysis), and setting F _fac The equation of the plane where the point is located is:

Ax+By+Cz+D＝0 (1)

wherein A, B and C are three normal vectors of a plane and satisfy A ² +B ² +C ² ＝1；

F _fac Last point p _i Distance d to the plane _i Comprises the following steps:

d _i ＝|Ax _i +By _i +Cz _i +D| (2)

in the formula, x _i ，y _i And z _i Is a point p _i The coordinate values of (a);

objective function

The following conditions should be satisfied:

in the formula, N is F _fac The number of points in;

solving an objective function using a lagrange multiplier method

The minimum value, first of all, constitutes the following function:

in the formula, lambda is a parameter to be determined;

and (5) calculating the partial derivative of D by the formula (4) and making the derivative result be 0 to obtain the following formula:

then there are:

in the formula,

formula (2) may be rewritten as follows:

by using equation (4) to calculate the partial derivatives for A, B and C, the following relationship can be obtained:

in the formula,

order to

Due to M _3×3 Is always symmetrically and semi-positively determined, such that

Are each M _3×3 Of the feature vector λ ₁ ≥λ ₂ ≥λ ₃ Are respectively M _3×3 Characteristic value of (1), M _3×3 The decomposition is as follows:

calculate F _fac Each point p in _i To by

And F _fac Determines the distance t of the plane _i Then the standard deviation σ is:

if t is _i If > 2 sigma, the corresponding point is selected from F _fac Removing, recalculating with the remaining points

And a centroid point; setting the number of iterations to epsilon ₁ Through epsilon ₁ The precise values of the plane normal vectors A, B and C are solved through the secondary iteration, and then the value of D is solved through the formula (6);

f is to be _fac The points are projected onto a fitting plane, and the projection equation is as follows:

in the formula, x _i 、y _i And z _i Is a coordinate value, x, before projection _j 、y _j And z _j Are the coordinate values after projection.

4. The method for rapidly calculating the surface area of a complex surface based on point cloud data of claim 1, wherein:

wherein, in step 2, F _porj The central point of (a) is:

in the formula, x _porj ⁱ ，y _porj ⁱ And z _porj ⁱ Is F _porj The coordinate value of any point in the coordinate system, N is F _porj The number of points in;

next, F _porj The local coordinate system above is defined as follows:

in the formula,

and

respectively representing the directions of three coordinate axes of a local coordinate system, x ₁ ，y ₁ And z ₁ Is represented by F _porj Coordinate values of the first point;

the target coordinate axis system is defined as follows:

in the formula,

let M _3×3 Is a rotation matrix, then M _3×3 The calculation is as follows:

let p _i (x _i ,y _i ,z _i ) Is F _porj Middle point, p _j (x _j ,y _j ,z _j ) Is to p _i (x _i ,y _i ,z _i ) Point obtained after rigid body transformation, then p _j (x _j ,y _j ,z _j ) The calculation is as follows:

5. the method for rapidly calculating the surface area of a complex surface based on point cloud data of claim 1, wherein:

in step 4, F is first calculated _bou Central point p of _cen ：

In the formula, x _bou ⁱ And y _bou ⁱ Are respectively F _bou The coordinate value of a middle point, M is F _bou The number of points in;

order to

For the starting direction, for F _bou At any point p in _bou From p _cen And p _bou Determined normal vector

The calculation is as follows:

in the formula, i is more than or equal to 1 and less than or equal to M;

and with

Angle theta therebetween _i Comprises the following steps:

in the formula,

the calculated angle theta _ii Ordered clockwise or counterclockwise and stored at theta _sort In which i is more than or equal to 1 and less than or equal to M, the order of the corresponding points is also adjusted and stored in P _sort In (1).

6. The method for rapidly calculating the surface area of a complex surface based on point cloud data of claim 1, wherein:

wherein, in step 4, the edge points P are ordered _sort Interpolation is carried out using a modified 3-degree B-spline and the 3-degree B-spline is passed through 4 given control points P ₁ ～P ₄ ：

In the formula, P _i For controlling characteristic points of the curve, P ₁ ～P ₄ Respectively correspond to the current edge point P ₀ Four adjacent points, P ₀ ＝P _sort (ii) a Fi (t) represents the basis function of the 3-th-order B-spline curve; t is F _i (t) taking any real number between 0 and 1 as an independent variable;

for a sparse edge point, encrypting the sparse edge point by using improved 3-time B-spline curve interpolation to ensure that the distance between every two interpolation points is not greater than a threshold value, reordering the interpolated points to obtain a new edge point set P _interp 。

7. The method of claim 1 for fast computation of complex surface area based on point cloud data, wherein:

wherein, in step 5, let p _interp ⁱ And p _interp ^j Is P _interp If the ratio of the structural features of the target surface plane to be calculated reaches more than 50%, calculating by adopting a triangle, and calculating by using p _interp ⁱ 、p _interp ^j And p _cen The area S of the triangle _area Comprises the following steps:

in the formula,

in the formula (x) _interp ⁱ ,y _interp ⁱ ) And (x) _interp ^j ,y _interp ^j ) Are each p _interp ⁱ And p _interp ^j Point coordinates of (a);

if the ratio of the structural features of the curved surface of the target surface to be calculated exceeds 50%, calculating by adopting a fan shape, and calculating by using p _interp ⁱ 、p _interp ^j And p _cen Formed sector area S _area The calculation is as follows:

wherein,

surface area S of the resulting complex surface _all The calculation is as follows:

wherein S is a patch F _fac R is the current mth panel F _fac P of _interp The number of points in (a) is,

is the current mth face sheet F _fac The nth triangle or sector area.

Images