High-precision method for calculating volume of silo based on multi-scale grid reconstruction and constraint protection

Abstract

The application provides a high-precision silo volume calculation method based on multi-scale grid reconstruction and constraint protection, and relates to the technical field of material volume measurement, and comprises the following steps: pre-processing a material surface point cloud; performing point cloud edge completion on the pre-processed material surface point cloud, so as to obtain edge points after completion as an accurate material surface boundary; splicing the segmented partial template point cloud and the material surface point cloud after edge completion to obtain spliced point cloud; performing convex hull calculation and grid reconstruction on the spliced point cloud, and obtaining a boundary edge and a boundary point of the spliced point cloud; updating the boundary edge and the boundary point based on the distance from the edge point to the grid edge; performing progressive optimization reconstruction on the grid; performing local convex hull approximation repair on the boundary edge and the boundary point of the reconstructed grid to obtain an optimized closed grid; and performing volume calculation based on the optimized closed grid; and the application solves the problems of lacking a multi-scale processing mechanism and lacking a boundary protection mechanism in the prior art.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of material volume measurement, and particularly relates to a high-precision silo volume calculation method based on multi-scale grid reconstruction and constraint protection. BACKGROUND

[0002] At present, the existing silo volume calculation methods include the following technologies:

[0003] 1. Traditional convex hull volume calculation technology

[0004] This technology calculates the convex hull directly from the original point cloud and uses the convex hull volume as the approximate volume. However, there are significant defects: first, the convex hull volume is much larger than the actual volume, and the error can reach 30-50% when there is a depression in the material surface, causing systematic overestimation; second, it cannot accurately process the depression area of the material surface, resulting in a deviation of the volume calculation result from the true value; third, it is extremely sensitive to boundary loss, and when the point cloud boundary is incomplete, the error increases sharply, which cannot meet the industrial precision requirements; finally, it lacks a boundary protection mechanism, and its reliability is poor in complex industrial environments.

[0005] 2. Grid reconstruction volume calculation technology

[0006] This technology uses Poisson reconstruction, Delaunay triangulation and other methods to reconstruct the surface grid and then calculates the volume. However, there are obvious shortcomings: first, holes and cracks are easily generated at the boundary, resulting in an unsealed grid, which directly affects the accuracy of volume calculation; second, it is sensitive to noise and outliers, requiring a large amount of manual repair and parameter tuning, increasing the complexity of use; third, the calculation complexity is high, and the real-time performance is poor when processing large-scale point clouds, making it difficult to meet the demand of rapid measurement in industrial sites; finally, it lacks a boundary constraint protection mechanism, and the boundary integrity is easily damaged during grid optimization.

[0007] 3. Slice integration volume calculation technology

[0008] This technology slices the point cloud, calculates the area of each layer, and then integrates to obtain the volume. However, there are key defects: first, the slice area calculation is wrong when the boundary is incomplete, resulting in systematic errors; second, it cannot handle complex topological structures, and the precision decreases sharply when there are overhangs or complex shapes in the material surface; third, the precision is significantly affected by the slice thickness, and it is difficult to balance precision and calculation efficiency; finally, it lacks an adaptive mechanism and cannot adjust the calculation strategy according to different point cloud densities and completeness.

[0009] 4. Deep learning surface reconstruction technology

[0010] The technology calculates the volume after predicting the complete surface through a neural network. However, there are fundamental limitations: first, a large amount of high-quality labeled data is needed for training, and the data acquisition cost is high; second, the model is over-fitted to a specific scene, and the generalization ability is poor in a new scene, which needs to be retrained; third, the demand for computing resources is large, and the inference process needs GPU support, which is not suitable for edge device deployment; and finally, the "black box" characteristics of the model result in a lack of explainability of the results, making it difficult to perform error analysis and optimization adjustment.

[0011] In summary, the existing technology generally has the following fundamental problems: the lack of boundary protection mechanism leads to boundary deformation or loss during grid optimization and deformation; and the lack of multi-scale processing mechanism cannot find a suitable balance between precision and efficiency. SUMMARY

[0012] The purpose of the present application is to solve the problems of lack of multi-scale processing mechanism and lack of boundary protection mechanism in the background art, and to provide a high-precision silo volume calculation method based on multi-scale grid reconstruction and constraint protection.

[0013] The technical scheme of the present application is as follows:

