CN115760954A - Method for rapidly calculating surface area of complex surface based on point cloud data - Google Patents
Method for rapidly calculating surface area of complex surface based on point cloud data Download PDFInfo
- Publication number
- CN115760954A CN115760954A CN202211256840.5A CN202211256840A CN115760954A CN 115760954 A CN115760954 A CN 115760954A CN 202211256840 A CN202211256840 A CN 202211256840A CN 115760954 A CN115760954 A CN 115760954A
- Authority
- CN
- China
- Prior art keywords
- point
- points
- formula
- bou
- porj
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Landscapes
- Image Generation (AREA)
Abstract
The invention provides a method for quickly calculating the surface area of a complex surface based on point cloud data, which comprises the following steps: step 1, dividing original point cloud of a complex surface to be calculated into a series of small surface patches, and projecting points to a fitting surface; step 2, projecting the point cloud F porj Performing specially designed rigid body transformation to obtain point cloud F tran (ii) a Step 3. Extraction of F tran Edge point F in (1) bou (ii) a Step 4, calculating F bou Central point p of cen For F bou At an arbitrary point p in bou Calculate by p cen And p bou Determined normal vector
Then calculate the included angle theta i Sorted and stored at θ sort In the method, aiming at edge points which are sequenced, reordering is carried out after interpolation to obtain a new edge point set; and 5, calculating the area of a triangle or a sector surrounded by the center point for any two adjacent points in the edge point set, and accumulating all the surface patches to obtain the surface area of the complex surface.
Description
Technical Field
The invention belongs to the technical field of three-dimensional laser scanning, and particularly relates to a method for rapidly calculating the surface area of a complex surface based on point cloud data.
Background
In actual work, the surface area of an object is often required to be calculated, and the calculation of the surface area plays an important role in the fields of medical treatment, land survey, agriculture and the like. The surface area of a regular object can be calculated by simple mathematics, while the surface area of an irregular object is calculated by more complicated methods, and the results obtained by different methods are not as same as each other. Conventionally, the surface area is calculated mainly using various measuring instruments such as a GPS, a total station, a laser tracker, and a laser range finder, and the area is calculated by a curved surface fitting method, a curve fitting method, an irregular network method, and the like, based on the obtained discrete points and side lengths. However, the precision of single-point measurement is generally limited by the precision of laser ranging and angle measurement of scales, and the effects of ranging and angle measurement errors are difficult to reduce. These methods are inefficient, have poor accuracy, and are computationally expensive, especially for large scale, rough, curved target objects that are difficult to access.
Disclosure of Invention
The invention aims to solve the problems of low efficiency, poor precision, high cost and the like in the surface area calculation of large-scale complex surfaces, and aims to provide a method for quickly calculating the surface area of the complex surface based on point cloud data, which is particularly suitable for target objects which are difficult to contact, large in scale, rough and bent.
In order to achieve the purpose, the invention adopts the following scheme:
the invention provides a method for rapidly calculating the surface area of a complex surface based on point cloud data, which is characterized by comprising the following steps of: step 1, adopting a self-adaptive hyper-voxel segmentation algorithm to carry out point cloud P of 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;
The method is mainly used for calculating the surface area on the unorganized point cloud, is irrelevant to the source mode of point cloud data, and can be used for acquiring the required data in the steps. For example, for some target objects which are difficult to access, large-scale, rough and curved, the accurate surface area of the target objects is difficult to obtain by a common measurement mode, a three-dimensional laser scanner is firstly used for obtaining the point cloud of the surface of the object, and then the surface area is calculated by the method.
Further, the method for rapidly calculating the surface area of the complex surface based on the point cloud data provided by the invention also has the following characteristics: in the process of super-voxel segmentation in step 1, the number epsilon of adjacent points of each point is determined ) Set to 20 if P orig If the structural feature ratio of the midplane 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.
Further, the method for rapidly calculating the surface area of the complex surface based on the point cloud data provided by the invention can also have the following characteristics: in step 1, for each patch F generated fac Fitting the plane by adopting a robust plane fitting method based on PCA (Principal Component Analysis), and setting F fac The plane of the point in (1) is as follows:
Ax+By+Cz+D＝0 (1)
wherein A, B and C are three normal vectors of a plane and satisfy A 2 +B 2 +C 2 ＝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 Coordinate values of (2);
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 pending parameter;
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 (I), the compound is shown in the specification,
formula (2) can 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 (I), the compound is shown in the specification,
order to
Due to M 3×3 Is always symmetrical and semi-positive, such that
Are each M 3×3 Of the feature vector λ 1 ≥λ 2 ≥ λ 3 Are respectively M 3×3 Characteristic value of (D), M 3×3 The decomposition is as follows:
to obtain more accurate fitting parameters, F is calculated fac Each point p in (1) i To by
And F fac Determines the distance t of the plane i Then the standard deviation σ is:
if t i If > 2 sigma, the corresponding point is changed from F fac Removing the residual points, and recalculating with the residual points
And a centroid point; setting the number of iterations to epsilon 1 Through epsilon 1 The precise values of the plane normal vectors A, B and C are obtained through the secondary iteration, and then the value of D is obtained through the formula (6);
f is to be fac The point on the projection line is 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.
Further, the method for rapidly calculating the surface area of the complex surface based on the point cloud data provided by the invention also has the following characteristics: in step 1, raw point cloud data P is given orig Managing P with kd-Tree orig And establishing corresponding index, and utilizing self-adaptive hyper-voxel segmentation algorithm to divide P orig Splitting into a series of small patches, where a patch involves two parameters to be adjusted: number of neighboring points of each point ε k And a distance scale ε d And is of d Can be regarded as the average side length, epsilon, of a patch k Typically set at 20. Epsilon d Cannot be determined by a fixed value, but depends on the particular situation. If the original point cloud data consists mainly of planar structures, ε should be determined d Set to be slightly larger (e.g., 0.5 m), if the point cloud contains rich surface features, ε should be adjusted d Set slightly smaller (e.g., 0.25 m).
Further, the method for rapidly calculating the surface area of the complex surface based on the point cloud data provided by the invention can also have the following characteristics: in step 2, F porj The central point of (a) is:
in the formula, x porj i ，y porj i And z porj i 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 (I), the compound is shown in the specification,
and
respectively representing the directions of three coordinate axes of a local coordinate system, x 1 ，y 1 And z 1 Is shown as F porj Coordinate value of the first point;
the target coordinate axis system is defined as follows:
in the formula (I), the compound is shown in the specification,
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:
further, the method for rapidly calculating the surface area of the complex surface based on the point cloud data provided by the invention can also have the following characteristics: in step 3, assume F porj The point cloud obtained after rigid body transformation is F tran After rigid body transformation, a plane with any attitude angle in space will be parallel to the X-Y 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, F is extracted using a two-dimensional alpha-shape algorithm tran Edge points (only x and y coordinate values participate in the calculation). In general, the α -shape algorithm is a method for extracting edge points in a discrete point set, and the principle can be regarded as that a radius is epsilon 2 Is surrounded by F tran Rolling outside, the ring passing through F tran Any two points in the table if there are no other points in the tableIn the circle, these two points are then considered as edge points, which are iteratively obtained in this way. Let the obtained edge point be F bou . It should be noted that, since the average point densities of different point cloud patches may not be the same, and sometimes the difference is large, in order to extract the edge points better, epsilon is set 2 When this parameter is (epsilon) 2 Representing the extraction F of alpha-shape algorithm tran Radius value of the circle at the middle edge point), which is usually made to match F tran Average dot spacing d of aver Correlation, e.g. set to F porj 10 times the average dot spacing. For F tran Any point p in (1) i Performing a k-neighbor search to find the corresponding n neighbor points (denoted as p) i 1 ,p i 2 …,p i n ) Then d is aver The following can be calculated:
in the formula, N is F tran The number of the points is that i is more than or equal to 1 and less than or equal to N. n may be set to 2.
Further, the method for rapidly calculating the surface area of the complex surface based on the point cloud data provided by the invention can also have the following characteristics: in step 4, the key to the area calculation is to be represented by F bou The formed polygon is divided into a series of small triangles or small sectors, and the total area can be obtained by counting the sum of the areas of the small triangles or the small sectors; therefore, before area calculation, the edge points need to be regularized; first, F is calculated bou Central point p of (2) cen ：
In the formula, x bou i And y bou i 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;
in the formula (I), the compound is shown in the specification,
the calculated angle theta i 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).
Further, the method for rapidly calculating the surface area of the complex surface based on the point cloud data provided by the invention can also have the following characteristics: in step 4, the edge points P are ordered sort If the distance between two adjacent edge points is greater than a certain threshold (e.g. 0.01 m), it is necessary to interpolate and encrypt them, and interpolate the 3-time B-spline curve by using the improved 3-time B-spline curve and making the 3-time B-spline curve pass through 4 given pointsControl point P 1 ～P 4 ；
In the formula, P i For controlling characteristic points of the curve, P 1 ～P 4 Respectively correspond to the current edge point P 0 Four adjacent points, P 0 ＝P sort (ii) a Fi (t) represents the basis function of the 3-degree B-spline curve; t is F i (t) argument, t ∈ [0,1 ]]I.e. t takes any real number between 0 and 1;
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, and reordering the interpolated points to obtain a new edge point set P interp 。
Further, the method for rapidly calculating the surface area of the complex surface based on the point cloud data provided by the invention can also have the following characteristics: in step 5, let p interp i And p interp j Is P interp If more than 50% of the target surface to be calculated consists of planar structures (the proportion of planar structure features is more than 50%), then a triangle is used for calculation, and the target surface to be calculated consists of three or more planar structuresp interp i 、p interp j And p cen The area S of the triangle area Comprises the following steps:
in the formula (I), the compound is shown in the specification,
in the formula (x) interp i ,y interp i ) And (x) interp j ,y interp j ) Are each p interp i And p interp j The point coordinates of (2).
If more than 50% of the target surface to be calculated is formed by the curved surface (the structural feature of the curved surface accounts for more than 50%), calculating by adopting a fan shape and p interp i 、p interp j And p cen Sector area S of formation area The calculation is as follows:
wherein, the first and the second end of the pipe are connected with each other,
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.
Action and effects of the invention
The method for quickly calculating the surface area of the complex surface based on the point cloud data firstly proposes that the surface area calculation is directly carried out on the unorganized point cloud without three-dimensional reconstruction, and the method is suitable for various complex surfaces.
Drawings
FIG. 1 is a flow chart of a method for rapidly calculating a surface area of a complex surface based on point cloud data according to an embodiment of the present invention;
FIG. 2 is a schematic diagram of voxel segmentation according to an embodiment of the present invention, in which (a) is original point cloud data and (b) is a voxel segmentation result;
FIG. 3 is a comparison diagram before and after projection, wherein (a) is before projection and (b) is after projection, according to an embodiment of the present invention;
FIG. 4 is a schematic representation of rigid body transformation according to 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 according to an embodiment of the present invention, wherein (a) is before extraction and (b) is after extraction;
FIG. 7 is a schematic plan projection diagram of an embodiment of the present invention, wherein (a) is before projection and (b) is after projection;
FIG. 8 is a schematic diagram of the effect of improved 3-degree B-spline interpolation on edge point encryption according to an embodiment of the present invention, where (a) is before encryption and (B) is after encryption;
FIG. 9 is a diagram illustrating the results of area calculations according to an embodiment of the present invention;
FIG. 10 is a simulated point cloud data, wherein (a) is a cube, (b) is a sphere, and (c) is a circle, according to an embodiment of the present invention;
FIG. 11 is a calculation result of the proposed method at different sampling rates according to an embodiment of the present invention, where (a) is calculation time and (b) is relative accuracy;
FIG. 12 is a graph of the results of a greedy triangle algorithm at different sampling rates, where (a) is computation time and (b) is relative precision;
fig. 13 is a calculation result of a poisson surface reconstruction algorithm according to an embodiment of the present invention at different sampling rates, where (a) is calculation time, and (b) is relative accuracy;
FIG. 14 shows the calculation results of the three methods according to the embodiment of the present invention under different Gaussian noises, wherein (a) is the method according to the present invention, (b) is a greedy triangle algorithm, and (c) is a Poisson surface reconstruction algorithm;
FIG. 15 shows a real experiment 1 according to an embodiment of the present invention, in which (a) is a "3S" sculpture, (b) is a target area 1, (c) is a stadium of Hubei scientific institute, and (d) is a target area 2;
FIG. 16 is a diagram of a real experiment 2 according to an embodiment of the present invention, in which (a 1) is an original geometric model of Horse, (b 1) is point cloud data of Horse, and (c 1) is a hyper-voxel segmentation result of Horse; (a2) The original geometric model of the Skeleton Hand, (b 2) the point cloud data of the Skeleton Hand, and (c 2) the super voxel segmentation result of the Skeleton Hand.
Detailed Description
The following describes in detail a specific embodiment of the method for rapidly calculating a surface area of a complex surface based on point cloud data according to the present invention with reference to the accompanying drawings.
< example >
As shown in fig. 1, the method for rapidly calculating the surface area of a complex surface based on point cloud data provided by this embodiment includes the following steps:
step 1To show the calculation process of the method of the present invention, a simulation point cloud of a cube shown in fig. 2 is first selected for example, and a simulation point cloud P is obtained by sampling six faces of a cube with a side length of 1m orig Each facet contains 26,896 points, with a total of 161,376 points. The hyper-voxel segmentation involves two parameters: nearest neighbor number epsilon of each point k And a distance scale ε d ，ε d Can be viewed as the average edge length of a patch. In our experiments,. Epsilon. k Set to 20. Epsilon d Cannot be determined by a fixed value, but as the case may be, if the original point cloud data contains abundant (more than 50%) planar features, ε should be d Set to be slightly larger (e.g., 0.5 m), if the point cloud data contains abundant (over 50%) surface features, ε should be increased d Set slightly smaller (e.g., 0.25 m). The proposed method involves two parameters: epsilon 1 And ε 2 In this example,. Epsilon 1 Set to 200, epsilon 2 Is set to F porj 10 times the average dot spacing. Unless otherwise stated, the above parameter settings will be used for all experiments referred to in this example without adjustment. By hyper-voxel segmentation, 28 surface patches F are formed fac . For each patch F generated fac Fitting by using a robust plane fitting method based on PCA (Principal Component Analysis), setting the iteration number to be 200, and fitting F fac All points above are projected onto the fitting surface, as shown in fig. 3, and the small patch with noise becomes very flat after projection, which is very beneficial for the extraction of edge points and the area calculation.
And 3, as shown in fig. 5, after rigid body transformation, a plane with any attitude angle in the space is parallel to the X-Y plane, which means that the z coordinate value of each point is the same, so that the three-dimensional data is reduced to two-dimensional data. Next, an improved adaptive alpha-shape is usede algorithm extraction of F tran Edge point F in (1) bou (only the x and y coordinate values participate in the calculation), as shown in fig. 6, the edge points of the patch are successfully extracted.
And 4, as shown in the figure 7, taking the north direction as the positive direction, sequencing the edge points clockwise, calculating the Euclidean distance between every two points after sequencing, and if the Euclidean distance is more than 0.01m, interpolating by using an improved 3-time B-spline curve to ensure that the distance between any two adjacent points is not more than 0.01m, and as shown in the figure 8, the Euclidean distance is an interpolation graph of sparse edge points.
As shown in FIG. 10, to further illustrate the reliability of the method, in addition to the cube, two additional simulated point clouds were selected for testing, where the radius of the Sphere is 1m, the total number of points is 160,000, ε d Set at 0.25m, torus inner diameter of 0.3m, outer diameter of 1m, total number of dots of 160,000,. Epsilon d Set to 0.5m. The point cloud is typically very dense and, to improve efficiency, it is often necessary to downsample the point cloud prior to segmentation. To test the effect 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 by the present invention on three simulated point clouds, from which the following conclusions can be drawn: the proposed surface area calculation method is very efficient and reliable in results, and can obtain calculation results highly consistent with theoretical values on both planar structures and curved surfaces (such as spheres and rings).
To test the effect of different sampling rates on the calculated results, three simulated point clouds were downsampled at five different density levels. As shown in fig. 11 (a) and 11 (b), proper sampling does not necessarily reduce relative accuracy compared to performing surface area calculation directly on the original point cloud, and as the sampling rate increases, the calculation efficiency greatly increases, and higher reliability and accuracy can still be maintained.
TABLE 1
To further verify the performance of the method, the sampled point cloud is reconstructed using a greedy triangle algorithm and a poisson surface reconstruction algorithm (both of which are prior art methods), and then the surface area of each mesh grid is calculated using the off-the-shelf software CloudCompare. As shown in fig. 12, the greedy triangle algorithm is insensitive to the sampling rate when the surface is smooth and the point clouds are evenly distributed, in which case the area result calculated by the greedy triangle algorithm is highly reliable. The poisson surface reconstruction method may fit and approximate all points to generate a closed surface with water tightness and good geometric characteristics. Although poisson surface reconstruction is a good overall reconstruction method, it is only applicable to smooth closed surfaces without sharp features. For the sake of fairness, only the reconstruction results of the first two scenarios were used for comparative experiments. As shown in fig. 13 (a) and 13 (b), although the calculated surface area based on the poisson surface reconstruction has a higher accuracy for a specific object with a smooth surface, the method is very time-consuming, and the method proposed by the present invention is faster by an order of magnitude (more than 10 times) than the method based on the poisson surface reconstruction. In fact, the surface of a point cloud is often saturated with a large amount of random noise, influenced by various factors. In order to test the robustness of the proposed method, area calculation is performed from the noisy point cloud, as shown in fig. 14 (a), even if there is a lot of noise, the proposed method of the present invention can perform area calculation well, which shows that the method has good noise resistance. Since the greedy triangle algorithm is very sensitive to noise, only a few small gaussian noise values are set here, but nevertheless, the calculation accuracy rapidly decreases as the noise increases (see fig. 14 (b)). As shown in fig. 14 (c), the poisson surface reconstruction algorithm is insensitive to gaussian noise, but the algorithm is not applicable to all types of surfaces, and precisely, is applicable only to certain specific surfaces. Compared with the prior art, the result shows that the method not only greatly improves the calculation speed, but also has better noise resistance and calculation precision, and is suitable for all types of surfaces.
As shown in FIG. 15 (a), in 12 months of 2020, we used Faro Focus at Wuhan university s The 150 laser scanner scanned a sculpture named "3S" with a sampling rate of 1/20, the sculpture consisting of a flat and curved portion, the middle portion being a quadrangular prism. To facilitate comparison and avoid errors caused by point cloud registration, only point cloud data acquired from one scan angle is used here, as shown in fig. 15 (b), only the middle flat portion is retained, with the number of points 130,308. Then, a Leica AT960 laser tracker is used for obtaining the precise coordinates of six angles of the target area, the point measurement precision of the laser tracker is 15 +/-6 microns, and the measured area is 6.9218m 2 . Note that the parameter ε d Set to 0.5m, the area is calculated using triangles. As shown in fig. 15 (c), a three-dimensional color point cloud of the northHubei academy of science was collected using a low altitude unmanned aerial vehicle based on the Changzhou new wing X650 flight platform and the Regel Mini300 laser scanner, and a relatively flat football field was selected as a target area (see fig. 15 (d)), with the number of target area points being 2,974,086. The coordinates of the four corners of the football field were then measured using a Topcon total station, which gave an area 7154.722m 2 . Note that the parameter ε is the result of sparse point clouds in the target region d Set to 5m, the area is calculated using a fan. The area calculation results of the two target areas are listed in table 2, and it can be seen that the method successfully calculates the area of the test point cloud data within an acceptable time, which is very close to an actual measurement value, thereby greatly improving the working efficiency.
TABLE 2
As shown in fig. 16, two geometric models named Horse (Horse) and Skeleton Hand (skeeleton Hand) were selected for detailed display in order to test the effect of the proposed method on complex surfaces. Both models were made by the large geometry model archive of the Georgia college of Physician, the first model being enlarged by a factor of 10 to approximate the true scale, and the second modelThe two models remain unchanged. Reference values were calculated from cloudbiare. Parameter epsilon d Set to 0.2m, the area is calculated using a sector. As shown in table 3, the method of the present invention can achieve good results even for very complex curved surfaces.
TABLE 3
The above embodiments are merely illustrative of the technical solutions of the present invention. The method for rapidly calculating the surface area of a complex surface based on point cloud data according to the present invention is not limited to the contents described in the above embodiments, but is subject to the scope defined by the claims. Any modification or supplement or equivalent replacement made by a person skilled in the art on the basis of this embodiment is within the scope of the invention as claimed in the claims.
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) i And p j All calculate by p i 、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 2 +B 2 +C 2 ＝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 (I), the compound is shown in the specification,
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 (I), the compound is shown in the specification,
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 λ 1 ≥λ 2 ≥λ 3 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 1 Through epsilon 1 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 i ，y porj i And z porj i 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 (I), the compound is shown in the specification,
and
respectively representing the directions of three coordinate axes of a local coordinate system, x 1 ，y 1 And z 1 Is represented by F porj Coordinate values of the first point;
the target coordinate axis system is defined as follows:
in the formula (I), the compound is shown in the specification,
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 i And y bou i 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;
in the formula (I), the compound is shown in the specification,
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 1 ～P 4 ：
In the formula, P i For controlling characteristic points of the curve, P 1 ～P 4 Respectively correspond to the current edge point P 0 Four adjacent points, P 0 ＝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 i 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 i 、p interp j And p cen The area S of the triangle area Comprises the following steps:
in the formula (I), the compound is shown in the specification,
in the formula (x) interp i ,y interp i ) And (x) interp j ,y interp j ) Are each p interp i 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 i 、p interp j And p cen Formed sector area S area The calculation is as follows:
wherein, the first and the second end of the pipe are connected with each other,
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.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211256840.5A CN115760954A (en) | 2022-10-14 | 2022-10-14 | Method for rapidly calculating surface area of complex surface based on point cloud data |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211256840.5A CN115760954A (en) | 2022-10-14 | 2022-10-14 | Method for rapidly calculating surface area of complex surface based on point cloud data |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115760954A true CN115760954A (en) | 2023-03-07 |
Family
ID=85351408
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211256840.5A Pending CN115760954A (en) | 2022-10-14 | 2022-10-14 | Method for rapidly calculating surface area of complex surface based on point cloud data |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115760954A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116839511A (en) * | 2023-07-28 | 2023-10-03 | 江苏建筑职业技术学院 | Irregular surface area measuring method based on photogrammetry technology |
-
2022
- 2022-10-14 CN CN202211256840.5A patent/CN115760954A/en active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116839511A (en) * | 2023-07-28 | 2023-10-03 | 江苏建筑职业技术学院 | Irregular surface area measuring method based on photogrammetry technology |
| CN116839511B (en) * | 2023-07-28 | 2024-03-19 | 江苏建筑职业技术学院 | Irregular surface area measuring method based on photogrammetry technology |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105654483B (en) | The full-automatic method for registering of three-dimensional point cloud | |
| CN107610061B (en) | Edge-preserving point cloud hole repairing method based on two-dimensional projection | |
| CN111696210A (en) | Point cloud reconstruction method and system based on three-dimensional point cloud data characteristic lightweight | |
| CN101751695A (en) | Estimating method of main curvature and main direction of point cloud data | |
| CN111415379B (en) | Three-dimensional point cloud data registration method based on cuckoo optimization | |
| WO2021203711A1 (en) | Isogeometric analysis method employing geometric reconstruction model | |
| CN110807781A (en) | Point cloud simplification method capable of retaining details and boundary features | |
| CN113628263A (en) | Point cloud registration method based on local curvature and neighbor characteristics thereof | |
| CN116543117B (en) | High-precision large-scene three-dimensional modeling method for unmanned aerial vehicle images | |
| CN103700135B (en) | A kind of three-dimensional model local spherical mediation feature extracting method | |
| CN115760954A (en) | Method for rapidly calculating surface area of complex surface based on point cloud data | |
| CN113538689A (en) | Three-dimensional model mesh simplification method based on feature fusion of neural network | |
| CN113012063A (en) | Dynamic point cloud repairing method and device and computer equipment | |
| Chen et al. | Denoising of point cloud data for computer-aided design, engineering, and manufacturing | |
| CN110688440B (en) | Map fusion method suitable for less sub-map overlapping parts | |
| CN114494641B (en) | Three-dimensional model light weight method and device | |
| CN109345571B (en) | Point cloud registration method based on extended Gaussian image | |
| CN113343328B (en) | Efficient closest point projection method based on improved Newton iteration | |
| CN111765883B (en) | Robot Monte Carlo positioning method, equipment and storage medium | |
| CN107818578B (en) | Rapid face model reconstruction algorithm and system based on registration method | |
| Maurya et al. | Performance of greedy triangulation algorithm on reconstruction of coastal dune surface | |
| CN108226926A (en) | A kind of three-dimensional scattering distribution reconstructing method based on radar network | |
| CN108387897B (en) | Projectile body positioning method based on improved Gauss Newton-genetic hybrid algorithm | |
| Dyshkant | Comparison of point clouds acquired by 3d scanner | |
| CN116447977B (en) | Round hole feature measurement and parameter extraction method based on laser radar |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |