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

Abstract

The application provides a method for quickly calculating the surface area of a complex surface based on point cloud data, comprising the following steps: 1. segmenting the original point cloud of the complex surface to be calculated into a series of small facets, and projecting the points to a fitting surface; 2. performing a rigid body transformation specially designed on the projected point cloud F porj to obtain a point cloud F tran ; 3. extracting edge points F tran in F bou ; 4. calculating the center point p bou of F cen , and for any point p bou in F bou , calculating a normal vector determined by p cen and p bou , then calculating an included angle θ i , sorting and storing θ sort , reordering the edge points after interpolation, and obtaining a new edge point set; and 5. for any two adjacent points in the edge point set, calculating the triangular or sector area surrounded by the center point, and accumulating all the facet areas to obtain the surface area of the complex surface.

Description

Technical Field

[0001] This invention belongs to the field of three-dimensional laser scanning technology, specifically relating to a method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data. Background Technology

[0002] In practical work, it is often necessary to calculate the surface area of ​​objects. Surface area calculation plays an important role in fields such as medicine, land surveys, and agriculture. The surface area of ​​regular objects can be calculated using simple mathematical calculations, while the calculation of the surface area of ​​irregular objects is more complex, and different methods yield different results. Previously, the calculation of land surface area mainly used various measuring instruments such as GPS, total stations, laser trackers, and laser rangefinders. Based on the obtained discrete points and side lengths, methods such as surface fitting, curve fitting, and irregular network methods were used to calculate the area. However, the accuracy of single-point measurements is generally limited by the accuracy of laser ranging and graduated angle measurement, and the influence of ranging and angle measurement errors is difficult to mitigate. These methods are inefficient, have poor accuracy, and are computationally expensive, especially for inaccessible, large-scale, rough, and curved targets. Summary of the Invention

[0003] This invention addresses the problems of low efficiency, poor accuracy, and high cost in calculating the surface area of ​​large-scale complex surfaces. Its purpose is to provide a method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data, which is particularly suitable for target objects that are difficult to access, large in scale, rough, and curved.

[0004] To achieve the above objectives, the present invention employs the following solution:

[0005] This invention provides a method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data, characterized by the following steps: Step 1. Using an adaptive supervoxel segmentation algorithm to segment the original point cloud P of the complex surface to be calculated... orig Divided into a series of small pieces F fac For each facet F fac Perform a fit and project the points on the patch onto the corresponding precisely fitted plane;

[0006] Step 2. Let F porj For F fac The projected point cloud, for F porj Perform specially designed rigid body transformations: using F porj The center point on the graph is the origin of the local coordinate system, and the coordinates of the center point and F are... porj The normal vector of the first point is the local coordinate system X-axis, with F porj The normal vector is the Z-axis of the local coordinate system, and the cross product of the X-axis and Z-axis of the local coordinate system is the Y-axis of the local coordinate system. A rigid body transformation is performed to make F...porj The normal vector is parallel to the Z-axis of the global coordinate system;

[0007] Step 3. Let F porj The point cloud obtained after rigid body transformation is F tran After rigid body transformation, a plane with arbitrary orientation angles in space will be parallel to the XY plane, and the z-coordinate value of each point will be the same. Then, F is extracted. tran edge point F in bou ;

[0008] Step 4. Calculate F bou center point p cen ,make As the starting direction, for F bou Any point p in bou Calculate the result from p cen and p bou A defined normal vector Then calculate and The included angle θ i The calculated angle θ i Sort and store in a predetermined order (e.g., clockwise or counterclockwise) at θ sort In the process, the order of the corresponding points is also adjusted and stored in P. sort In the middle, 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), interpolation is performed to refine the points, and the interpolated points are reordered to obtain a new set of edge points P. interp ;

[0009] Step 5. For P interp Any two adjacent points p in i and p j All were calculated from p i p j and p cen The area S of the triangle or sector formed area The surface area S of the complex surface is obtained by adding up the areas of all triangles or sectors within all the facets. all .

