Image reconstruction method for incomplete data of three-dimensional point cloud of stocktaking coal in coal yard

Abstract

The invention provides an image reconstruction method for three-dimensional point cloud incomplete data of coal stocktaking in a coal yard, and relates to the technical field of coal pile reconstruction. The method specifically comprises the following steps: acquiring initial three-dimensional point cloud data of a coal pile by using a laser radar; removing noise from the initial three-dimensional point cloud data by adopting a point cloud filtering algorithm and carrying out point cloud simplification by adopting a self-adaptive point cloud simplification algorithm in sequence, and constructing a triangular network topological structure among points in the simplified three-dimensional point cloud data by adopting an improved greedy triangulation algorithm; identifying all holes in the triangular network topological structure by using a hole identification algorithm, and adding patch curved surfaces at all the identified holes by using a minimum angle grid filling algorithm; and carrying out subdivision processing on the triangular network topological structure by using a grid subdivision algorithm, and carrying out smoothing processing on the triangular network topological structure after subdivision processing by using a grid smoothing algorithm to obtain repaired and reconstructed three-dimensional point cloud data, so that the image reconstruction requirement of the current power plant coal inventory three-dimensional point cloud incomplete data is met.

Description

Technical Field

[0001] The invention relates to the technical field of coal pile reconstruction, and in particular to an image reconstruction method for incomplete three-dimensional point cloud data of coal piles in a coal yard. Background Art

[0002] Thermal power generation is the main source of electricity production in my country. Coal, as the main energy source in thermal power plants, accounts for the vast majority of the cost of thermal power plants. The management of coal yards is directly related to the cost control and coal blending and combustion capabilities of thermal power plants. Research on the intelligence of coal yards has shown that improving the coal blending and combustion capabilities and simplifying the coal pile inventory process can effectively reduce costs. At present, thermal power plants mainly use manual methods to inventory coal yards. Although this method is simple, it is inefficient. The method of using bucket wheel excavator airborne laser scanning is highly accurate, but it is time-consuming and costly.

[0003] The 3D reconstruction technology based on point cloud is currently widely used in reverse engineering, data visualization, vision and medical technology. After long-term development, this technology has achieved remarkable results. According to the different methods of point cloud 3D reconstruction, it can be divided into parametric surface reconstruction, implicit surface reconstruction and mesh surface reconstruction. Among them, the parametric surface reconstruction method mainly includes B-spline surface reconstruction and NURBS surface reconstruction. Since the parametric surface reconstruction method requires parameterization of point cloud data, the parameterization effect of scattered point cloud data is poor, which may cause holes in the model during reconstruction. Therefore, the parametric surface reconstruction method is suitable for the reconstruction of objects with smooth surfaces and small texture changes. The implicit surface reconstruction methods mainly include Poisson algorithm and marching cube algorithm. Implicit surface reconstruction methods can better fit objects with simple structures into smooth surfaces, but it is difficult to achieve ideal reconstruction effects for objects with complex surface structures, and the reconstruction process takes a long time. The mesh surface reconstruction method is the most widely used in point cloud 3D reconstruction, and is often used to reconstruct objects with complex structures. The greedy projection triangulation algorithm is one of the main representative algorithms of the mesh surface reconstruction algorithm. Although the greedy projection triangulation algorithm can well reconstruct objects with complex topological structures, it is time-consuming and has the problem of holes in the model when reconstructing a large amount of point cloud data. It is necessary to use a point cloud hole repair algorithm for optimization.

[0004] At present, point cloud hole repair algorithms are divided into two categories: one is the hole repair for triangular mesh patch models, which completes the re-measurement of hole data by post-processing the triangular mesh model after triangulation of the scattered point cloud; the other is based on the scattered point cloud model for repair. The curved hole repair algorithm based on radial basis function increases the hole data filling points by continuously triangulating the hole area, but this method is more complicated and computationally intensive. Summary of the invention

[0005] In view of the above-mentioned deficiencies in the prior art, the present invention proposes an image reconstruction method for incomplete three-dimensional point cloud data of coal in coal yards based on an improved greedy triangulation algorithm and a hole repairing algorithm, aiming to solve the actual needs of image reconstruction of incomplete three-dimensional point cloud data of coal in power plants.

[0006] The present invention proposes an image reconstruction method for incomplete three-dimensional point cloud data of coal in a coal field, the method comprising the following steps:

[0007] Step 1: Collect initial 3D point cloud data of the coal pile using LiDAR;

[0008] Step 2: Use a point cloud filtering algorithm to remove noise from the initial three-dimensional point cloud data to obtain denoised three-dimensional point cloud data;

[0009] Step 3: Adopt an adaptive point cloud simplification algorithm to simplify the denoised three-dimensional point cloud data to obtain simplified three-dimensional point cloud data;

[0010] Step 4: Use the improved greedy triangulation algorithm to construct the triangular network topology structure between the points in the simplified three-dimensional point cloud data;

[0011] Step 5: Use hole identification algorithm to identify all holes in the triangular network topology structure;

[0012] Step 6: Use the minimum angle mesh filling algorithm to add patch surfaces to all identified holes in the triangulated network topology;

[0013] Step 7: Use a mesh subdivision algorithm to subdivide the added patch surface and the original surface in the triangular network topology structure;

[0014] Step 8: Use a mesh smoothing algorithm to smooth the triangular network topology structure after subdivision to obtain the repaired and reconstructed three-dimensional point cloud data;

[0015] The step 2 further comprises:

[0016] Step 2.1: Use straight-through filtering to remove invalid point cloud data in the initial three-dimensional point cloud data to obtain three-dimensional point cloud data after straight-through filtering;

[0017] Step 2.2: Using a statistical filtering method based on statistical principles to remove large-scale noise in the three-dimensional point cloud data after direct filtering, to obtain three-dimensional point cloud data after statistical filtering;

[0018] Step 2.2.1: For any point P in the three-dimensional point cloud data after direct filtering n , calculate point P n With point P nThe average distance and standard deviation of the K nearest neighbor points in the neighborhood; where n is the nth point in the 3D point cloud data, and n = 1,...,m; m represents the number of points in the 3D point cloud data;

[0019] Step 2.2.2: Calculate the average value of the distance between points and the variance of the distance between points in the three-dimensional point cloud data after the straight-through filtering by using the average distance and standard deviation of the K nearest neighbors of each point in the three-dimensional point cloud data after the straight-through filtering and the K nearest neighbors in the neighborhood of the point;

[0020] Step 2.2.3: Calculate the distance threshold according to the average value of the distance between points and the variance of the distance between points in the three-dimensional point cloud data after the straight-through filtering;

[0021] Step 2.2.4: Determine the average distance μ between each point in the three-dimensional point cloud data after direct filtering and the K nearest neighbor points in the neighborhood of the point n Is it greater than the distance threshold S? If so, remove the point; otherwise, keep the point to obtain the three-dimensional point cloud data after statistical filtering;

[0022] Step 2.3: Use the moving least squares method to remove small-scale noise in the three-dimensional point cloud data after statistical filtering to obtain denoised three-dimensional point cloud data;

[0023] Step 2.3.1: Construct the fitting function;

[0024] Step 2.3.2: Select the cubic spline function as the basis function of the fitting function;

[0025] Step 2.3.3: Perform least squares fitting on the neighboring points in the neighborhood of the fixed point u, and smooth the three-dimensional point cloud data after statistical filtering by calculating the local curvature and normal vector of the fixed point to obtain the denoised three-dimensional point cloud data;

[0026] The step 3 further comprises:

[0027] Step 3.1: For any point in the denoised 3D point cloud data, take the point as the target point, use KD-Tree to perform K nearest neighbor search on the target point, and obtain K nearest neighbor points of the target point;

[0028] Step 3.2: Use the least squares method to fit the target point and its K neighboring points, calculate the curvature of the target point, and obtain the curvature of all points in the denoised three-dimensional point cloud data;

[0029] Step 3.3: Set the curvature threshold, and use the improved octree voxel filtering algorithm to simplify the point cloud of all points whose curvature is greater than the curvature threshold in the denoised three-dimensional point cloud data to obtain simplified three-dimensional point cloud data;

[0030] The step 3.2 further comprises:

[0031] Step 3.2.1: Construct the surface parametric equation;

[0032] Step 3.2.2: For any point (α0, β0) on the surface Q, construct the second basic formula of the surface, and then calculate the distance between the point (α0, β0) and its K nearest neighbors;

[0033] The second basic formula of the surface is expressed as:

[0034]