[0014] A high-precision silo volume calculation method based on multi-scale grid reconstruction and constraint protection, comprising the following steps:

[0015] S1. Preprocessing the material surface point cloud data;

[0016] S2. Completing the edges of the preprocessed material surface point cloud to obtain the accurate material surface boundary;

[0017] S3. Segmenting the template point cloud through the edge-completed material surface point cloud, and splicing the segmented template point cloud and the edge-completed material surface point cloud to obtain a spliced point cloud;

[0018] S4. Convex hull calculation and grid reconstruction are performed on the spliced point cloud, and the boundary edges and boundary points of the spliced point cloud are obtained;

[0019] S5. Updating the boundary edges and boundary points based on the distance from the edge points to the grid edges;

[0020] S6. Gradually optimizing and reconstructing the grid;

[0021] S7. Repairing the local convex hull approximation in the influence range of the boundary edges and boundary points of the reconstructed grid to obtain an optimized closed grid;

[0022] S8. Volume calculation based on the optimized closed grid.

[0023] The overall scheme of the present application starts from extracting the material surface boundary, then obtains the boundary of the spliced point cloud, then obtains the accurate boundary of the spliced point cloud, and finally optimizes the grid based on the boundary. These are essentially to protect the accurate presentation of the model at the boundary, that is, boundary protection.

[0024] Specifically, the S2 comprises:

[0025] S2.1. Projecting the pre-processed material surface point cloud to the X-Y plane;

[0026] S2.2. In the X-Y plane, evenly divide a 360-degree circle into K sectors with the point cloud geometric center as the center;

[0027] S2.3. In each sector, traverse all points, calculate the polar coordinate angle and find the farthest edge point F i , whose distance from the center is r i , where i represents the i-th sector.

[0028] S2.4. For the edge point sector, adopt the nearest neighbor sector replacement completion strategy to complete the edge point, and add the completed edge point to the material surface point cloud and update the sector;

[0029] S2.5. Obtaining the accurate material surface boundary according to the edge points of all sectors.

[0030] Specifically, in the S2.4, for any edge point sector, the nearest neighbor sector replacement completion strategy to complete the edge point is specifically:

[0031] Traverse all sectors with edge points, calculate the angle difference between each sector and the current edge point sector, sort the near neighbor sequence in ascending order of angle difference, and the sector with the smallest angle difference in the near neighbor sequence is p, and the edge point is F p , the distance of the edge point from the center is r p , let r o equal to the radius r o of the current edge point sector, and calculate the edge point, the specific calculation process is as follows:

[0032] The x and y coordinates of the edge point of the current edge point sector are respectively,

[0033] x=x c +r o cosθ

[0034] y=y c +r o sinθ

[0035] Wherein, x c , y cθ is the angle of the current edgeless point sector, and

[0036] The z coordinate of the edge point of the current edgeless point sector is calculated by the nearest neighbor inverse distance weighted interpolation,

[0037] z = (d q / (d p +d q ))·z p + (d p / (d p +d q ))·z q .

[0038] wherein z p , z q are the z coordinates of the edge points F p , F q of the sectors p and q with the smallest and the second smallest angle difference in the neighbor sequence, and d p , d q are the distances from the X-Y plane to the edge points F p , F q .

[0039] Specifically, the S3 comprises:

[0040] S3.1. Construct a KD tree index of the material surface point cloud, and perform a nearest neighbor search on each point in the template point cloud to find the nearest material surface point and calculate the height difference value of the two points.

[0041] S3.2. Traverse all the template points, if the height of a template point is lower than that of its nearest material surface point, the template point is retained; if the height of a template point is higher than that of its nearest material surface point, but the height difference value does not exceed the height tolerance threshold ε, the template point is also retained, and the height of the template point is adjusted to the height of its nearest material surface point; the remaining template points other than the above two kinds of template points are not retained.

[0042] S3.3. Merge the material surface point cloud and the retained template points into a spliced point cloud.

[0043] Specifically, the S4 comprises:

[0044] S4.1. Calculate the initial convex hull mesh of the spliced point cloud.

[0045] S4.2. Identify boundary edges and boundary points of the initial convex hull mesh; traverse each mesh edge e, take half edge h and the two corresponding sides f1 and f2 of its opposite half edge, calculate the normal vectors n1 and n2 of the two sides, and consider the edge as a boundary edge when the angle between the normal vectors n1 and n2 is greater than 45°; the two vertices of the boundary edge are the boundary points.

[0046] S4.3 performs network reconstruction using the isotropic_remeshing function.

[0047] Specifically, S5 includes:

[0048] S5.1. Construct a KD-tree index for the edge point set. For each edge e in the grid, take the three-dimensional coordinates of its two vertices v1 and v2, query the Euclidean distances d1 and d2 from the two vertices to the nearest edge point, and set the distance d from the edge to the nearest edge point to min(d1, d2).

[0049] S5.2 For all boundary edges, if the distance d from the edge to the nearest edge point is less than the removal threshold, then the edge and its corresponding vertex are removed from the boundary edge set and the boundary point set; for each edge e in the mesh, if the Euclidean distances d1 and d2 from its two vertices to the nearest edge point are both less than the protection threshold, then the edge and its corresponding vertex are added to the boundary edge set and the boundary point set.

[0050] Specifically, S6 is as follows:

[0051] A progressive optimization strategy is adopted to update the KD tree index of the material surface point cloud and then reconstruct the mesh in multiple stages.

[0052] In the first phase, all grids are reconstructed using the first side length;

[0053] In the second stage, reconstruction is performed through multiple iterations with smaller side lengths. In each iteration with smaller side lengths:

[0054] (1) Divide the mesh surface into material surface area surface and non-material surface area surface. The criterion is: traverse all vertices of the surface. If the distance from any vertex to the nearest point of the material surface point cloud is less than D, then the mesh surface is divided into material surface area surface ... min If the distance is 1 meter, then the surface is classified as a material surface area; otherwise, it is classified as a non-material surface area area. Where, D... min The distance threshold used to distinguish between material surface areas and non-material surface areas;

[0055] (2) Differentiated reconstruction is performed on the two types of mesh surface sets respectively: the material surface area is reconstructed using the second side length, and the non-material surface area is reconstructed using the third side length;

[0056] The length of the first side is greater than the length of the second side, and the length of the second side is greater than the length of the third side.

[0057] Specifically, S7 involves: constructing a KD-tree index for the constraint vertex set using the constraint vertex set formed by the boundary edge vertices identified in S5, and identifying the D-values ​​near the boundary edge. ne Within the range of mesh vertices, the convex hull is calculated only for the local point set formed by these vertices. Vertices located inside the convex hull in the reconstructed mesh after S6 are moved along the convex hull normal direction to the convex hull surface, achieving precise filling of local depressions. Where D... ne It is the threshold of the approximation range when approximating the convex hull.

[0058] Preferably, the D min The first side length ranges from 0.5 to 0.6 meters; the second side length ranges from 0.05 to 0.1 meters; and the third side length ranges from 0.1 to 0.3 meters.

[0059] After adopting the above solution, the beneficial effects of the present invention are as follows:

[0060] (1) S2 performed boundary completion operation, providing accurate boundary data source for boundary edge identification and accurate boundary basis for subsequent point cloud volume calculation; the angle partitioning edge detection and interpolation completion algorithm proposed in this step addresses the technical defect of traditional convex hull algorithm in generating large errors when dealing with irregular material surface boundaries.

[0061] (2) The process of acquiring spliced ​​point clouds realizes the precise alignment, segmentation and splicing of the template point cloud and the material surface point cloud, avoiding edge misalignment caused by direct segmentation.

[0062] (3) When performing convex hull calculation and mesh reconstruction on the spliced ​​point cloud in S4, the boundary edges and boundary points of the spliced ​​point cloud were obtained. In S5, the boundary edges and boundary points were updated again. The inaccurate parts of the boundary edges and boundary points in S4 were removed using the accurate point set obtained in S2, and a relatively complete and accurate set of effective boundary edges was obtained.

[0063] (4) After network reconstruction in S4, a two-stage reconstruction and incremental optimization strategy was adopted again in S6, which greatly improved computational efficiency while ensuring the accuracy of key areas.

[0064] (5) The local convex hull approximation repair algorithm proposed in this invention within the influence range of the boundary edge can further solve the problem that the local depressions that may occur after mesh deformation affect the accuracy of volume calculation. Attached Figure Description

[0065] Figure 1 This is an overall flowchart of a specific implementation of the present invention;

[0066] Figure 2 This is a schematic diagram of the material surface point cloud and template point cloud in a specific embodiment of the present invention;

[0067] Figure 3 This refers to the boundary result obtained in S4 of a specific implementation of the present invention;

[0068] Figure 4 This is a schematic diagram showing the location of the final effective set of boundary edges in the point cloud in a specific implementation of the present invention;

[0069] Figure 5 This is a local display result of the closed mesh after directly using the general isotropic remeshing technique in a specific implementation of the present invention;

[0070] Figure 6 The image shows a partial view of the closed mesh after using the algorithm of this invention. The markings in the image are: 1 - the boundary where the material surface point cloud and the template point cloud are joined in the spliced ​​point cloud; 2 - edge-face boundary. Detailed Implementation

[0071] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0072] In the description of this invention, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Thus, a feature defined as "first" or "second" may explicitly or implicitly include one or more of the stated features. In the description of this invention, "a plurality of" means two or more, unless otherwise explicitly specified.

[0073] In the description of this invention, the term "for example" is used to mean "used as an example, illustration, or description." Any embodiment described as "for example" in this invention is not necessarily to be construed as being more preferred or advantageous than other embodiments. The following description is provided to enable any person skilled in the art to make and use the invention. Details are set forth in the following description for purposes of explanation. It should be understood that those skilled in the art will recognize that the invention can be made without using these specific details. In other instances, well-known structures and processes will not be described in detail to avoid obscuring the description of the invention with unnecessary detail. Therefore, the invention is not intended to be limited to the embodiments shown, but is consistent with the broadest scope of the principles and features disclosed herein.

[0074] A high-precision silo volume calculation method based on multi-scale mesh reconstruction and constraint protection, such as Figure 1 As shown, it includes the following steps:

[0075] S1. Preprocess the material surface point cloud data; perform data preprocessing and format standardization on the acquired material surface point cloud data, including skipping invalid data in the file header, uniformly processing tab and space delimiter formats, and detecting and filtering NaN outliers; data preprocessing can ensure that the input data meets the requirements of subsequent processing, and the specific preprocessing details are existing technologies.

[0076] S2. The point cloud of the pre-processed material surface is filled in by polar coordinate partitioning, and the edge points after filling in are used as the accurate boundary of the material surface.

[0077] S2.1. Project the pre-processed material surface point cloud (3D point cloud) onto the XY plane;

[0078] S2.2. In the XY plane, with the geometric center of the point cloud as the center, the 360-degree circle is evenly divided into K sectors. In practice, common values ​​for K are 36, 72, 90, 180, and 360. The higher the density of the material point cloud, the larger the value of K. At the same time, the larger the value of K, the finer the division, and the more accurate the subsequent edge completion. However, this will increase the amount of calculation. Therefore, in this specific implementation, the 360-degree circle is divided into 180 2-degree sectors, which is based on a comprehensive consideration of the amount of calculation and accuracy.

[0079] S2.3. Within each sector, traverse all points (these points are on a two-dimensional plane), calculate their polar coordinate angles, and find the edge point F that is farthest from the center. i The distance from its point to the center of the circle is r. i Where i represents the i-th sector; in each sector, there is usually only one edge point farthest from the center. If there are multiple edge points, one can be selected.

[0080] S2.4. For sectors without edge points, the nearest neighbor sector replacement completion strategy is adopted to complete their edge points. The completed edge points are added to the material surface point cloud and the sector is updated. Since the point density of the input material surface point cloud is uncontrollable, the point cloud density is different in different regions. There are some regions where the input material surface point cloud is relatively sparse. For example, in the 180 2-degree sectors in this specific implementation, due to the large number of sectors, it is common for the entire sector to have no points.

[0081] S2.5. Obtain the precise material surface boundary based on the edge points of all sectors;

[0082] For any sector without edge points, the nearest neighbor sector replacement completion strategy completes its edge points as follows:

[0083] Iterate through all sectors with edge points, calculate the angle difference between each sector and the current sector without edge points, and sort them in ascending order of angle difference to obtain the nearest neighbor sequence. The sector with the smallest angle difference from the current sector without edge points in the nearest neighbor sequence is p, and its edge point is F. p The distance from the edge point to the center of the circle is r. p , let r o Equal to the radius r of the current sector without edge points o The edge points are calculated, and the specific calculation process is as follows:

[0084] The x and y coordinates of the edge point of the current edgeless sector are as follows:

[0085] x=x c +r o cosθ

[0086] y=y c +r o sinθ

[0087] Where, x c , y c Let θ be the two-dimensional coordinates of the center of the circle, and θ be the angle of the current sector without edge points.

[0088] The z-coordinates of the edge points of the current edgeless sector are calculated using nearest neighbor inverse distance weighted interpolation.

[0089] z = (d q / (d p +d q ))·z p + (d p / (d p +d q ))·z q ;

[0090] Among them, z p z q These are the edge points F of sector p with the smallest angle difference and sector q with the second smallest angle difference in the nearest neighbor sequence. p F q z-coordinate, d p d q They are edge points F p F q Distance from the edge point of the current edgeless sector in the X–Y plane;

[0091] This completion provides an accurate boundary data source for boundary edge identification and a precise boundary basis for subsequent point cloud volume calculation. The angle partitioning edge detection and interpolation completion algorithm proposed in this step addresses the technical shortcomings of traditional convex hull algorithms, which produce large errors when dealing with irregular material surface boundaries.

[0092] S3. The template point cloud is segmented using the material surface point cloud after edge completion. The segmented template point cloud and the edge-completed material surface point cloud are then stitched together to obtain a stitched point cloud. The template point cloud is a standard template generated by acquiring the key model parameters of the silo, used to represent the silo point cloud in the empty state. The key model parameters of the silo include the radius of the cylindrical part, the height of the cylindrical part, and the height of the conical part. The template point cloud is combined with and segmented with the material surface point cloud in the material-containing state to generate a volume model of the raw material within the silo and to enable subsequent calculations.

[0093] The process of segmenting the template point cloud using the edge-completed material surface point cloud, and then stitching the segmented template point cloud with the edge-completed material surface point cloud to obtain a stitched point cloud includes:

[0094] S3.1. Construct a KD tree index for the material surface point cloud, perform a nearest neighbor search for each point in the template point cloud, find the nearest material surface point to each template point and calculate the height difference between the two on the z-axis;

[0095] It's important to understand that a KD-tree is a spatial indexing structure. It doesn't rely on the identity of the point cloud; instead, it searches based on spatial coordinates. Regardless of which point cloud a point comes from, as long as they exist in the same 3D space, the same KD-tree can be used for nearest neighbor search. For example, if point cloud A contains point a(1, 2, 3), A can use the KD-tree of point cloud B (essentially a structure describing the adjacency relationships of points in B) to find the nearest point on B to a. In simpler terms, B's KD-tree is a map of B; A uses this map to find the nearest point, regardless of whether the two points belong to the same point cloud.

[0096] S3.2. Traverse all template points. If the height of a template point on the z-axis is lower than its nearest material surface point, retain the template point. If the height of a template point on the z-axis is higher than its nearest material surface point, but the height difference does not exceed the height tolerance threshold ε, also retain the template point and adjust its z-axis height to the z-axis height of its nearest material surface point. In this specific implementation, ε is set to 0.02 meters. Template points other than the above two types are not retained.

[0097] S3.3. Combine the material surface point cloud and the retained template points into a single spliced ​​point cloud;

[0098] The above operations achieve precise alignment, segmentation, and stitching of the template point cloud and the material surface point cloud edges, avoiding edge misalignment caused by direct segmentation. For example... Figure 2As shown, the template point cloud is divided into two parts, gray and orange (upper and lower parts), according to the edge of the material surface point cloud. The gray and orange parts (upper and lower parts) together form the complete template point cloud, and the pink (the point cloud in the middle part) is the material surface point cloud. After S3, the orange and pink (the point cloud in the lower part and the middle part) are spliced ​​together to form a new point cloud. This is the spliced ​​point cloud, which is used to calculate the subsequent volume.

[0099] S4. Perform convex hull calculation and mesh reconstruction on the stitched point cloud, and obtain the boundary edges and boundary points of the stitched point cloud; S4 includes:

[0100] S4.1. The initial convex hull mesh for stitching point clouds is calculated using the CGAL computational geometry algorithm library convex_hull_3;

[0101] S4.2. Identify boundary edges and boundary points of the initial convex hull mesh; traverse each mesh edge e, take half edge h and the two corresponding sides f1 and f2 of its opposite half edge, calculate the normal vectors n1 and n2 of the two sides, and consider the edge as a boundary edge when the angle between the normal vectors n1 and n2 is greater than 45°; the two vertices of the boundary edge are the boundary points.

[0102] To clearly understand the solution of this invention, the related concepts such as mesh edge and half-edge are explained here. Mesh edge: A line segment connecting two vertices; Half-edge: Each mesh edge is decomposed into two directed half-edges, pointing in opposite directions; Opposite half-edge: The other half-edge corresponding to any half-edge; Function of half-edge: Half-edges record the orientation information of the face, and each half-edge belongs to a specific face.

[0103] S4.3 performs mesh reconstruction using the isotropic_remeshing function. In mesh reconstruction, the target edge length is set to 0.1 meters, the number of iterations is 1, and boundary protection is enabled with explicit passing of the boundary edge mapping (parameters in the isotropic_remeshing function) to ensure that the aforementioned boundary edges remain intact during reconstruction and relaxation, thus obtaining a regular and shape-preserving initial mesh structure, providing a reliable basis for subsequent mesh deformation optimization (S6). The boundary edge results obtained in S4 are as follows... Figure 3 As shown.

[0104] S5. Update the boundary edges and boundary points based on the distance from the edge points to the mesh edges;

[0105] S5.1. To address the issue of boundary features being easily lost during mesh optimization, a KD-tree index for the edge point set is constructed. For a mesh edge, the three-dimensional coordinates of its two vertices v1 and v2 are taken, and the Euclidean distances d1 and d2 from the two vertices to the nearest edge point are queried respectively. The distance d from the edge to the nearest edge point is set to min(d1, d2).

[0106] S5.2 For all boundary edges, if the distance d from the edge to the nearest edge point is less than the removal threshold, then the edge and its corresponding vertex are removed from the boundary edge set and the boundary point set; for each edge e in the mesh, if the Euclidean distances d1 and d2 from its two vertices to the nearest edge point are both less than the protection threshold, then the edge and its corresponding vertex are added to the boundary edge set and the boundary point set; in this specific implementation, the removal threshold is 0.5 meters and the protection threshold is 0.1 meters.

[0107] After S5.1 and S5.2, the final set of valid boundary edges consists of the following two types of edges:

[0108] (1) All edges on the grid whose vertices are within the protection threshold of the nearest edge point are the boundary 1 at the junction of the material surface point cloud and the template point cloud in the spliced ​​point cloud, such as Figure 4 As shown.

[0109] (2) The set of boundary edges remaining after removing the boundary in (1) is the boundary of the template point cloud model itself, which is edge-face boundary 2, and not the boundary at the junction of the material surface point cloud and the template point cloud, such as Figure 4 As shown.

[0110] The set of valid boundary points is the set of vertices corresponding to the set of valid boundary edges.

[0111] like Figure 4 As shown, the boundary edges obtained in S4 include the material surface point cloud boundary and the remaining geometric structure boundaries (edge-face intersections). In fact, the resulting boundary edge point set and boundary point set are in set form. Figure 4 The distinction has been made solely for ease of understanding; therefore, this invention will... Figure 4 A portion of the material surface point cloud outside the two boundaries needs to be removed. In S2, the completed material surface point cloud is obtained, which is a relatively accurate point set. This accurate point set is used to remove the material surface point cloud boundary portion from the boundary edges obtained in S4. After this removal, only the geometric structure boundary is obtained. Based on the geometric structure boundary and the material surface point cloud boundary corresponding to the accurate point set, a relatively complete and accurate set of effective boundary edges is obtained.

[0112] S6. Perform incremental optimization and reconstruction on the mesh; adopt an incremental optimization strategy to update the KD tree index of the material surface point cloud, and then perform multi-stage reconstruction on the mesh;

[0113] To overcome the technical bottleneck that fixed mesh density cannot adapt to complex material surface shapes, this invention achieves mesh density control based on material surface distance discrimination. Specifically, this invention employs a two-stage reconstruction method, which is still based on CGAL's isotropic remeshing technology, a common technique in this field.

[0114] In the first stage, all grids are reconstructed using the first side length, iterating 3 times; the first stage is actually the early stage with large side lengths; the value range of the first side length is 0.1-0.5 meters, and in this specific implementation, the first side length is 0.15 meters;

[0115] In the second stage, reconstruction is performed through multiple iterations of smaller side lengths; this second stage is actually the later stage of smaller side length reduction. In each iteration of smaller side length reduction:

[0116] (1) Divide the mesh surface into material surface area surface and non-material surface area surface. The discrimination criterion is: traverse all vertices of the surface. If the distance from any vertex to the nearest point of the material surface point cloud (the material surface point cloud obtained in S2) is less than D min If the value is in meters, then the surface is classified as a material surface area; otherwise, it is classified as a non-material surface area area. Where D... min The value range is 0.5-0.6 meters; in this specific implementation, D min The value is 0.5;

[0117] (2) Differentiated reconstruction is performed on the two types of mesh surfaces respectively: the material surface area is reconstructed using the second side length, and the non-material surface area is reconstructed using the third side length. Both are set with an iteration number of N21 and a relaxation step of 3. The second side length ranges from 0.05 to 0.1 meters, and the third side length ranges from 0.1 to 0.3 meters. In this specific implementation, the second side length is 0.08 meters and the third side length is 0.15 meters.

[0118] By employing a progressive optimization strategy of using a large side length in the early stage (0.15 meters × 3 times) and a small side length in the later stage (0.08 meters for the material surface + 0.15 meters × 1 time for the non-material surface), the computational efficiency is significantly improved while ensuring the accuracy of key areas, and the boundary edges are protected from damage throughout the process.

[0119] S7. Perform convex hull approximation repair on the boundary region of the reconstructed mesh to obtain an optimized closed mesh;

[0120] To address the issue of local depressions affecting the accuracy of volume calculations after mesh deformation, this invention proposes a local convex hull approximation repair algorithm within the influence range of boundary edges. A KD-tree index is constructed using the boundary vertex set identified in S5, identifying mesh vertices within a 0.5-meter radius of the boundary edge. The convex hull is calculated only for the local point set formed by these vertices. Vertices located inside the convex hull in the reconstructed mesh after S6 are moved along the convex hull normal direction to the convex hull surface, thus achieving accurate filling of local depressions.

[0121] S8. Perform volume calculations based on the optimized closed mesh; achieving a complete closed loop from raw point cloud to high-precision volume measurement. Please participate. Figure 5 andFigure 6 , Figure 5 This is a local display of the closed mesh after directly using the general isotropic remeshing technique. Figure 6 The image shows a local view of the closed mesh after using the algorithm of this invention. As can be seen, the edge accuracy is significantly improved after using the algorithm of this invention.

[0122] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention.

[0123] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A high-precision silo volume calculation method based on multi-scale grid reconstruction and constraint protection, characterized in that, Includes the following steps: S1. Preprocess the point cloud data of the material surface; S2. Perform edge completion on the pre-processed material surface point cloud, and use the completed edge points as the accurate material surface boundary; S3. The template point cloud is segmented by the material surface point cloud after edge completion, and the segmented template point cloud and the material surface point cloud after edge completion are spliced ​​together to obtain the spliced ​​point cloud. S4. Perform convex hull calculation and mesh reconstruction on the stitched point cloud, and obtain the boundary edges and boundary points of the stitched point cloud; S5. Update the boundary edges and boundary points based on the distance from the edge points to the mesh edges; S6. Perform progressive optimization and reconstruction of the mesh; S7. The optimized closed mesh is obtained by approximating and repairing the local convex hull within the influence range of the boundary edges and boundary points of the reconstructed mesh; S8. Perform volume calculation based on the optimized closed mesh; S2 includes: S2.

1. Project the pre-processed material surface point cloud onto the XY plane; S2.

2. In the XY plane, with the geometric center of the point cloud as the center, the 360-degree circle is evenly divided into K sectors; S2.

3. Within each sector, traverse all points, calculate their polar coordinate angles, and find the edge point F that is farthest from the center. i The distance from its point to the center of the circle is r. i Where i represents the i-th sector; S2.

4. For sectors without edge points, the nearest neighbor sector replacement completion strategy is adopted to complete their edge points, and the completed edge points are added to the material surface point cloud and the sector is updated. S2.

5. Obtain the precise material surface boundary based on the edge points of all sectors; In S2.4, for any sector without edge points, the nearest neighbor sector replacement completion strategy completes its edge points as follows: Iterate through all sectors with edge points, calculate the angle difference between each sector and the current sector without edge points, and sort them in ascending order of angle difference to obtain the nearest neighbor sequence. The sector with the smallest angle difference from the current sector without edge points in the nearest neighbor sequence is p, and its edge point is F. p The distance from the edge point to the center of the circle is r. p , let r o Equal to the radius r of the current sector without edge points o The edge points are calculated, and the specific calculation process is as follows: The x and y coordinates of the edge point of the current edgeless sector are as follows: x=x c +r o cosθ; y=y c +r o sinθ; Where, x c , y c Let θ be the two-dimensional coordinates of the center of the circle, and θ be the angle of the current sector without edge points. The z-coordinates of the edge points of the current edgeless sector are calculated using nearest neighbor inverse distance weighted interpolation. z = (d q / (d p +d q ))·z p + (d p / (d p +d q ))·z q ; Among them, z p z q These are the edge points F of sector p with the smallest angle difference and sector q with the second smallest angle difference in the nearest neighbor sequence. p F q z-coordinate, d p d q They are edge points F p F q Distance from the edge point of the current edgeless sector in the X–Y plane; S3 includes: S3.

1. Construct a KD tree index for the material surface point cloud, perform a nearest neighbor search for each point in the template point cloud, find the nearest material surface point to each template point and calculate the height difference between them; S3.

2. Traverse all template points. If the height of a template point is lower than its nearest material surface point, then retain the template point. If the height of a template point is higher than its nearest material surface point, but the height difference does not exceed the height tolerance threshold ε, then retain the template point as well, and adjust the height of the template point to the height of its nearest material surface point. Other template points besides the above two types are not retained. S3.

3. Combine the material surface point cloud and the retained template points into a single spliced ​​point cloud; S4 includes: S4.

1. Calculate the initial convex hull mesh for the stitched point cloud; S4.

2. Identify boundary edges and boundary points of the initial convex hull mesh; traverse each mesh edge e, take half edge h and the two corresponding sides f1 and f2 of its opposite half edge, calculate the normal vectors n1 and n2 of the two sides, and consider the edge as a boundary edge when the angle between the normal vectors n1 and n2 is greater than 45°; the two vertices of the boundary edge are the boundary points. S4.3 performs network reconstruction using the isotropic_remeshing function; S5 includes: S5.

1. Construct a KD-tree index for the edge point set. For each edge e in the grid, take the three-dimensional coordinates of its two vertices v1 and v2, query the Euclidean distances d1 and d2 from the two vertices to the nearest edge point, and set the distance d from the edge to the nearest edge point to min(d1, d2). S5.2 For all boundary edges, if the distance d from the edge to the nearest edge point is less than the removal threshold, then the edge and its corresponding vertex are removed from the boundary edge set and the boundary point set; for each edge e in the mesh, if the Euclidean distances d1 and d2 from its two vertices to the nearest edge point are both less than the protection threshold, then the edge and its corresponding vertex are added to the boundary edge set and the boundary point set. Specifically, S6 is: A progressive optimization strategy is adopted to update the KD tree index of the material surface point cloud and then reconstruct the mesh in multiple stages. In the first phase, all grids are reconstructed using the first side length; In the second stage, reconstruction is performed through multiple iterations with smaller side lengths. In each iteration with smaller side lengths: (1) Divide the mesh surface into material surface area surface and non-material surface area surface. The criterion is: traverse all vertices of the surface. If the distance from any vertex to the nearest point of the material surface point cloud is less than D, then the mesh surface is divided into material surface area surface ... min If the distance is less than 1 meter, then the surface is classified as a material surface area; otherwise, it is classified as a non-material surface area area. (2) Differentiated reconstruction is performed on the two types of mesh surface sets respectively: the material surface area is reconstructed using the second side length, and the non-material surface area is reconstructed using the third side length; The length of the first side is greater than the length of the second side, and the length of the second side is greater than the length of the third side.

2. The high-precision silo volume calculation method based on multi-scale grid reconstruction and constraint protection according to claim 1, characterized in that, Specifically, S7 involves constructing a KD-tree index for the constraint vertex set using the constraint vertex set formed by the boundary edge vertices identified in S5, and identifying the D-values ​​near the boundary edge. ne Within the range of mesh vertices, the convex hull is calculated only for the local point set formed by these vertices. Vertices located inside the convex hull in the reconstructed mesh after S6 are moved to the surface of the convex hull along the convex hull normal direction to achieve accurate filling of local depressions.

3. The high-precision silo volume calculation method based on multi-scale grid reconstruction and constraint protection according to claim 1, characterized in that, The D min The value range of the first side length is 0.5-0.6 meters; the value range of the second side length is 0.05-0.1 meters; and the value range of the third side length is 0.1-0.3 meters.

Images