[0010] This method primarily calculates surface area on unorganized point clouds, regardless of the source of the point cloud data. As long as the data required in our steps is available, it will suffice. For example, this method is suitable for point clouds obtained using oblique photography or artificially synthesized simulated point clouds, and is particularly suitable for non-contact acquisition of precise 3D information about object surfaces using LiDAR technology. For instance, for some inaccessible, large-scale, rough, or curved target objects, ordinary measurement methods may struggle to obtain accurate surface areas. In such cases, a 3D laser scanner can be used to first obtain a point cloud of the object's surface, and then this method can be used to calculate the surface area.

[0011] Furthermore, the method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data provided by this invention may also have the following characteristics: during the supervoxel segmentation process in step 1, the number ε of neighboring points of each point is... k Set to 20, if P orig If the proportion of planar structural features reaches more than 50%, then the distance scale ε will be... d Set it to 0.4–0.6 m, otherwise change the distance scale ε. d Set to 0.3~0.2m.

[0012] Furthermore, the method for rapidly 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 1, for each generated surface patch F fac A robust plane fitting method based on PCA (Principal Component Analysis) is used to fit it, and F is set to... fac The equation of the plane containing the point in the middle is:

[0013] Ax + By + Cz + D = 0 (1)

[0014] In the formula, A, B, and C are three normal vectors of the plane and satisfy A 2 +B 2 +C 2 =1;

[0015] F fac A little bit above p i Distance d to the plane i for:

[0016] d i =|Ax i +By i +Cz i +D| (2)

[0017] In the formula, x i y i and z i For point p i The coordinate values;

[0018] objective function The following conditions must be met:

[0019]

[0020] In the formula, N is F fac The number of points in the middle;

[0021] Solving the objective function using the Lagrange multiplier method The minimum value is first determined by the following function:

[0022]

[0023] In the formula, λ is a parameter to be determined;

[0024] Taking the partial derivative of equation (4) with respect to D and setting the result to 0, we obtain the following equation:

[0025]

[0026] Then we have:

[0027]

[0028] In the formula,

[0029] Equation (2) can be rewritten as follows:

[0030]

[0031] Using equation (4) with respect to A, B, and C respectively, we can obtain the following relationship:

[0032]

[0033] In the formula,

[0034] make

[0035]

[0036] Because M 3×3 It is always symmetric positive semi-definite, let M respectively 3×3 The eigenvectors, λ1≥λ2≥λ3, are M respectively. 3×3 Eigenvalues, M 3×3 Decomposed into:

[0037]

[0038] To obtain more accurate fitting parameters, F is calculated. fac Each point p ini From and F fac The distance t of the plane is determined by the centroid. i Then the standard deviation σ is:

[0039]

[0040] If t i If the value is greater than 2σ, then the corresponding point will be moved from F. fac Remove from the middle, and recalculate using the remaining points. And the centroid; set the number of iterations to ε1, and after ε1 iterations, find the exact values ​​of the plane normal vectors A, B and C, and then use equation (6) to find the value of D;

[0041] F fac The projection equation of a point on the plane onto the fitted plane is:

[0042]

[0043] In the formula, x i y i and z i The coordinates before projection, x j y j and z j These are the projected coordinate values.

[0044] Furthermore, the method for rapidly 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 1, given the original point cloud data P orig Use kd-tree to manage P orig And establish its corresponding index, and use the adaptive supervoxel segmentation algorithm to segment P orig The process involves dividing the surface into a series of smaller faces, which requires adjusting two parameters: the number of neighboring points ε for each point. k and distance scale ε d , and ε d It can be considered as the average side length of the face, ε k It is usually set to 20. ε d It cannot be determined by a fixed value, but rather depends on the specific circumstances. If the original point cloud data mainly consists of planar structures, then ε should be... d Set it to a slightly larger value (e.g., 0.5m). If the point cloud contains rich surface features, then ε should be set higher. d Set it to a slightly smaller value (e.g., 0.25m).

[0045] Furthermore, the method for rapidly 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, F porj The center point is:

[0046]

[0047] In the formula, x porj i y porj i and z porj i For F porj The coordinates of any point in N, where N is F porj The number of points in the middle;

[0048] Next, F porj The local coordinate system on is defined as follows:

[0049]

[0050] In the formula, and These represent the directions of the three coordinate axes of the local coordinate system, x1, y1, and z1, respectively, and represent F. porj The coordinates of the first point in the middle;

[0051] The target coordinate system is defined as follows:

[0052]

[0053] In the formula,

[0054]

[0055]

[0056] Let M 3×3 If M is a rotation matrix, then 3×3 The calculation is as follows:

[0057]

[0058] Let p i (x i ,y i ,z i ) is F porj At the middle point, p j (x j ,y j ,z j ) for p i (x i ,y i ,z i The point obtained after performing a rigid body transformation, then p j (x j ,y j ,z j The calculation is as follows:

[0059]

[0060] Furthermore, the method for rapidly 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, assuming F porj The point cloud obtained after rigid body transformation is F tran After rigid body transformation, a plane with arbitrary orientation angles in space will be parallel to the XY plane. This means that the z-coordinate value of each point is the same, thus reducing the three-dimensional data to two-dimensional data. Next, the two-dimensional α-shape algorithm is used to extract F. tran Edge points in a discrete set (only x and y coordinates are used in the calculation). Generally, the α-shape algorithm is a method for extracting edge points from a discrete set of points. Its principle can be viewed as a circle with radius ε² surrounding F. tran The ring rolls outside, passing through F. tran For any two points in the circle, if no other points are in the circle, these two points are considered edge points. Edge points are obtained iteratively in this way. Let the obtained edge points be F. bou It should be noted that since the average point density between different point cloud patches may not be the same, and sometimes the difference can be significant, in order to better extract edge points, the ε2 parameter should be set (ε2 represents the value of the α-shape algorithm for extracting F-shape points). tran The radius of the annulus is usually set to the value of F at the edge point of the annulus. tran Average point spacing d aver Related, for example, set to F porj Ten times the average point spacing. For F tran Any point p in i Perform a k-nearest neighbor search to find the corresponding n nearest neighbors (denoted as p). i 1 ,p i 2 …,p i n ), then d aver The following can be calculated:

[0061]

[0062] In the formula, N is F tran The number of points, 1 ≤ i ≤ N. n can be set to 2.

[0063] Furthermore, the method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data provided by this invention may also have the following characteristics: In step 4, the key to area calculation is to divide the area calculated by F... bouThe resulting polygon is divided into a series of smaller triangles or sectors. The total area can be obtained by summing the areas of these smaller triangles or sectors. Therefore, before calculating the area, the edge points need to be regularized. First, calculate F. bou center point p cen :

[0064]

[0065] In the formula, x bou i and y bou i F respectively bou The coordinates of a point in the middle, M is F bou The number of points in the middle;

[0066] make As the starting direction, for F bou Any point p in bou , by p cen and p bou A defined normal vector The calculation is as follows:

[0067]

[0068] In the formula, 1≤i≤M;

[0069] and The included angle θ i for:

[0070]

[0071] In the formula,

[0072]

[0073] The calculated angle θ i Sort by clockwise or counterclockwise and store in θ sort In the case where 1≤i≤M, the order of the corresponding points is adjusted accordingly and stored in P. sort middle.

[0074] Furthermore, the method for rapidly 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, for the sorted edge points P sort If the distance between two adjacent edge points is greater than a certain threshold (e.g., 0.01m), then interpolation densification is required. An improved cubic B-spline curve is used for interpolation calculation, and the cubic B-spline curve passes through four given control points P1 to P4.

[0075]

[0076]

[0077]

[0078]

[0079]

[0080] In the formula, P i As feature points controlling the curve, P1 to P4 correspond to the four points adjacent to the current edge point P0, respectively, where P0 = P sort Fi(t) represents the basis function of a cubic B-spline curve; t is the basis function of F. i The independent variable in (t) is t∈[0,1], that is, t takes any real number between 0 and 1;

[0081] For a sparse edge point, an improved cubic B-spline interpolation is used to densify it, ensuring that the distance between any two interpolated points is no greater than a threshold. The interpolated points are then reordered to obtain a new set P of edge points. interp .

[0082] Furthermore, the method for rapidly 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 5, let p interp i and p interp j It is P interp If more than 50% of the target surface to be calculated consists of planar structures (planar structural features account for more than 50%), then triangles are used for calculation between two adjacent points in the p matrix. interp i p interp j and p cen The area S of the triangle formed area for:

[0083]

[0084] In the formula,

[0085]

[0086] In the formula, (x interp i ,y interp i ) and (x interp j ,y interp j ) are p respectivelyinterp i and p interp j The coordinates of the point.

[0087] If more than 50% of the target surface to be calculated is composed of curved surfaces (curved surface structural features account for more than 50%), then a sector is used for calculation, determined by p. interp i p interp j and p cen The area S of the sector formed area The calculation is as follows:

[0088]

[0089] in,

[0090]

[0091] The final surface area S of the complex surface all The calculation is as follows:

[0092]

[0093] In the formula, S represents the surface F. fac The total number, R is the current m-th face F fac P interp The number of points in the middle, For the current m-th face F fac The area of ​​the nth triangle or sector in the equation.

[0094] The role and effect of invention

[0095] The method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data proposed in this invention is the first to propose performing surface area calculation directly on unorganized point clouds without the need for 3D reconstruction. This method is applicable to a wide variety of complex surfaces. By automatically and rapidly calculating the surface area of ​​complex surfaces using point cloud data, this invention addresses the problems of low efficiency, poor accuracy, and high computational cost inherent in most existing methods. Its advantages are particularly evident when calculating the surface area of ​​inaccessible, large-scale, rough, and curved targets. It not only significantly improves work efficiency but also demonstrates high reliability when compared with measurement data collected manually using reliable equipment. Attached Figure Description

[0096] Figure 1 This is a flowchart of a method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data, as described in an embodiment of the present invention.

[0097] Figure 2This is a schematic diagram of supervoxel segmentation according to an embodiment of the present invention, where (a) is the original point cloud data and (b) is the supervoxel segmentation result;

[0098] Figure 3 These are before-and-after comparison images of the projection involved in the embodiments of the present invention, wherein (a) is before projection and (b) is after projection;

[0099] Figure 4 This is a schematic diagram of rigid body transformation involved in an embodiment of the present invention;

[0100] Figure 5 This is a schematic diagram of a point cloud patch after rigid body transformation according to an embodiment of the present invention;

[0101] Figure 6 This is a schematic diagram of edge point extraction according to an embodiment of the present invention, where (a) is before extraction and (b) is after extraction;

[0102] Figure 7 This is a schematic diagram of a planar projection involved in an embodiment of the present invention, wherein (a) is before projection and (b) is after projection;

[0103] Figure 8 This is a schematic diagram illustrating the edge point encryption effect of the improved cubic B-spline curve interpolation according to an embodiment of the present invention, where (a) is before encryption and (b) is after encryption.

[0104] Figure 9 This is a schematic diagram of the area calculation results involved in an embodiment of the present invention;

[0105] Figure 10 The above refers to the simulated point cloud data involved in the embodiments of the present invention, wherein (a) is a cube, (b) is a sphere, and (c) is a ring;

[0106] Figure 11 The above are the calculation results of the method proposed by the present invention under different sampling rates in the embodiments of the present invention, where (a) is the calculation time and (b) is the relative accuracy;

[0107] Figure 12 The results of the greedy triangle algorithm according to the embodiments of the present invention at different sampling rates are shown, where (a) is the calculation time and (b) is the relative accuracy.

[0108] Figure 13 The results of the Poisson surface reconstruction algorithm according to the embodiments of the present invention at different sampling rates are shown, where (a) is the calculation time and (b) is the relative accuracy.

[0109] Figure 14 The following are the calculation results of three methods involved in the embodiments of the present invention under different Gaussian noise, wherein (a) is the method proposed in the present invention, (b) is the greedy triangle algorithm, and (c) is the Poisson surface reconstruction algorithm;

[0110] Figure 15 This is a real experiment 1 involved in the embodiment of the present invention, wherein (a) is a “3S” sculpture, (b) is target area 1, (c) is the stadium of Hubei University of Science and Technology, and (d) is target area 2;

[0111] Figure 16 This is a real experiment 2 involved in the embodiments of the present invention, where (a1) is the original geometric model of Horse, (b1) is the point cloud data of Horse, and (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, and (c2) is the supervoxel segmentation result of Skeleton Hand. Detailed Implementation

[0112] The following describes in detail, with reference to the accompanying drawings, the specific implementation scheme of the method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data involved in this invention.

[0113] <Example>

[0114] like Figure 1 As shown in the figure, the method for quickly calculating the surface area of ​​complex surfaces based on point cloud data provided in this embodiment includes the following steps:

[0115] Step 1. To facilitate the demonstration of the calculation process of the method of the present invention, we first select as follows: Figure 2 The example shown is a simulated point cloud of a cube. A simulated point cloud P is obtained by sampling the six faces of a cube with a side length of 1m. orig Each face contains 26,896 points, for a total of 161,376 points. Supervoxel segmentation involves two parameters: the number of nearest neighbors ε for each point. k and distance scale ε d , ε d It can be considered as the average side length of the facet. In our experiment, ε k Set to 20. ε d It 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, then ε should be set to... d Set it to a slightly larger value (e.g., 0.5m). If the point cloud data contains abundant (more than 50%) surface features, then ε should be set to a larger value. d The value is set to a slightly smaller value (e.g., 0.25m). The proposed method involves two parameters: ε1 and ε2. In this embodiment, ε1 is set to 200 and ε2 is set to F. porjTen times the average pixel spacing. Unless otherwise stated, the above parameter settings will be used for all experiments involved in this embodiment without adjustment. After hypervoxel segmentation, a total of 28 facets F were generated. fac For each generated facet F fac All fittings were performed using a robust plane fitting method based on PCA (Principal Component Analysis), with the number of iterations set to 200. fac All points on the surface are projected onto the fitting surface, such as Figure 3 As shown, after projection, the noisy small facets become very flat, which is very beneficial for edge point extraction and area calculation.

[0116] Step 2. As Figure 4 As shown, for the projected surface F porj Perform a rigid body transformation so that its normal vector is parallel to the Z-axis (i.e., the transformed plane is parallel to the XY plane), and then obtain the point cloud patch F after the rigid body transformation. bou .

[0117] Step 3. As Figure 5 As shown, after rigid body transformation, a plane with arbitrary orientation angles in space will be parallel to the XY plane, meaning that the z-coordinate value of each point is the same, thus reducing the three-dimensional data to two-dimensional data. Next, an improved adaptive α-shape algorithm is used to extract F... tran edge point F in bou (Only x and y coordinate values ​​are involved in the calculation), such as Figure 6 As shown, the edge points of the facet were successfully extracted.

[0118] Step 4. Figure 7 As shown, taking the north direction as the positive direction, the edge points are sorted clockwise. The Euclidean distance between any two points after sorting is calculated. If it is greater than 0.01m, an improved cubic B-spline curve is used for interpolation to ensure that the distance between any two adjacent points is no greater than 0.01m. Figure 8 As shown, this is an inset of sparse edge points.

[0119] Step 5. Figure 9 As shown, calculate the area S of each facet. area The total area S is obtained by summing the results. all .

[0120] like Figure 10 As shown, to further illustrate the reliability of this method, in addition to the cube, two other simulated point clouds were selected for testing. The sphere has a radius of 1m and a total of 160,000 points, ε... dSet to 0.25m, Torus inner diameter 0.3m, outer diameter 1m, total number of points 160,000, ε d The value was set to 0.5m. Point clouds are generally very dense, and to improve efficiency, they usually need to be downsampled before segmentation. To test the impact of different sampling rates on the calculation results, three simulated point clouds were downsampled at five different density levels. Table 1 shows the statistical results of the surface area calculation method proposed in this invention on the three simulated point clouds. From this, we can draw the following conclusions: the proposed surface area calculation method is very efficient and reliable, and can obtain calculation results that are in high agreement with theoretical values ​​on both planar structures and curved surfaces (such as spheres and tori).

[0121] To test the impact of different sampling rates on the calculation results, three simulated point clouds were downsampled at five different density levels. For example... Figure 11 (a) and Figure 11 As shown in (b), compared with the surface area calculation directly on the original point cloud, appropriate sampling does not necessarily reduce the relative accuracy. With the increase of sampling rate, the calculation efficiency is greatly improved, while still maintaining high reliability and accuracy.

[0122] Table 1

[0123]

[0124] To further verify the performance of this method, the sampled point cloud was reconstructed using a greedy triangle algorithm and a Poisson surface reconstruction algorithm (both existing techniques). Then, the readily available software CloudCompare was used to calculate the surface area of ​​each mesh. Figure 12 As shown, when the surface is smooth and the point cloud is uniformly distributed, the greedy triangle algorithm is insensitive to the sampling rate; in this case, the area calculated by the greedy triangle algorithm is highly reliable. The Poisson surface reconstruction method can fit and approximate all points to generate a closed surface with watertightness and good geometric features. Although Poisson surface reconstruction is a good overall reconstruction method, it is only suitable for smooth, closed surfaces without sharp features. For fairness, only the reconstruction results of the first two scenarios are used for comparative experiments. Figure 13 (a) and Figure 13 As shown in (b), although the surface area calculation based on Poisson surface reconstruction has high accuracy for specific objects with smooth surfaces, this method is very time-consuming. The method proposed in this invention is an order of magnitude (more than 10 times) faster than the method based on Poisson surface reconstruction. In reality, due to various factors, the surface of point clouds is usually filled with a large amount of random noise. To test the robustness of the proposed method, area calculation is then performed from noisy point clouds, as follows... Figure 14As shown in (a), even with significant noise, the method proposed in this invention can perform area calculations well, demonstrating its good noise resistance. Since the greedy triangle algorithm is highly sensitive to noise, only a few small Gaussian noise values ​​were set here. Nevertheless, as the noise increases, the calculation accuracy decreases rapidly (see...). Figure 14 (b)). For example Figure 14 As shown in (c), the Poisson surface reconstruction algorithm is insensitive to Gaussian noise, but it is not applicable to all types of surfaces; more precisely, it is only applicable to certain specific surfaces. Comparison with existing technologies shows that the method of this invention not only significantly improves computational speed but also has better noise resistance and computational accuracy, and is applicable to all types of surfaces.

[0125] like Figure 15 As shown in (a), in December 2020, we used Faro Focus at Wuhan University. s A 150 laser scanner was used to scan a sculpture named "3S" at a sampling rate of 1 / 20. The sculpture consists of planar and curved sections, with a central quadrangular prism. To facilitate comparison and avoid errors caused by point cloud registration, only point cloud data acquired from one scanning angle is used here, such as... Figure 15 As shown in (b), only the flat middle section is retained, with a total of 130,308 points. Next, a Leica AT960 laser tracker is used to obtain the precise coordinates of the six corners of the target area. The laser tracker's point measurement accuracy is 15 ± 6 micrometers, and the measured area is 6.9218 m². 2 Note the parameter ε. d Set the height to 0.5m and use triangles to calculate the area. For example... Figure 15 As shown in (c), a low-altitude UAV based on the Changzhou Xinyi X650 flight platform and a Regel Mini300 laser scanner was used to collect three-dimensional color point clouds of Hubei University of Science and Technology, and a relatively flat football field was selected as the target area (see...). Figure 15 (d) The target area has 2,974,086 points. Then, using a Topcon total station, the coordinates of the four corners of the football field were measured, yielding an area of ​​7154.722 m². 2 Note that because the point cloud in the target region is relatively sparse, the parameter ε d The area was set to 5m, and a sector was used to calculate the area. The area calculation results for the two target regions are listed in Table 2. It can be seen that the method proposed in this invention 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.

[0126] Table 2

[0127]

[0128] like Figure 16 As shown, to test the effectiveness of the proposed method on complex surfaces, two geometric models named Horse and Skeleton Hand are presented in detail. Both models were created by the Large Geometric Model Archive at Georgia Institute of Technology. The first model is magnified 10 times to approximate its true scale, while the second model remains unchanged. Reference values ​​are calculated using CloudCompare. Parameter ε d The area was calculated using a sector, with a value set to 0.2m. As shown in Table 3, the method proposed in this invention still achieves good results even for very complex curved surfaces.

[0129] Table 3

[0130]

[0131] The above embodiments are merely illustrative examples of the technical solutions of the present invention. The method for rapidly 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 by the claims. Any modifications, additions, or equivalent substitutions made by those skilled in the art based on these embodiments are within the scope of protection claimed by the claims of the present invention.

Claims

1. A method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data, characterized in that, Includes the following steps: Step 1. Obtain the surface point cloud of the target object to be measured as the original point cloud P. orig An adaptive supervoxel segmentation algorithm is used to segment the original point cloud P. orig Divided into a series of small pieces F fac For each facet F fac Perform a fitting operation and project the points on the patch onto the corresponding precisely fitted plane; obtain the original point cloud P. orig The methods include, but are not limited to, any one of: oblique photography and 3D laser scanning; Step 2. Let F porj For F fac The projected point cloud, for F porj Perform specially designed rigid body transformations: using F porj The center point on the graph is the origin of the local coordinate system, and the coordinates of the center point and F are... porj The normal vector of the first point is the local coordinate system X-axis, with F porj The normal vector is the Z-axis of the local coordinate system, and the cross product of the X-axis and Z-axis of the local coordinate system is the Y-axis of the local coordinate system. A rigid body transformation is performed to make F... porj The normal vector is parallel to the Z-axis of the global coordinate system; Step 3. Let F porj The point cloud obtained after rigid body transformation is F tran After rigid body transformation, a plane with arbitrary orientation angles in space will be parallel to the XY plane, and the z-coordinate value of each point will be the same. Then, F is extracted. tran edge point F in bou ; Step 4. Calculate F bou center point p cen ,make = (0, 1) is the starting direction, for F bou Any point p in bou Calculate the result from p cen and p bou A defined normal vector Then calculate and The angle between , calculate Sort and store in a predetermined order In the process, the order of the corresponding points is also adjusted and stored in P. sort In the middle, for the sorted edge points P sort If the distance between two adjacent edge points is greater than a certain threshold, then interpolation is performed to encrypt the points, and the interpolated points are reordered to obtain a new set of edge points P. interp ; Step 5. For P interp Any two adjacent points p in i and p j All were calculated from p i p j and p cen The area S of the triangle or sector formed area The surface area S of the target object is obtained by adding up the areas of all triangles or sectors within all the facets. all .

2. The method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data according to claim 1, characterized in that: in, During the supervoxel segmentation process in step 1, the number of neighboring points of each point is calculated. Set to 20, if P orig If the proportion of planar structural features reaches more than 50%, then the distance scale will be... Set it to 0.4~0.6m, otherwise change the distance scale. Set to 0.2~0.3m.

3. The method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data according to claim 1, characterized in that: in, In step 1, for each generated patch F fac A robust plane fitting method based on PCA is used to fit it, and F is set to... fac The equation of the plane containing the point in the middle is: In the formula, A, B, and C are three normal vectors of the plane and satisfy... ; F fac A little bit above p i Distance d to the plane i for: In the formula, x i y i and z i For point p i The coordinate values; objective function The following conditions must be met: In the formula, N is F fac The number of points in the middle; Solving the objective function using the Lagrange multiplier method To find the minimum value, first construct the following function: In the formula, These are parameters to be determined. Taking the partial derivative of equation (4) with respect to D and setting the result to 0, we obtain the following equation: Then we have: In the formula ， , , ; Equation (2) can be rewritten as follows: Using equation (4) with respect to A, B, and C respectively, we can obtain the following relationship: In the formula, , , ; make because It is always symmetric positive semi-definite, let , , They are respectively eigenvectors, ≥ ≥ They are respectively eigenvalues, Decomposed into: Calculate F fac Each point p in i From and F fac The distance t of the plane is determined by the centroid. i Then the standard deviation for: If t i >2 Then the corresponding point will be moved from F. fac Remove from the middle, and recalculate using the remaining points. And the centroid; set the number of iterations to ,go through The exact values ​​of the plane normal vectors A, B and C are obtained through several iterations, and then the value of D is obtained using equation (6). F fac The projection equation of a point on the plane onto the fitted plane is: In the formula, x i y i and z i The coordinates before projection, x j y j and z j These are the projected coordinate values.

4. The method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data according to claim 1, characterized in that: in, In step 2, F porj The center point is: In the formula, x porj i y porj i and z porj i For F porj The coordinates of any point in N, where N is F porj The number of points in the middle; Next, F porj The local coordinate system on is defined as follows: In the formula, , and These represent the directions of the three coordinate axes of the local coordinate system, x1, y1, and z1, respectively, and represent F. porj The coordinates of the first point in the middle; The target coordinate system is defined as follows: In the formula, ， ； Let M 3×3 If M is a rotation matrix, then 3×3 The calculation is as follows: Let p i (x i , y i , z i ) is F porj At the middle point, p j (x j , y j , z j ) for p i (x i , y i , z i The point obtained after performing a 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 ​​complex surfaces based on point cloud data according to claim 1, characterized in that: in, In step 4, F is first calculated. bou center point p cen : In the formula, x bou i and y bou i F respectively bou The coordinates of a point in the middle, M is F bou The number of points in the middle; make = (0, 1) is the starting direction, for F bou Any point p in bou , by p cen and p bou A defined normal vector The calculation is as follows: In the formula, 1≤i≤M; and The angle between for: In the formula, ； The calculated angle Sort by clockwise or counterclockwise and store In the case where 1≤i≤M, the order of the corresponding points is adjusted accordingly and stored in P. sort middle.

6. The method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data according to claim 1, characterized in that: in, In step 4, for the sorted edge points P sort An improved cubic B-spline curve is used for interpolation calculation, and this cubic B-spline curve passes through four given control points P1~P4: In the formula, P i As feature points controlling the curve, P1~P4 correspond to the four points adjacent to the current edge point P0, respectively, where P0=P sort Fi(t) represents the basis function of a cubic B-spline curve; t For F i The independent variable in (t) takes any real number between 0 and 1; For a sparse edge point, an improved cubic B-spline interpolation is used to densify it, ensuring that the distance between any two interpolated points is no greater than a threshold. The interpolated points are then reordered to obtain a new set P of edge points. interp .

7. The method for rapidly calculating the surface area of ​​complex surfaces based on point cloud data according to claim 1, characterized in that: in, In step 5, let p interp i and p interp j It is P interp If two adjacent points in the target surface planar structure feature account for more than 50% of the total area, a triangle is used for calculation, based on p. interp i p interp j and p cen The area S of the triangle formed area for: In the formula, In the formula, (x interp i , y interp i ) and (x interp j , y interp j ) are p respectively interp i and p interp j The coordinates of the point; If the surface structure features of the target surface to be calculated account for more than 50%, a sector is used for calculation, based on p. interp i p interp j and p cen The area S of the sector formed area The calculation is as follows: in, ， The surface area S of the target object all The calculation is as follows: In the formula, S represents the surface F. fac The total number, R is the current m-th face F fac P interp The number of points in the middle, For the current m-th face F fac The area of ​​the nth triangle or sector in the equation.

Images