[0035] Where L, M and N are the second-class basic parameters of the surface Q; α and β are variables; t αα (α0, β0) is the value of the second-order partial derivative of the surface parameter equation with respect to α at the point (α0, β0); t ββ (α0, β0) is the value of the second-order partial derivative of the surface parameter equation with respect to β at the point (α0, β0); t αβ (α0, β0) is the value of the second-order partial derivative of the surface parameter equation with respect to α and β at the point (α0, β0); h is the normal vector of the point (α0, β0); is the distance between the point (α0, β0) and its K nearest neighbors; Δα and Δβ are the increments of α and β respectively; o(Δα 2 ,Δβ 2 ) is an infinitesimal term;

[0036] Step 3.2.3: Calculate the mean curvature of the surface;

[0037] The average curvature of the surface is expressed as:

[0038]

[0039] Where H represents the average curvature of the surface; G, E and F are the first-class basic parameters of the surface Q;

[0040] Step 3.2.4: Determine the curvature of the target point based on the average curvature of the surface;

[0041] Step 3.2.5: Repeat steps 3.2.1-3.2.4 to obtain the curvature of all points in the denoised 3D point cloud data;

[0042] The step 3.3 further comprises:

[0043] Step 3.3.1: Set the curvature threshold and filter out all points whose curvature is greater than the curvature threshold from the denoised 3D point cloud data;

[0044] Step 3.3.2: Construct the first cube based on the maximum size of the point cloud space containing all the filtered points. The length, width, and height of the cube are recorded as X 长 , Y 宽 , Z 高 , the volume is V;

[0045] Step 3.3.3: Divide the first cube into C sub-cubes, calculate the centroid of each sub-cube, and mark the centroid coordinates as (x0, y0, z0); where x0, y0, and z0 are the horizontal, vertical, and vertical coordinates of the centroid respectively;

[0046] Step 3.3.4: Calculate the distance between each point in each sub-cube and the center of gravity of the sub-cube, and only keep the point closest to the center of gravity and the point farthest from the center of gravity to obtain the simplified three-dimensional point cloud data;

[0047] The coordinates of the point closest to the center of gravity and the point farthest from the center of gravity are expressed as:

[0048]

[0049] where x max The horizontal coordinate of the point farthest from the center of gravity; y max The ordinate of the point farthest from the center of gravity; z max represents the vertical coordinate of the point farthest from the centroid; η represents the number of points in the subcube; γ represents the γth point in the subcube; represents the volume of each point; x γ represents the horizontal coordinate of the γth point; l max Indicates the farthest distance from the center of gravity; y γ represents the ordinate of the γth point; z γ represents the vertical coordinate of the γth point; x min The horizontal coordinate of the point closest to the center of gravity; y min The ordinate of the point closest to the center of gravity; z min Indicates the vertical coordinate of the point closest to the center of gravity; l min Indicates the closest distance to the center of gravity;

[0050] The step 4 further comprises:

[0051] Step 4.1: Perform a planar projection on the simplified three-dimensional point cloud data to obtain a point set on a two-dimensional plane;

[0052] Step 4.1.1: For any point p in the simplified 3D point cloud data, construct the tangent plane equation of the k-nearest neighboring points of point p according to the tangent plane normal vector of the surface fitted by point p and its k-nearest neighboring points at point p;

[0053] Step 4.1.2: Normalize the tangent plane normal vector and calculate the i-th nearest neighbor point p of point p i Distance to the tangent plane;

[0054] Step 4.1.3: Based on the least squares principle, calculate the sum of the squares of the distances from the k nearest neighbor points of point p to the tangent plane. When the sum of the squares is the smallest, solve the tangent plane equation of the k-nearest neighbor points of point p.

[0055] Step 4.1.4: Project the k nearest neighbor points of point p onto the tangent plane according to the solved tangent plane equation to obtain a point set on the two-dimensional plane;

[0056] Step 4.2: Use the Delaunay triangulation algorithm based on the traditional point-by-point insertion method to triangulate the point set on the two-dimensional plane to obtain the triangular network topology structure between each point in the simplified three-dimensional point cloud data;

[0057] Step 4.2.1: Construct a super triangle containing all the points in the point set on the two-dimensional plane, and determine the extreme points x of the point set on the two-dimensional plane in the positive and negative directions of the coordinate axis of the two-dimensional coordinate system. max1 , x min1 ,y max1 and min1 , and add a positive offset value offset, and offset>0, to construct the minimum bounding box of the two-dimensional point set; where x max1 is the extreme point of the point set on the two-dimensional plane in the positive direction of the X axis; min1 is the extreme point of the point set on the two-dimensional plane in the negative direction of the X axis; y max1 is the extreme point of the point set on the two-dimensional plane in the positive direction of the Y axis; min1 It is the extreme point of the point set on the two-dimensional plane in the negative direction of the Y axis;

[0058] Step 4.2.2: Traverse all points in the point set on the two-dimensional plane, and insert each point into an initial triangular grid in turn to obtain a new triangular grid;

[0059] Step 4.2.3: Locally optimize the new triangular grid to obtain an optimized triangular grid;

[0060] Step 4.2.4: Delete the edges and vertices of the super triangle in the optimized triangular mesh to obtain the triangular network topology structure between the points in the simplified three-dimensional point cloud data;

[0061] The step 5 further comprises:

[0062] Step 5.1: Obtain the edges of all triangles in the triangular network topology and the face neighbor relationship of each edge, and divide all the obtained edges into mesh internal edges and mesh boundary edges according to the face neighbor relationship of each edge;

[0063] Step 5.2: Call the PolygonMesh data structure to obtain the connection relationship between the vertices and edges of all triangles in the triangular network topology structure, traverse all vertices and filter out the vertices that are not referenced by any triangle, and identify the holes in the triangular network topology structure based on the connection relationship between the filtered vertices and edges;

[0064] The step 6 further comprises:

[0065] Step 6.1: Obtain the coordinates of the mesh boundary edges and boundary vertices of the hole, calculate the Euclidean distance between two boundary vertices in each mesh boundary edge in space to obtain the length of the mesh boundary edge, and calculate the average length of all mesh boundary edges;

[0066] Step 6.2: For each vertex of the mesh boundary edge, calculate the adjacent edge angle between the two adjacent edges corresponding to the vertex;

[0067] Step 6.3: Filter out the boundary vertex with the smallest adjacent edge angle value, and calculate the distance D between the two boundary vertices adjacent to the boundary vertex;

[0068] Step 6.4: Add triangles to the triangular network topology structure according to the distance D;

[0069] Step 6.5: Update the coordinates of the boundary vertices of the hole and delete the boundary vertex Pi that has been repaired;

[0070] Step 6.6: Determine whether there are any unrepaired boundary vertices. If so, repeat steps 6.2-6.5; if not, the hole filling is completed;

[0071] The step 7 further comprises:

[0072] Step 7.1: Add a new vertex P' to the triangular network topology structure after the patch surface is added, and calculate the position of the new vertex P' according to the positions of all vertices in the triangular network topology structure after the patch surface is added;

[0073] Step 7.2: Add new edge points to the triangular network topology structure after adding the patch surface;

[0074] Step 7.3: For each new vertex P', connect point P' to all new vertices adjacent to point P' to obtain new edges; for each triangle, connect the three new edge points in the triangle to form a new edge;

[0075] Step 7.4: Every three interconnected new edge points and new vertices are combined into a new triangle to complete the subdivision of the triangular network topology structure;

[0076] The step 8 further comprises:

[0077] Step 8.1: Establish energy function;

[0078] The energy function is expressed as:

[0079]

[0080] Where E(f) represents the energy function; (v,w) represents an integration interval; f ss represents the second-order derivative; s is the independent variable;

[0081] Step 8.2: Use variational calculus to solve the functional function f that minimizes the energy function;

[0082] Step 8.2.1: Assume that f is the minimum value of the energy function E(f), and set is any function in the interval (v,w), and Satisfies first-order differentiability in the interval (v,w) and satisfies

[0083] Step 8.2.2: Assume there exists a function Where ζ is the coefficient;

[0084] Step 8.2.3: Assume that f(s) is the minimum value of the energy function E(f). When ζ = 0, the energy function reaches its minimum value, and the derivative of the energy function must also be zero, which is recorded as:

[0085]

[0086] Formula (33) is transformed into:

[0087]

[0088] in express f represents f(s);

[0089] Step 8.2.4: Calculate the derivative with respect to ζ;

[0090] The derivative with respect to ζ is expressed as:

[0091]

[0092] in express The second derivative of ss represents the second-order derivative of f;

[0093] Step 8.2.5: When ζ = 0, the right side of the equation must satisfy the condition of being equal to zero in the derivative with respect to ζ, which is written as:

[0094]

[0095] Step 8.2.6: Integrate formula (36) step by step to obtain:

[0096]

[0097] where f sss represents the third-order derivative of f(s), express

[0098] Step 8.2.7: Combine Simplifying formula (37) yields:

[0099]

[0100] Step 8.2.8: Integrate formula (38) step by step to obtain:

[0101]

[0102] where f ssss represents the fourth-order derivative of f(s);

[0103] Step 8.2.9: According to Simplifying formula (39) yields:

[0104]

[0105] Step 8.2.10: When formula (38) is satisfied, And the only case where the function is differentiable is:

[0106] f ssss =0 (41)

[0107] Step 8.2.11: Equivalently transform formula (39) into formula (41) to obtain:

[0108] ΔΔf=Δ 2 f=0 (42)

[0109] where Δ 2 is a double Laplacian operator;

[0110] Step 8.2.12: For the refined patch surface in the triangular network topology structure, the Laplace-Beltrami algorithm in the libigl library is called to solve formula (40) to obtain the smoothed triangular network topology structure and complete the repair and reconstruction of the three-dimensional point cloud data.

[0111] The beneficial effects of adopting the above technical solution are:

[0112] (1) Due to the complex scene environment of a thermal power plant, the point cloud data obtained by the lidar often carries a large amount of feature information, and there are also many points with repeated features and no features in the obtained point cloud data. Therefore, the method of the present invention proposes a lightweight greedy triangulation algorithm by improving the greedy triangulation algorithm, which can better restore the original contour of the point cloud data.

[0113] (2) In order to solve the problems of complex scenes around the target and uneven distribution of measurement data in the laser radar measurement of thermal power plants, the method of the present invention pre-processes the incomplete data of the obtained pan coal point cloud to obtain point cloud data with obvious smooth features.

[0114] (3) In view of the problems of large computational complexity and distortion of the simplified model in many current point cloud simplification algorithms, and considering that coal is often piled into regular trapezoids in actual applications with a relatively regular shape, the method of the present invention proposes an improved octree voxel filtering algorithm for simplifying the point cloud data after preliminary processing, thereby reducing the complexity of the point cloud data while trying to maintain the characteristics of the original data and reducing the computational complexity.

[0115] (4) In order to solve the problem that in the process of triangulating a planar region using a greedy triangulation algorithm, the computer performance is severely occupied and the calculation time is long because all points need to be traversed multiple times when using a region growing algorithm for triangulation, the method of the present invention proposes to project the three-dimensional point cloud data onto a two-dimensional plane, then use a point-by-point interpolation method to triangulate the two-dimensional plane, and map the topological structure on the two-dimensional plane back to three-dimensional space according to the depth value of the point cloud data.

[0116] (5) In order to solve the problem that the curvature change at the interface between the patch surface and the original surface is not taken into account in the traditional hole patching algorithm, resulting in a large difference in curvature between the patched surface and the original surface and poor patching effect, the method of the present invention proposes an improved hole patching algorithm. First, a hole recognition algorithm is used to identify and classify all holes. Secondly, a minimum angle grid filling algorithm is used to add patch surfaces to all holes. Then, a grid subdivision algorithm is used to subdivide the patch surfaces to obtain more operable vertices. Finally, a grid smoothing algorithm is used to make the transition of the patch surface smoother, solving the problem of large curvature change at the interface.

[0117] (6) The hole patches generated by the minimum angle mesh filling algorithm often produce triangular facets with narrow and long shapes and uneven sizes, which are inconsistent with the surface contour of the original model and have obvious transition problems. The method of the present invention adopts a dual Laplace operator deformation method to change the mesh shape at the hole patch, so that the added surface is smoother and closer to the surface contour of the original model.

[0118] In summary, the method of the present invention proposes an image reconstruction method and a hole repairing technology based on the residual point cloud data of the coal yard, that is, a lightweight greedy triangulation surface reconstruction algorithm and an improved hole repairing algorithm are proposed, which not only makes the greedy triangulation algorithm more lightweight and can better restore the original contour of the point cloud data, but also can automatically identify the edges of the holes and realize parallel repairing of multiple holes, which greatly improves the efficiency of hole repairing and has foreseeable greater economic value and social value. BRIEF DESCRIPTION OF THE DRAWINGS

[0119] Figure 1 This is a flow chart of an image reconstruction method for incomplete three-dimensional point cloud data of coal in a coal field in this embodiment;

[0120] Figure 2 This is a flow chart of a lightweight greedy triangulation algorithm in this implementation mode;

[0121] Figure 3 1 is a comparison diagram of the effects before and after the point cloud filtering algorithm is applied in this embodiment; wherein Figure (a) is the effect diagram before the point cloud filtering algorithm is applied; Figure (b) is the effect diagram after the point cloud filtering algorithm is applied;

[0122] Figure 4 This is a flow chart of the point cloud simplification algorithm in this implementation;

[0123] Figure 5 This is a schematic diagram of the two-dimensional projection of the three-dimensional point cloud data of Pan Coal in this implementation method;

[0124] Figure 6 is a schematic diagram of triangulating a point set on a two-dimensional plane in this implementation mode;

[0125] Figure 7 is a flow chart of the hole repair algorithm in this implementation mode;

[0126] Figure 8 This is a triangulated flow chart of the present embodiment;

[0127] Fig. 9 A flowchart of grid filling in this implementation mode;

[0128] Fig.10 Schematic diagram of adding triangular facets to the triangular network topology structure in this embodiment; FIG. (a) is a diagram satisfying D<1.2×lavg Schematic diagram of adding triangles to the triangular network topology structure; Figure (b) does not satisfy D < 1.2 × l avg Schematic diagram of adding triangles to the triangular network topology;

[0129] Fig.11 Schematic diagram of mesh subdivision in this implementation manner; wherein Figure (a) is a schematic diagram of adding new vertices in mesh subdivision; Figure (b) is a schematic diagram of adding new edge points in mesh subdivision; and Figure (c) is a schematic diagram of adding new triangles in mesh subdivision. DETAILED DESCRIPTION

[0130] In order to facilitate the understanding of the present application, the specific embodiments of the present invention are further described in detail below in conjunction with the accompanying drawings and embodiments. The following embodiments are used to illustrate the present invention, but are not intended to limit the scope of the present invention. On the contrary, the purpose of providing these embodiments is to make the disclosure of the present application more thoroughly understood.

[0131] The present embodiment provides an image reconstruction method for incomplete three-dimensional point cloud data of coal in a coal field, such as Figure 1 As shown, the method comprises the following steps:

[0132] Step 1: Use LiDAR to collect initial 3D point cloud data of the coal pile.

[0133] Step 2: Use the point cloud filtering algorithm to remove the noise in the initial 3D point cloud data to obtain denoised 3D point cloud data.

[0134] Step 2.1: Use straight-through filtering to remove invalid point cloud data in the initial three-dimensional point cloud data to obtain three-dimensional point cloud data after straight-through filtering.

[0135] In this embodiment, invalid point cloud data mainly refers to point cloud data that is not related to the coal pile, that is, three-dimensional point cloud data that does not belong to the coal pile. For example, the conveyor belt next to the coal pile. When collecting point cloud data, since the positions of the conveyor belt and the coal pile are close, it is inevitable to collect data from the conveyor belt when collecting data, so it is necessary to remove this part of invalid point cloud data first.

[0136] In this embodiment, if Figure 2As shown, for the collected initial three-dimensional point cloud data, the coordinate range of the initial three-dimensional point cloud data is obtained in the two specified directions x and y respectively, and the threshold range of the point cloud data coordinates belonging to the coal pile in the x and y directions is determined according to the obtained coordinate range. Each point in the initial three-dimensional point cloud data is traversed to determine whether the coordinate value of the point in a specified direction is within the threshold range. If so, the point is retained; if not, the point is deleted; in this way, all point cloud data that obviously does not belong to the coal pile are removed. The points left after the traversal are used as filtered point cloud data. For example, in the x direction, the coordinate range of the initial three-dimensional point cloud data obtained is 1 to 200, but the coordinates of the point cloud data belonging to the conveyor belt are roughly within 150 to 200, then the threshold range of the point cloud data coordinates belonging to the coal pile can be set to 1 to 150.

[0137] Step 2.2: Use a statistical filtering method based on statistical principles to remove large-scale noise in the three-dimensional point cloud data after direct filtering to obtain three-dimensional point cloud data after statistical filtering.

[0138] Step 2.2.1: For any point P in the three-dimensional point cloud data after direct filtering n , calculate point P n With point P n The K nearest neighbor average distance and standard deviation of the K nearest neighbor points in the neighborhood; where n is the nth point in the 3D point cloud data, and n=1,...,m; m represents the number of points in the 3D point cloud data.

[0139] The K nearest neighbors average distance is expressed as:

[0140]

[0141] where μ n Indicates point P n With point P n The average distance of the K nearest neighbors in the neighborhood; i is the i-th nearest neighbor; d i Indicates point P n With point P n The Euclidean distance between the i-th nearest neighbor in the neighborhood.

[0142] In this embodiment, the main idea of ​​the statistical filtering method based on statistical principles is to assume that the distance between points in the point cloud satisfies the Gaussian distribution, and its linearity depends on the mean μ and the standard deviation σ. Then, a distance threshold can be determined based on the mean and the standard deviation, and the distance threshold is controlled by a multiple of the standard deviation for filtering; when the average distance between a point and its nearest neighboring points is greater than the distance threshold, the point is judged as an outlier and removed.

[0143] The standard deviation is expressed as:

[0144]

[0145] where σ n is the standard deviation.

[0146] Step 2.2.2: Calculate the mean value of the distance between points and the variance of the distance between points in the three-dimensional point cloud data after straight-through filtering using the average distance and standard deviation of the K nearest neighbors of each point in the three-dimensional point cloud data after straight-through filtering and the K nearest neighbors in the neighborhood of the point.

[0147] The average value of the distance between the points is expressed as:

[0148]

[0149] Where U is the average distance between points.

[0150] The variance of the distance between the points is expressed as:

[0151]

[0152] Where δ is the variance of the distance between points.

[0153] Step 2.2.3: Calculate the distance threshold according to the average value of the distance between points and the variance of the distance between points in the three-dimensional point cloud data after the straight-through filtering.

[0154] The distance threshold is expressed as:

[0155] S=U±s td δ (5)

[0156] Where S is the distance threshold; s td is the threshold control coefficient.

[0157] Step 2.2.4: Determine the average distance μ between each point in the three-dimensional point cloud data after direct filtering and the K nearest neighbor points in the neighborhood of the point n Is it greater than the distance threshold S? If so, remove the point; otherwise, retain the point to obtain the three-dimensional point cloud data after statistical filtering.

[0158] Step 2.3: Use the moving least squares method to remove small-scale noise in the three-dimensional point cloud data after statistical filtering to obtain denoised three-dimensional point cloud data.

[0159] Step 2.3.1: Construct the fitting function.

[0160] The fitting function is expressed as:

[0161]

[0162] Where f(u) is the fitting function; u is the coordinate of the fixed point; a(u) is the coefficient of the fitting function to be determined, and a(u) = [a1(u)a2(u)...a l (u)]; p(u) is the basis function, and p(u)=[p1(u)p2(u)...p l (u)] T ; l is the number of basis functions; a q (u) is the coefficient of the qth term to be determined; p q (u) is the qth basis function.

[0163] Step 2.3.2: Select the cubic spline function as the basis function of the fitting function;

[0164] The cubic spline function is:

[0165]

[0166] Where ω(R) represents the cubic spline function, R = |r| / βh i , where r = uu i ,u i represents the coordinates of the i-th neighboring point in the neighborhood of the fixed point u; β is a constant, and its value range is usually 1.5 to 2.5. In this implementation, β is 2; h i is the support radius of the i-th neighbor point.

[0167] Step 2.3.3: Perform least squares fitting on the neighboring points in the neighborhood of the fixed point u, and smooth the three-dimensional point cloud data after statistical filtering by calculating the local curvature and normal vector of the fixed point to obtain the denoised three-dimensional point cloud data.

[0168] In this embodiment, the local curvature and normal vector of each point are estimated by performing least squares fitting on the neighborhood around each point, thereby eliminating noise points and outliers and smoothing the data. Figure 3 As shown in Figure 2, the effects before and after applying the point cloud filtering algorithm are compared.

[0169] Step 3: Use an adaptive point cloud simplification algorithm to simplify the denoised three-dimensional point cloud data to obtain simplified three-dimensional point cloud data.

[0170] This embodiment proposes an adaptive point cloud simplification algorithm for processing denoised three-dimensional point cloud data.

[0171] Step 3.1: For any point in the denoised 3D point cloud data, take the point as the target point, use KD-Tree to perform K nearest neighbor search on the target point, and obtain the K nearest neighbor points of the target point.

[0172] In this embodiment, the K nearest neighbor points obtained by performing K nearest neighbor search on each point in the denoised three-dimensional point cloud data in step 3.1 are essentially the same as the K nearest neighbor points in the neighborhood of each point in step 2, both of which find K points near the point, but for the same point, the nearest neighbor points determined in the above two steps are not necessarily the same point.

[0173] Step 3.2: Use the least squares method to fit the target point and its K nearest neighboring points, calculate the curvature of the target point, and obtain the curvature of all points in the denoised three-dimensional point cloud data.

[0174] Step 3.2.1: Construct the surface parametric equations.

[0175] The surface parameter equation is expressed as:

[0176] Q:t=t(α,β) (8)

[0177] Where Q represents the surface; t is a variable; (α, β) represents the coordinates of any point on the surface, and α and β are both variables.

[0178] Then the first basic formula of the surface is:

[0179]

[0180] Among them, G, E and F are the first-class basic parameters of the surface Q; t β is the partial derivative of the surface Q with respect to β; t α is the partial derivative of the surface Q with respect to α.

[0181] Step 3.2.2: For any point (α0, β0) on the surface Q, calculate the distance between the point (α0, β0) and its K nearest neighbors.

[0182] The distance between the point (α0, β0) and its K nearest neighbor points is expressed as:

[0183]

[0184] in is the distance between the point (α0, β0) and its K nearest neighboring points; Δα and Δβ are the increments of α and β respectively; h is the normal vector of the point (α0, β0); t(α0+Δα, β0+Δβ) is the surface equation at the point (α0+Δα, β0+Δβ); t(α0, β0) is the surface equation at the point (α0, β0).

[0185] Expand the right side of equation (10) using Taylor's formula:

[0186]

[0187] t α (α0, β0) is the value of the first-order partial derivative of the surface parameter equation with respect to α at the point (α0, β0); t β (α0, β0) is the value of the first-order partial derivative of the surface parameter equation with respect to β at the point (α0, β0); t αβ (α0, β0) is the value of the second-order partial derivative of the surface parameter equation with respect to α and β at the point (α0, β0).

[0188] and Therefore, by rewriting formula (11), we get the second basic formula of the surface, which is expressed as:

[0189]

[0190] Where L, M and N are the second type of basic parameters of the surface Q; t αα (α0, β0) is the value of the second-order partial derivative of the surface parameter equation with respect to α at the point (α0, β0); t ββ (α0, β0) is the value of the second-order partial derivative of the surface parameter equation with respect to β at the point (α0, β0).

[0191] In this embodiment, the second basic formula of the surface can be used to directly calculate the distance between a point (α0, β0) and its K neighboring points.

[0192] Step 3.2.3: Calculate the mean curvature of the surface.

[0193] The average curvature of the surface is expressed as:

[0194]

[0195] Where H represents the mean curvature of the surface.

[0196] Step 3.2.4: Determine the curvature of the target point based on the average curvature of the surface.

[0197] In this embodiment, the average curvature of the curved surface is approximated as the curvature of the target point.

[0198] Step 3.2.5: Repeat steps 3.2.1-3.2.4 to obtain the curvature of all points in the denoised 3D point cloud data.

[0199] Step 3.3: Set the curvature threshold, and use the improved octree voxel filtering algorithm to simplify the point cloud of all points whose curvature is greater than the curvature threshold in the denoised three-dimensional point cloud data to obtain simplified three-dimensional point cloud data.

[0200] In this embodiment, if Figure 4As shown in Figure 2, all points that meet the curvature requirements are selected from the denoised 3D point cloud data, that is, the curvature threshold H is set. st , if the curvature of the target point is less than the curvature threshold H st , then no point cloud simplification operation is performed. If it is greater than the curvature threshold H st The improved octree voxel filtering algorithm is used to perform point cloud simplification operations. The octree voxel filtering algorithm is an algorithm for point cloud processing. It divides the point cloud into a series of voxels, that is, small cubes in three-dimensional space, to achieve functions such as space partitioning and fast neighbor search. In the improved octree voxel filtering algorithm, a cube is constructed, and all points in the sub-cube of the cube are replaced by the point cloud data points closest to the center of gravity in the cube space. However, since all points in the space are not the same feature, the points that may be the distance from the center of gravity of another feature are also retained to fully retain the features of the original point cloud data, thereby obtaining smoothed point cloud data and completing the point cloud simplification operation.

[0201] Step 3.3.1: Set the curvature threshold and filter out all points whose curvature is greater than the curvature threshold from the denoised 3D point cloud data.

[0202] Step 3.3.2: Construct the first cube based on the maximum size of the point cloud space containing all the filtered points. The length, width, and height of the cube are recorded as X 长 , Y 宽 , Z 高 , the volume is V.

[0203] Step 3.3.3: Divide the first cube into C sub-cubes, calculate the centroid of each sub-cube, and mark the centroid coordinates as (x0, y0, z0); where x0, y0, and z0 are the horizontal, vertical, and vertical coordinates of the centroid, respectively.

[0204] The side length of the sub-cube is expressed as:

[0205]

[0206] Where L is the side length of the subcube; χ is the scale coefficient; ε is the scale factor; and N is the number of points in the point cloud space.

[0207]

[0208] The ceil() function is used to round up to an integer.

[0209] Step 3.3.4: Calculate the distance between each point in each sub-cube and the center of gravity of the sub-cube respectively, and only retain the point closest to the center of gravity and the point farthest from the center of gravity to obtain the simplified three-dimensional point cloud data.

[0210] In this embodiment, the distance l from each point in each sub-cube to the centroid of the sub-cube is calculated respectively. γ , γ represents the γth point in the subcube, and finds the closest distance l to the centroid min and the farthest distance l max , respectively corresponding to the point closest to the center of gravity (x min ,y min ,z min ) and the point farthest from the centroid (x max ,y max ,z max ), use point (x min ,y min ,z min ) and point (x max ,y max ,z max ) replaces all points in the sub-cube, that is, only the two points closest to and farthest from the center of gravity are retained to achieve the point cloud simplification and obtain the simplified 3D point cloud data. Among them:

[0211]

[0212] where η represents the number of points in the subcube, and represents the volume of each point, and

[0213] Step 4: Use the improved greedy triangulation algorithm to construct the triangular network topology structure between the points in the simplified 3D point cloud data.

[0214] Step 4.1: Perform a planar projection on the simplified three-dimensional point cloud data to obtain a point set on a two-dimensional plane.

[0215] Step 4.1.1: For any point p in the simplified 3D point cloud data, construct the tangent plane equation of the k-nearest neighbor points of point p based on the tangent plane normal vector of the surface fitted by point p and its k-nearest neighbor points at point p.

[0216] The tangent plane equation is:

[0217] ax+by+cz+e=0 (17)

[0218] Where a, b and c represent the components of the normal vector of the tangent plane in the directions of the X-axis, Y-axis and Z-axis respectively; e is a constant; x, y and z are the horizontal coordinate, vertical coordinate and vertical coordinate of any point in space respectively;

[0219] Step 4.1.2: Normalize the tangent plane normal vector and calculate the i-th nearest neighbor point p of point p i Distance to the tangent plane;

[0220] The i-th neighbor point p of the point p i The distance to the tangent plane is expressed as:

[0221] j i =|ax i +bx i +cz i +e| (18)

[0222] where j i Represents the i-th neighbor point p of point p i Distance to the tangent plane; (x i ,y i ,z i ) represents the neighboring point p i The coordinates of .

[0223] In this embodiment, the unitization operation of the tangent plane normal vector includes calculating the modulus of the normal vector, and then dividing each component of the normal vector by the modulus to obtain the unit normal vector. i where i=1,2,…,k.

[0224] Step 4.1.3: Based on the least squares principle, calculate the sum of the squares of the distances from the k nearest neighbor points of point p to the tangent plane. When the sum of the squares is the smallest, solve the tangent plane equation of the k-nearest neighbor points of point p.

[0225] The sum of the squares of the distances of the k nearest neighbors of the point p to the tangent plane is expressed as:

[0226]

[0227] Where D d Represents the sum of the squares of the distances from the k nearest neighbors of point p to the tangent plane.

[0228] In this embodiment, when When it is minimum, the tangent plane of the k-nearest neighbors of point p can be obtained. Therefore, the Lagrange multiplier method is used to calculate the extreme value of the sum of the squares of the distances from the k nearest neighbors of point p to the tangent plane.

[0229] The extreme value function constructed based on the Lagrange multiplier method is:

[0230]

[0231] Where T is the extreme value; λ is the Lagrange coefficient.

[0232] According to formula (20), i Find the partial derivative and get the partial derivative as:

[0233]

[0234] Let the partial derivative be zero, and combine formula (18) and (21) to obtain:

[0235]

[0236] in and They represent the mean of the horizontal coordinate, the mean of the vertical coordinate, and the mean of the vertical coordinate of the k neighboring points of point p respectively.

[0237] Substituting formula (22) into formula (20), we can calculate the partial derivatives of a, b and c respectively, and set the partial derivatives to zero to obtain:

[0238]

[0239] Simplifying formula (23) yields:

[0240]

[0241] And transform formula (24) into matrix form.

[0242]

[0243] in

[0244] Step 4.1.4: According to the solved tangent plane equation, project the k neighboring points of point p onto the tangent plane to obtain a point set on the two-dimensional plane.

[0245] The points obtained by projecting the k neighboring points of the point p onto the tangent plane are expressed as:

[0246]

[0247] in represents the coordinates of the projected point; p i Represents the coordinates before projection, and i represents the i-th point.

[0248] Step 4.2: Use the Delaunay triangulation algorithm based on the traditional point-by-point insertion method to triangulate the point set on the two-dimensional plane.

[0249] Step 4.2.1: Construct a super triangle containing all the points in the point set on the two-dimensional plane, and determine the extreme points x of the point set on the two-dimensional plane in the positive and negative directions of the coordinate axis of the two-dimensional coordinate system. max1 , x min1 ,y max1 and min1 , and add a positive offset value offset, and offset>0, to construct the minimum bounding box of the two-dimensional point set; where x max1is the extreme point of the point set on the two-dimensional plane in the positive direction of the X axis; min1 is the extreme point of the point set on the two-dimensional plane in the negative direction of the X axis; y max1 is the extreme point of the point set on the two-dimensional plane in the positive direction of the Y axis; min1 It is the extreme point of the point set on the two-dimensional plane in the negative direction of the Y axis.

[0250] In this embodiment, if Figure 5 As shown in the figure, the three-dimensional point cloud data is transformed into a point set in a two-dimensional coordinate system through a plane projection transformation, and the extreme point x of the point set in the positive and negative directions of the coordinate axis is found. max1 , x min1 ,y max1 and min1 , and on this basis add a positive offset value offset, and offset>0, to construct the minimum bounding box of the two-dimensional point set.

[0251]

[0252] Where (X max1 ,Y max1 ), (X max1 ,Y min1 ), (X min1 ,Y max1 ) and (X min1 ,Y min1 ) are the coordinates of the four vertices of the minimum bounding box.

[0253] In this embodiment, a large triangle, ie, a super triangle, is constructed which includes all points in the point set on the two-dimensional plane, and it is ensured that all points are inside the triangle.

[0254] Step 4.2.2: Traverse all points in the point set on the two-dimensional plane, and insert each point into an initial triangular grid in turn to obtain a new triangular grid.

[0255] The method for inserting each point into an initial triangular grid is as follows: for the currently traversed point, a triangle containing the point is determined according to the coordinates of the point; if the point is located inside the triangle, the three vertices of the triangle are connected to the point respectively, thereby dividing the triangle into three new triangles; if the point is located on the edge of the triangle, the vertices of the two triangles adjacent to the edge are connected respectively, thereby dividing the two triangles adjacent to the edge into four new triangles.

[0256] In this embodiment, if Figure 6 and Figure 7As shown, all points in the two-dimensional point set are traversed in storage order and each point is inserted into the initial triangular grid in turn. For each newly inserted point, determine which triangle the point is in based on the coordinates of the point, and find the triangle containing the point; if the point is inside a triangle, the triangle will be split into three new triangles by the current point; if the point is on an edge, the two triangles connected by the edge will be split into four new triangles.

[0257] Step 4.2.3: Perform local optimization on the new triangular grid to obtain an optimized triangular grid.

[0258] The local optimization method is: for all new triangles in the new triangular grid, determine whether the edge of each new triangle satisfies the Delaunay property, if so, keep it; if not, flip the edge that does not satisfy the Delaunay property.

[0259] In this embodiment, in order to maintain the properties of Delaunay triangulation, all new triangles in the initial triangular mesh are checked after all points are inserted, and the edges that do not meet the Delaunay property are flipped.

[0260] Step 4.2.4: Delete the edges and vertices of the super triangle in the optimized triangular mesh to obtain the triangular network topology structure between the points in the simplified three-dimensional point cloud data.

[0261] In this embodiment, if Figure 8 As shown, after all points are inserted into the initial triangle mesh and all triangles are locally optimized, the initially added super triangles and their associated edges and vertices are deleted.

[0262] Step 5: Use the hole identification algorithm to identify all holes in the triangular network topology structure.

[0263] In this embodiment, the three-dimensional point cloud data of the coal yard will have holes due to occlusion, uneven point cloud density, etc. Even if the improved greedy triangulation algorithm is used for surface reconstruction, there will still be many holes. Therefore, it is necessary to use a hole recognition algorithm to identify all the holes and classify the holes.

[0264] Step 5.1: Obtain the edges of all triangles in the triangular network topology and the face neighbor relationship of each edge, and divide all the obtained edges into mesh internal edges and mesh boundary edges according to the face neighbor relationship of each edge.

[0265] In this embodiment, the main goal of triangular mesh hole recognition is to identify and fill holes in triangular meshes. Hole recognition can usually be achieved by detecting the boundaries of the mesh. First, an edge data structure is constructed to extract all edges from the triangular network topology and record the face neighbor relationship of each edge. According to the face neighbor relationship of each edge, all acquired edges are divided into mesh internal edges and mesh boundary edges. In general, the edges connecting two triangular facets are called mesh internal edges, and the edges connected to only one triangular facet are called mesh boundary edges. Since mesh boundary edges are edges that belong to only one triangle, these boundary edges can be found by traversing all edges. A depth-first search (DFS) or breadth-first search (BFS) algorithm can be used to construct one or more boundary rings starting from the found boundary edges.

[0266] Step 5.2: Call the PolygonMesh data structure to obtain the connection relationship between the vertices and edges of all triangles in the triangular network topology structure, traverse all vertices and filter out the vertices that are not referenced by any triangles, and identify the holes in the triangular network topology structure based on the connection relationship between the filtered vertices and edges.

[0267] In this embodiment, if Figure 7 As shown in the figure, when realizing hole recognition, the key is to identify those vertices that are not covered by any triangles. These vertices are usually located inside or outside the model and are not associated with any edges or faces. These vertices can be regarded as potential hole vertices. In order to identify these vertices, it is necessary to traverse and analyze the data structure of the triangular mesh model. In order to realize hole boundary recognition, the PolygonMesh data structure in PCL is used to store the triangular face information in the triangular network topology structure. The PolygonMesh data structure is an existing structure in the widely used open source library PCL (PointCloudLibrary). The data structure contains the following lists: a list for storing vertex coordinates; a list for storing indexes of points constituting triangular facets and a list for storing information such as normal vectors. By calling the information of vertices and edges in the list contained in the data structure, the index list of points for the boundary recognition of the triangular mesh model and the list for storing information such as normal vectors are realized, that is, by calling the information of vertices and edges in the list contained in the data structure, the boundary recognition of the triangular mesh model is realized.

[0268] In this embodiment, all vertices are traversed to check whether each vertex is referenced by any triangle. If a vertex is not referenced by any triangle, then the vertex may be inside or outside the hole. And by analyzing the connection relationship of the edges, it can be further determined whether these unreferenced vertices actually form a hole. For example, it can be checked whether there is one or more unreferenced vertices surrounding a closed area, and if so, the closed area is used as an identified hole. Through the above process, holes in the triangular network topology can be effectively identified. Common holes include internal holes, non-structural boundary holes, and structural boundary holes.

[0269] Step 6: Use the minimum angle mesh filling algorithm to add patch surfaces to all identified holes in the triangulated network topology.

[0270] Step 6.1: Get the coordinates of the mesh boundary edges and boundary vertices of the hole, calculate the Euclidean distance between two boundary vertices in each mesh boundary edge in space to obtain the length of the mesh boundary edge, and calculate the average length lavg of all mesh boundary edges.

[0271] Step 6.2: For each boundary vertex of a mesh boundary edge, calculate the adjacent edge angle between the two adjacent edges corresponding to the boundary vertex.

[0272] Step 6.3: Filter out the boundary vertex Pi with the smallest adjacent edge angle value, and calculate the distance D between the two boundary vertices adjacent to the boundary vertex.

[0273] Step 6.4: Add triangles to the triangulated network topology according to the distance D.

[0274] In this embodiment, if Fig. 9 As shown, the method for adding a triangular face in the triangular network topology structure according to the distance D is: determine whether D<1.2×lavg is established; if established, add a triangular face by connecting the above three boundary vertices; if not established, when adding a triangular face, the edges of the two boundary vertices adjacent to the boundary vertex with the smallest angle value of the adjacent edge need to be connected, such as Fig.10 As shown, a new boundary vertex Pnew is added at the midpoint of the edge Pi-1Pi+1, and the elevation of the boundary vertex Pnew is also the median of the elevations of points Pi-1 and Pi+1. Two triangles are added by connecting points Pnew, Pi-1 and Pi+1; Pi-1 and Pi+1 are two boundary vertices adjacent to the boundary vertex with the smallest angle value with the adjacent edge.

[0275] Step 6.5: Update the coordinates of the boundary vertices of the hole and delete the boundary vertex Pi that has been repaired.

[0276] Step 6.6: Determine whether there are any unrepaired boundary vertices. If so, repeat steps 6.2-6.5; otherwise, the hole filling is completed.

[0277] In this embodiment, it is determined whether the hole is completely filled. If not, steps 6.2-6.5 are repeated.

[0278] Step 7: Use a mesh subdivision algorithm to subdivide the added patch surface and the original surface in the triangular network topology structure.

[0279] Step 7.1: Add a new vertex P' to the triangular network topology structure after the patch surface is added, and calculate the position of the new vertex P' according to the positions of all vertices in the triangular network topology structure after the patch surface is added.

[0280] In this embodiment, if Fig.11 As shown, first, a new vertex P' is added to the original triangle mesh, that is, the triangle network topology structure after adding the patch surface, and according to the vertex position of the neighborhood of the original triangle mesh, that is, P1+P2+…+P o Calculate the position of the new vertex P', expressed as:

[0281] P'=(1-oψ)P+ψ(P1+P2+…+P o ) (28)

[0282] Where P' is the new vertex; P represents the vertex; o is the degree of vertex P; ψ is the coefficient.

[0283] The coefficient ψ is defined by the following formula:

[0284]

[0285] Step 7.2: Add new edge points to the triangulated network topology after adding the patch surface.

[0286] In this embodiment, each edge in the triangular network topology structure after adding the patch surface is solved in turn to obtain a new edge point. Taking edge P1P2 as an example, the solution formula for the new edge point is as follows:

[0287]

[0288] Where P" is the new edge point; P3P4 are the other two vertices of the two triangles with P1P2 as the common edge, except P1P2;

[0289] Step 7.3: For each new vertex P', connect point P' to all new vertices adjacent to point P' to obtain new edges; for each triangle, connect the three new edge points in the triangle to form a new edge.

[0290] Step 7.4: Every three new edge points and new vertices connected to each other are combined into a new triangle to complete the subdivision of the triangular network topology structure.

[0291] In this embodiment, the first two steps of mesh subdivision are called geometric rules, which define the spatial position of the new point and are represented by B. The last two steps of mesh subdivision are called topological rules, which define the connection relationship between the new points and are represented by C. The process of the mesh subdivision method can be understood as follows:

[0292] 1B→1C→2B→2C→3B→3C→... (31)

[0293] Among them, 1B and 1C represent the rules applied to the triangular network topology structure before subdivision; 2B and 2C represent the rules applied to the triangular network topology structure obtained after the first subdivision; 3B and 3C represent the rules applied to the triangular network topology structure obtained after the second subdivision. For the triangular network topology structure before mesh refinement, the mesh density at the patch surface is less than that of the original triangular network topology structure, and the number of operable vertices is small. After mesh refinement, the number of vertices in the triangular network topology structure with the patch surface added doubles, and the number of points that can be meshed smoothed increases. In addition, the mesh of the original triangular network topology structure is also meshed subdivided, and the number of vertices of the original triangular network topology structure is also increased many times. In this way, after the surface at the patch is smoothed, the edge is easier to merge with the original model at the hole boundary.

[0294] Step 8: Use a mesh smoothing algorithm to smooth the triangular network topology structure after subdivision to obtain the repaired and reconstructed three-dimensional point cloud data.

[0295] Step 8.1: Establish the energy function.

[0296] In this embodiment, it is desired to find a function f(x) such that the second-order derivative is minimum at all vertices in the patch, so an energy function is established as follows:

[0297]

[0298] Where E(f) represents the energy function; (v,w) represents an integration interval; f ss represents the second-order derivative; s is the independent variable.

[0299] The above integral can be viewed as the integral of each vertex (f ss ) 2 If this integral is small, then the second-order derivative in the interval (v, w) will be small, which means that the function f meets our requirements and the created patch is a good patch surface.

[0300] Step 8.2: Use variational calculus to solve the functional function f that minimizes the energy function;

[0301] Step 8.2.1: Assume that f is the minimum value of the energy function E(f), and set is any function in the interval (v,w), and Satisfies first-order differentiability in the interval (v,w) and satisfies

[0302] Step 8.2.2: Assume there exists a function Where ζ is the coefficient.

[0303] In this embodiment, the function It is based on the independent variable f(s) of the energy function and adds and The size of depends on the coefficient ζ.

[0304] Step 8.2.3: Assume that f(s) is the minimum value of the energy function E(f). When ζ = 0, the energy function reaches its minimum value, and the derivative of the energy function must also be zero, which can be expressed as:

[0305]

[0306] For convenience, the formula is transformed as follows:

[0307]

[0308] in express f represents f(s).

[0309] Step 8.2.4: Calculate the derivative with respect to ζ;

[0310] The derivative with respect to ζ is expressed as:

[0311]

[0312] in express The second derivative of ss represents the second-order derivative of f.

[0313] Step 8.2.5: When ζ = 0, the right side of the equation must satisfy the condition of being equal to zero in the derivative with respect to ζ;

[0314] Because you need to follow:

[0315]

[0316] Step 8.2.6: Integrate formula (36) step by step to obtain:

[0317]

[0318] where f sss represents the third-order derivative of f(s), express

[0319] Step 8.2.7: Combine Simplifying formula (37) yields:

[0320]

[0321] Step 8.2.8: Integrate formula (38) step by step to obtain:

[0322]

[0323] where f ssss represents the fourth-order derivative of f(s).

[0324] Step 8.2.9: According to Simplifying formula (39) yields:

[0325]

[0326] Step 8.2.10: When formula (38) is satisfied, And the only case where the function is differentiable is:

[0327] f ssss =0 (41)

[0328] Step 8.2.11: Equivalently transform formula (39) into formula (41) to obtain:

[0329] ΔΔf=Δ 2 f=0 (42)

[0330] where Δ 2 It is often called the bi-Laplacian or biharmonic operator.

[0331] In this embodiment, since f ssss It must always be zero in the domain interval, which means that when E(f) reaches its minimum value, its fourth-order derivative is always zero in the domain. According to this rule, the minimum value f of E(f) is derived to minimize the curvature change, thus making the plane transition smoother. When the fourth-order derivative of function f is zero, its Laplace operator must be zero. Therefore, we can solve it by solving the Laplace operator. In the one-dimensional Laplace operator, there exists Δf=f ss .

[0332] Step 8.2.12: For the refined patch surface in the triangular network topology structure, the Laplace-Beltrami algorithm in the libigl library is called to solve formula (40) to obtain the smoothed triangular network topology structure and complete the repair and reconstruction of the three-dimensional point cloud data.

[0333] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some or all of the technical features therein. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope defined by the claims of the present invention.

Claims

1. An image reconstruction method for incomplete three-dimensional point cloud data of coal in a coal field, characterized in that: The method comprises the following steps: Step 1: Collect initial 3D point cloud data of the coal pile using LiDAR; Step 2: Use a point cloud filtering algorithm to remove noise from the initial three-dimensional point cloud data to obtain denoised three-dimensional point cloud data; Step 3: Adopt an adaptive point cloud simplification algorithm to simplify the denoised three-dimensional point cloud data to obtain simplified three-dimensional point cloud data; Step 4: Use the improved greedy triangulation algorithm to construct the triangular network topology structure between the points in the simplified three-dimensional point cloud data; Step 5: Use hole identification algorithm to identify all holes in the triangular network topology structure; Step 6: Use the minimum angle mesh filling algorithm to add patch surfaces to all identified holes in the triangulated network topology; Step 7: Use a mesh subdivision algorithm to subdivide the added patch surface and the original surface in the triangular network topology structure; Step 8: Use a mesh smoothing algorithm to smooth the triangular network topology structure after subdivision processing to obtain the repaired and reconstructed three-dimensional point cloud data.

2. According to claim 1, the image reconstruction method for incomplete three-dimensional point cloud data of coal in a coal field is characterized in that: The step 2 further comprises: Step 2.1: Use straight-through filtering to remove invalid point cloud data in the initial three-dimensional point cloud data to obtain three-dimensional point cloud data after straight-through filtering; Step 2.2: Using a statistical filtering method based on statistical principles to remove large-scale noise in the three-dimensional point cloud data after direct filtering, to obtain three-dimensional point cloud data after statistical filtering; Step 2.2.1: For any point P in the three-dimensional point cloud data after direct filtering n , calculate point P n With point P n The K nearest neighbors average distance and standard deviation of the K nearest neighbors in the neighborhood; where n is the nth point in the 3D point cloud data, and n=1,...,m; m represents the number of points in the 3D point cloud data; Step 2.2.2: Calculate the average value of the distance between points and the variance of the distance between points in the three-dimensional point cloud data after the straight-through filtering by using the average distance and standard deviation of the K nearest neighbors of each point in the three-dimensional point cloud data after the straight-through filtering and the K nearest neighbors in the neighborhood of the point; Step 2.2.3: Calculate the distance threshold according to the average value of the distance between points and the variance of the distance between points in the three-dimensional point cloud data after the straight-through filtering; Step 2.2.4: Determine the average distance μ between each point in the three-dimensional point cloud data after direct filtering and the K nearest neighbor points in the neighborhood of the point n Is it greater than the distance threshold S? If so, remove the point; otherwise, keep the point to obtain the three-dimensional point cloud data after statistical filtering; Step 2.3: Use the moving least squares method to remove small-scale noise in the three-dimensional point cloud data after statistical filtering to obtain denoised three-dimensional point cloud data; Step 2.3.1: Construct the fitting function; Step 2.3.2: Select the cubic spline function as the basis function of the fitting function; Step 2.3.3: Perform least squares fitting on the neighboring points in the neighborhood of the fixed point u, and smooth the three-dimensional point cloud data after statistical filtering by calculating the local curvature and normal vector of the fixed point to obtain the denoised three-dimensional point cloud data.

3. According to claim 2, the image reconstruction method for incomplete three-dimensional point cloud data of coal in a coal field is characterized in that: The step 3 further comprises: Step 3.1: For any point in the denoised 3D point cloud data, take the point as the target point, use KD-Tree to perform K nearest neighbor search on the target point, and obtain K nearest neighbor points of the target point; Step 3.2: Use the least squares method to fit the target point and its K neighboring points, calculate the curvature of the target point, and obtain the curvature of all points in the denoised three-dimensional point cloud data; Step 3.3: Set the curvature threshold, and use the improved octree voxel filtering algorithm to simplify the point cloud of all points whose curvature is greater than the curvature threshold in the denoised three-dimensional point cloud data to obtain simplified three-dimensional point cloud data.

4. According to claim 3, the image reconstruction method for incomplete three-dimensional point cloud data of coal in a coal field is characterized in that: The step 3.2 further comprises: Step 3.2.1: Construct the surface parametric equation; Step 3.2.2: For any point (α0, β0) on the surface Q, construct the second basic formula of the surface, and then calculate the distance between the point (α0, β0) and its K nearest neighbors; The second basic formula of the surface is expressed as: Where L, M and N are the second-class basic parameters of the surface Q; α and β are variables; t αα (α0, β0) is the value of the second-order partial derivative of the surface parameter equation with respect to α at the point (α0, β0); t ββ (α0, β0) is the value of the second-order partial derivative of the surface parameter equation with respect to β at the point (α0, β0); t αβ (α0, β0) is the value of the second-order partial derivative of the surface parameter equation with respect to α and β at the point (α0, β0); h is the normal vector of the point (α0, β0); is the distance between the point (α0, β0) and its K nearest neighbors; Δα and Δβ are the increments of α and β respectively; o(Δα 2 ,Δβ 2 ) is an infinitesimal term; Step 3.2.3: Calculate the mean curvature of the surface; The average curvature of the surface is expressed as: Where H represents the average curvature of the surface; G, E and F are the first-class basic parameters of the surface Q; Step 3.2.4: Determine the curvature of the target point based on the average curvature of the surface; Step 3.2.5: Repeat steps 3.2.1-3.2.4 to obtain the curvature of all points in the denoised 3D point cloud data.

5. The image reconstruction method for incomplete three-dimensional point cloud data of coal in a coal field according to claim 4, characterized in that: The step 3.3 further comprises: Step 3.3.1: Set the curvature threshold and filter out all points whose curvature is greater than the curvature threshold from the denoised 3D point cloud data; Step 3.3.2: Construct the first cube based on the maximum size of the point cloud space containing all the filtered points. The length, width, and height of the cube are recorded as X 长 , Y 宽 , Z 高 , the volume is V; Step 3.3.3: Divide the first cube into C sub-cubes, calculate the centroid of each sub-cube, and mark the centroid coordinates as (x0, y0, z0); where x0, y0, and z0 are the horizontal, vertical, and vertical coordinates of the centroid respectively; Step 3.3.4: Calculate the distance between each point in each sub-cube and the center of gravity of the sub-cube, and only keep the point closest to the center of gravity and the point farthest from the center of gravity to obtain the simplified three-dimensional point cloud data; The coordinates of the point closest to the center of gravity and the point farthest from the center of gravity are expressed as: where x max The horizontal coordinate of the point farthest from the center of gravity; y max The ordinate of the point farthest from the center of gravity; z max represents the vertical coordinate of the point farthest from the centroid; η represents the number of points in the subcube; γ represents the γth point in the subcube; represents the volume of each point; x γ represents the horizontal coordinate of the γth point; l max Indicates the farthest distance from the center of gravity; y γ represents the ordinate of the γth point; z γ represents the vertical coordinate of the γth point; x min The horizontal coordinate of the point closest to the center of gravity; y min The ordinate of the point closest to the center of gravity; z min Indicates the vertical coordinate of the point closest to the center of gravity; l min Indicates the closest distance to the center of gravity.

6. The image reconstruction method for incomplete three-dimensional point cloud data of coal in a coal field according to claim 5, characterized in that: The step 4 further comprises: Step 4.1: Perform a planar projection on the simplified three-dimensional point cloud data to obtain a point set on a two-dimensional plane; Step 4.1.1: For any point p in the simplified 3D point cloud data, construct the tangent plane equation of the k-nearest neighboring points of point p according to the tangent plane normal vector of the surface fitted by point p and its k-nearest neighboring points at point p; Step 4.1.2: Normalize the tangent plane normal vector and calculate the i-th nearest neighbor point p of point p i Distance to the tangent plane; Step 4.1.3: Based on the least squares principle, calculate the sum of the squares of the distances from the k nearest neighbor points of point p to the tangent plane. When the sum of the squares is the smallest, solve the tangent plane equation of the k-nearest neighbor points of point p. Step 4.1.4: Project the k nearest neighbor points of point p onto the tangent plane according to the solved tangent plane equation to obtain a point set on the two-dimensional plane; Step 4.2: Use the Delaunay triangulation algorithm based on the traditional point-by-point insertion method to triangulate the point set on the two-dimensional plane to obtain the triangular network topology structure between each point in the simplified three-dimensional point cloud data; Step 4.2.1: Construct a super triangle containing all the points in the point set on the two-dimensional plane, and determine the extreme points x of the point set on the two-dimensional plane in the positive and negative directions of the coordinate axis of the two-dimensional coordinate system. max1 , x min1 ,y max1 and min1 , and add a positive offset value offset, and offset>0, to construct the minimum bounding box of the two-dimensional point set; where x max1 is the extreme point of the point set on the two-dimensional plane in the positive direction of the X axis; min1 is the extreme point of the point set on the two-dimensional plane in the negative direction of the X axis; y max1 is the extreme point of the point set on the two-dimensional plane in the positive direction of the Y axis; min1 It is the extreme point of the point set on the two-dimensional plane in the negative direction of the Y axis; Step 4.2.2: Traverse all points in the point set on the two-dimensional plane, and insert each point into an initial triangular grid in turn to obtain a new triangular grid; Step 4.2.3: Locally optimize the new triangular grid to obtain an optimized triangular grid; Step 4.2.4: Delete the edges and vertices of the super triangle in the optimized triangular mesh to obtain the triangular network topology structure between the points in the simplified three-dimensional point cloud data.

7. The image reconstruction method for incomplete three-dimensional point cloud data of coal in a coal field according to claim 6 is characterized in that: The step 5 further comprises: Step 5.1: Obtain the edges of all triangles in the triangular network topology and the face neighbor relationship of each edge, and divide all the obtained edges into mesh internal edges and mesh boundary edges according to the face neighbor relationship of each edge; Step 5.2: Call the PolygonMesh data structure to obtain the connection relationship between the vertices and edges of all triangles in the triangular network topology structure, traverse all vertices and filter out the vertices that are not referenced by any triangles, and identify the holes in the triangular network topology structure based on the connection relationship between the filtered vertices and edges.

8. The image reconstruction method for incomplete three-dimensional point cloud data of coal in a coal field according to claim 7, characterized in that: The step 6 further comprises: Step 6.1: Obtain the coordinates of the mesh boundary edges and boundary vertices of the hole, calculate the Euclidean distance between two boundary vertices in each mesh boundary edge in space to obtain the length of the mesh boundary edge, and calculate the average length of all mesh boundary edges; Step 6.2: For each vertex of the mesh boundary edge, calculate the adjacent edge angle between the two adjacent edges corresponding to the vertex; Step 6.3: Filter out the boundary vertex with the smallest adjacent edge angle value, and calculate the distance D between the two boundary vertices adjacent to the boundary vertex; Step 6.4: Add triangles to the triangular network topology structure according to the distance D; Step 6.5: Update the coordinates of the boundary vertices of the hole and delete the boundary vertex Pi that has been repaired; Step 6.6: Determine whether there are any unrepaired boundary vertices. If so, repeat steps 6.2-6.5; if not, the hole filling is completed.

9. The image reconstruction method for incomplete three-dimensional point cloud data of coal in a coal field according to claim 8, characterized in that: The step 7 further comprises: Step 7.1: Add a new vertex P' to the triangular network topology structure after the patch surface is added, and calculate the position of the new vertex P' according to the positions of all vertices in the triangular network topology structure after the patch surface is added; Step 7.2: Add new edge points to the triangular network topology structure after adding the patch surface; Step 7.3: For each new vertex P', connect point P' to all new vertices adjacent to point P' to obtain new edges; for each triangle, connect the three new edge points in the triangle to form a new edge; Step 7.4: Every three new edge points and new vertices connected to each other are combined into a new triangle to complete the subdivision of the triangular network topology structure.

10. The image reconstruction method for incomplete three-dimensional point cloud data of coal in a coal field according to claim 9, characterized in that: The step 8 further comprises: Step 8.1: Establish energy function; The energy function is expressed as: Where E(f) represents the energy function; (v,w) represents an integration interval; f ss represents the second-order derivative; s is the independent variable; Step 8.2: Use variational calculus to solve the functional function f that minimizes the energy function; Step 8.2.1: Assume that f is the minimum value of the energy function E(f), and set is any function in the interval (v,w), and Satisfies first-order differentiability in the interval (v,w) and satisfies Step 8.2.2: Assume there exists a function Where ζ is the coefficient; Step 8.2.3: Assume that f(s) is the minimum value of the energy function E(f). When ζ = 0, the energy function reaches its minimum value, and the derivative of the energy function must also be zero, which is recorded as: Formula (33) is transformed into: in express f represents f(s); Step 8.2.4: Calculate the derivative with respect to ζ; The derivative with respect to ζ is expressed as: in express The second derivative of ss represents the second-order derivative of f; Step 8.2.5: When ζ = 0, the right side of the equation must satisfy the condition of being equal to zero in the derivative with respect to ζ, which is written as: Step 8.2.6: Integrate formula (36) step by step to obtain: where f sss represents the third-order derivative of f(s), express Step 8.2.7: Combine Simplifying formula (37) yields: Step 8.2.8: Integrate formula (38) step by step to obtain: where f ssss represents the fourth-order derivative of f(s); Step 8.2.9: According to Simplifying formula (39) yields: Step 8.2.10: When formula (38) is satisfied, And the only case where the function is differentiable is: f ssss ＝0 (41) Step 8.2.11: Equivalently convert formula (39) into formula (41) to obtain: ΔΔf=Δ 2 f=0 (42) Where Δ 2 is a double Laplacian operator; Step 8.2.12: For the refined patch surface in the triangular network topology structure, the Laplace-Beltrami algorithm in the libigl library is called to solve formula (40) to obtain the smoothed triangular network topology structure and complete the repair and reconstruction of the three-dimensional point cloud data.

Images