A method and system for extracting ground points from massive point cloud based on voxel morphology

Abstract

The application discloses a kind of mass point cloud ground point extraction method and system based on voxel morphology, belong to surveying technical field, by adaptive block and kernel external storage strategy, the efficient block processing of large-scale point cloud is realized, and the continuity between blocks is guaranteed by introducing overlapping area mechanism.Voxelization is carried out according to point cloud resolution, and attributes such as voxel barycenter and point number are calculated, the ground surface value is obtained by vertical occlusion analysis and neighborhood propagation, and then the occlusion voxel, projection planarization are removed, the fitting degree of voxel and real ground surface is accurately described.Edge thinning and contour line generation are used to extract non-ground range, and non-ground point labeling is completed.Using voxel barycenter point as boundary constraint, a fitting surface is constructed by using thin plate spline method and sampling filling, and complete and continuous ground point cloud is generated, which provides a reliable foundation for high-precision digital ground model construction.The application has good adaptability, robustness and precision, and is suitable for mass point cloud ground extraction in various complex scenes.

Description

Technical Field

[0001] This invention belongs to the field of surveying and mapping technology, specifically relating to a method and system for extracting ground points from massive point clouds based on voxel morphology. Background Technology

[0002] LiDAR, as an active remote sensing technology, can directly acquire high-precision three-dimensional spatial coordinate information of target surfaces, possessing advantages such as strong penetration, high sampling density, and insensitivity to lighting conditions. With the increasing prevalence of airborne, vehicle-mounted, and ground-based laser scanning systems, LiDAR point cloud data has been widely applied in fields such as digital elevation model (DEM) generation, urban 3D modeling, forest resource surveys, geological disaster monitoring, and high-precision map construction for autonomous driving. In these applications, accurately separating ground points from non-ground points (i.e., ground filtering) from the raw point cloud is the most fundamental and crucial preprocessing step. Point cloud data contains points of features such as buildings, vegetation, and vehicles; only by accurately extracting ground points can a DEM that truly reflects the terrain undulations be constructed, ensuring the reliability of subsequent analysis results.

[0003] Currently, existing ground filtering methods mainly include those based on mathematical morphology, progressive triangulation, moving surface fitting, and machine learning. Among them, filtering algorithms based on mathematical morphology are widely used due to their high computational efficiency and ease of implementation. Traditional morphological filtering performs opening operations on point clouds by setting structural elements of different scales, thereby filtering out non-ground points. However, these methods generally rely on manually set window parameters and slope thresholds. In complex scenarios such as steep terrain, dense vegetation, or mixed ground features, their adaptability is poor, easily leading to excessive removal of ground points or incorrect retention of non-ground points, making it difficult to meet the filtering accuracy requirements of high-precision DEM production. In addition, with the rapid development of LiDAR technology, the scale of point cloud data has increased dramatically, with single scans often generating millions or even hundreds of millions of point clouds. The efficient organization and processing of massive unstructured point clouds has become a major bottleneck restricting the real-time performance and scalability of ground filtering. Voxelization methods map discrete point clouds into regular three-dimensional meshes, which can effectively compress data scale and regularize spatial structure, providing a unified computational basis for neighborhood analysis and spatial filtering. Existing research has attempted to combine voxelization techniques with morphological filtering to identify ground voxels through voxel-level morphological operations. However, existing voxel morphology methods still suffer from shortcomings such as difficulty in adapting voxel scale to terrain changes and a lack of systematic design for multi-morphological structural element combination strategies, resulting in poor robustness when moving between different scenes. Therefore, there is an urgent need to develop a point cloud ground point extraction method that can balance terrain adaptability, filtering accuracy, and computational efficiency.

[0004] To address the problem of efficiently extracting ground points from massive point clouds (i.e., ground filtering), the invention patent "A Ground Filtering Method for Laser Point Clouds Based on Convex and Concave Terrain Classification" [Application No.: 202410624897.9] extracts ground points by classifying the land surface into concave and convex terrains and proposing different ground point discrimination formulas for different terrains, avoiding a large number of terrain parameters, but it lacks adaptability to terrain. The invention patent "Ground Filtering Method, Device, Electronic Equipment and Medium for Steep Environments" [Application No.: 202411088392.1] uses virtual points to occlude non-ground points at the edge of the filtering area, which can improve the filtering performance of ground filtering algorithms in steep areas with large slopes, but it depends on the generation effect of virtual points. The invention patent "Adaptive Point Cloud Thinning Method Based on Path Point Neighborhood Domain and Ground Filtering" [Application No.: 202110918552.0] proposes an adaptive point cloud thinning method based on path point neighborhood domain and ground filtering, but it depends on importing path point data retained by the acquisition device.

[0005] In summary, existing laser point clouds are characterized by massive amounts of data, complex terrain, and complex environments, which leads to problems such as environmental unsuitability and unstable extraction accuracy in existing traditional point cloud ground point extraction methods. Summary of the Invention

[0006] To address the issues of massive data volume, complex terrain, and complex environment in existing laser point clouds, which lead to problems such as unsuitability for the environment and unstable extraction accuracy in traditional point cloud ground point extraction methods.

[0007] In a first aspect, this invention proposes a method for extracting ground points from massive point clouds based on voxel morphology, comprising:

[0008] The massive point cloud is divided into blocks to form several point cloud blocks, and point cloud preprocessing is performed on each point cloud block.

[0009] Based on the resolution of the original massive point cloud, each preprocessed point cloud block is voxelized.

[0010] The voxelized point cloud blocks are calculated by combining vertical occlusion analysis and neighborhood propagation to determine the voxel surface value, and the occluded non-bottom voxels are removed.

[0011] Based on the voxel surface values, the effective voxels are vertically projected onto the horizontal plane to generate a two-dimensional raster image, and the invalid value regions in the two-dimensional raster image are filled through neighborhood information propagation and bilateral filtering.

[0012] The two-dimensional raster image is optimized by attribute constraints to fill in missing data, and edge thinning is performed to extract effective surface boundaries. The effective surface boundary edge pixels are connected into a closed vector contour line by an eight-neighbor tracing algorithm and simplified to obtain the ground contour line.

[0013] Based on the ground contour line, traverse each point in all point cloud blocks, and use the ray method to determine whether the point is inside the contour line. Points inside the contour line are marked as non-ground points, and points outside the contour line are marked as ground points.

[0014] Based on the marked non-ground points, the voxel centroids around the holes are used as boundary constraints. A thin plate spline algorithm is used to fit a smooth surface and virtual ground points are generated by sampling at DEM resolution. Boundary fusion and filtering are performed on all point cloud blocks to obtain a complete and continuous ground point cloud.

[0015] Preferably, the step of dividing the massive point cloud into blocks to form several point cloud blocks includes:

[0016] Obtain basic information of massive point clouds, and perform preliminary virtual segmentation of massive point clouds based on the basic information, dividing the original spatial range into several virtual spaces, each virtual space corresponding to a point cloud block;

[0017] During the virtual segmentation process, the number of points in the current point cloud block is counted in real time. After the virtual segmentation is completed, it is ensured that the number of points in each point cloud block reaches the preset segmentation threshold. For point cloud blocks with fewer than 0.5 times the threshold, they are merged according to the spatial neighborhood topology. For point cloud blocks with more than 1.5 times the threshold, a second spatial segmentation is performed to generate smaller point cloud blocks.

[0018] The point cloud blocks, which are divided into several virtual spaces, are traversed once. Each point in the massive point cloud is allocated and filled into a unique virtual space and stored according to a preset format.

[0019] Preferably, the point cloud preprocessing for each point cloud block includes:

[0020] For each point cloud block, noise removal, isolated point removal, and data repair are performed, and adjacent points in overlapping areas are included during processing to ensure continuity between point cloud blocks.

[0021] Preferably, the step of voxelizing each preprocessed point cloud block according to the resolution of the original massive point cloud includes:

[0022] Based on the resolution of the original point cloud and the task scale, the voxel edge length algorithm is used to determine the voxel edge length, and voxel segmentation is performed based on the voxel edge length to form a number of voxels.

[0023] Calculate the centroid, number of centroids, and elevation value for each non-empty voxel, and clean up voxels with attribute values ​​of zero, while retaining the voxel spatial topological relationships.

[0024] Preferably, the process of calculating voxel surface values ​​for voxelized point cloud blocks through vertical occlusion analysis and neighborhood propagation, and removing occluded non-bottom voxels, includes:

[0025] For a number of three-dimensional voxel grids, traverse each voxel column from top to bottom along the vertical direction;

[0026] For voxels in the same column, the surface value of the lower voxel is assigned based on the distance and number of points of the upper voxel; if the upper voxel exists, it will have an occlusion effect on the lower voxel.

[0027] The surface values ​​are smoothed by propagating through the neighborhood, and finally the occluded non-bottom voxels are removed, while the potential surface voxels are preserved.

[0028] Preferably, the step of vertically projecting effective voxels onto a horizontal plane to generate a two-dimensional raster image based on the voxel surface values, and filling the invalid value regions in the two-dimensional raster image through neighborhood information propagation and bilateral filtering, includes:

[0029] All effective voxels are vertically projected onto a horizontal plane to generate a two-dimensional raster image; the value of each pixel is assigned the corresponding voxel's surface value;

[0030] During projection, if there are no valid voxels at a certain plane position, the corresponding pixel is marked as invalid.

[0031] Invalid value regions are extracted. When the range of invalid value regions is less than 1 / 2 of the minimum sensitive area of ​​ground features, low values ​​are filled by the boundary values ​​around the void range. The surface values ​​of voxels in the current voxel spatial neighborhood are taken, and the surface values ​​of the current voxels are updated using bilateral filtering to enhance the neighborhood edge information, filter and smooth flat areas of the same nature, and consider the coherence of the information of the boundary link areas of different point cloud blocks.

[0032] Preferably, the step of performing attribute constraint optimization to fill missing data in the two-dimensional raster image, and performing edge refinement to extract effective surface boundaries; connecting the effective surface boundary edge pixels into a closed ground contour line using an eight-neighbor tracing algorithm includes:

[0033] For the two-dimensional raster image, attribute constraints are performed based on synchronously acquired image data and after heterogeneous data processing and alignment.

[0034] For the portion of the land surface value in the attribute-constrained 2D raster image that is higher than a preset threshold, edge refinement is performed by using a weighted average of the neighborhood spatial location and the land surface value to extract the effective surface boundary.

[0035] Remove the unclosed areas in the effective table boundary, and then determine whether the vertical coordinate value within the closed effective table boundary is greater than the corresponding boundary value. If it is greater, retain it; otherwise, delete the effective table boundary.

[0036] The processed effective surface boundary edge pixels are connected into a closed vector contour line using an eight-neighbor tracing algorithm, and then simplified using Douglas-Pockmark polyline to remove redundant nodes, generating a smooth and compact ground contour line.

[0037] Preferably, the marker-based non-ground points are used with the voxel centroids around the holes as boundary constraints, and a smooth surface is fitted using the thin-plate spline method. Virtual ground points are generated by sampling at DEM resolution. Boundary fusion and filtering are performed on all point cloud blocks to obtain a complete and continuous ground point cloud, including:

[0038] Using the set of voxel centroid points on the cavity boundary as constraints, the thin plate spline algorithm is used to fit the marked non-ground points to a smooth terrain surface;

[0039] Inside the cavity, samples are taken at equal intervals according to the DEM resolution, and the elevation value of each sampling point is calculated to generate virtual ground points.

[0040] Based on the virtual ground points, sampling and filling are performed, and boundary fusion and filtering are applied to the filled area to ensure a continuous transition between the virtual points and the original ground points, thereby obtaining a complete and continuous ground point cloud.

[0041] Furthermore, the formula for calculating the elevation value is as follows:

[0042]

[0043] in, This is the elevation value. For affine coefficients, The number of boundary points. For index value, To sample points at equal intervals according to DEM resolution, For the coordinates of the control points, As the base weight, As basis functions, The distance between the two points is the Euclidean distance.

[0044] Secondly, this invention proposes a massive point cloud ground point extraction system based on voxel morphology, comprising:

[0045] The block segmentation module is used to segment massive point clouds into blocks, forming several point cloud blocks, and to perform point cloud preprocessing on each point cloud block.

[0046] The voxelization module is used to voxelize each preprocessed point cloud block according to the resolution of the original massive point cloud.

[0047] The surface value calculation module is used to jointly calculate the voxel surface value of point cloud blocks after voxelization through vertical occlusion analysis and neighborhood propagation, and remove occluded non-bottom voxels.

[0048] The projection planarization module is used to vertically project the effective voxels onto the horizontal plane based on the voxel surface values ​​to generate a two-dimensional raster image, and fill the invalid value regions in the two-dimensional raster image through neighborhood information propagation and bilateral filtering.

[0049] The ground contour extraction module is used to perform attribute constraint optimization and fill in missing data on the two-dimensional raster image, and to perform edge refinement operations to extract effective surface boundaries; the effective surface boundary edge pixels are connected into closed vector contour lines through an eight-neighbor tracing algorithm and simplified to obtain the ground contour lines;

[0050] The non-ground point marking module is used to traverse each point in all point cloud blocks according to the ground contour line, and use the ray method to determine whether the point is inside the contour line. The internal points are marked as non-ground points, and the external points are marked as ground points.

[0051] The virtual hole filling module is used to fill non-ground points based on markers. It uses the centroid points of voxels around the holes as boundary constraints, fits smooth surfaces using a thin plate spline algorithm, and generates virtual ground points by sampling at DEM resolution. It then performs boundary fusion and filtering on all point cloud blocks to obtain a complete and continuous ground point cloud.

[0052] In another aspect, the present invention also provides a computing device, comprising: one or more processors;

[0053] A processor is used to execute one or more programs;

[0054] When the one or more programs are executed by the one or more processors, a method for extracting ground points from massive point clouds based on voxel morphology, as described above, is implemented.

[0055] In another aspect, the present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed, implements the above-described method for extracting ground points from massive point clouds based on voxel morphology.

[0056] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0057] This invention provides a method for extracting ground points from massive point clouds based on voxel morphology. First, the massive data is adaptively segmented and stored to reduce dependence on hardware computing performance and storage capacity, making the method universally applicable. Then, the segmented point cloud blocks are processed, including noise removal and voxelization. Based on the voxelized data, surface values ​​are calculated to guide subsequent projection planarization, generating a planar map that facilitates non-ground contour extraction. Furthermore, to ensure the integrity of remaining ground points, virtual hole filling is performed. Finally, the processing results of multiple point cloud blocks are merged and integrated. This approach avoids dependence on hardware devices (such as large memory and expensive graphics cards), efficiently achieving ground filtering from massive point clouds without relying on complex template control parameters and terrain limitations.

[0058] Furthermore, this invention employs an adaptive block partitioning strategy, performing initial and secondary spatial segmentation based on the average density of the point cloud and a preset block partitioning threshold. Combined with an external data storage mechanism, this achieves balanced block partitioning and efficient read / write operations for massive point clouds under typical hardware conditions. Simultaneously, voxelization converts irregular point clouds into regular voxel meshes, and centroid attribute settings further compress the data volume. These mechanisms effectively avoid memory overflow caused by loading large-scale point clouds all at once, significantly reducing the dependence of massive point cloud processing on hardware performance, improving computational efficiency and scalability, significantly lowering the hardware threshold, and enhancing the scalability and processing efficiency of the method in large-scale mapping tasks.

[0059] Furthermore, this invention defines voxel surface values, recursively calculates the occlusion effect from top to bottom along the vertical direction, and removes non-surface voxels containing lower voxels from bottom to top, retaining only the bottommost voxels. This method comprehensively considers vertical occlusion relationships and neighborhood information propagation, enabling accurate differentiation between ground and non-ground objects in complex scenarios such as continuous / discontinuous boundaries, faults, and laser penetration. This invention introduces voxel surface values ​​to comprehensively quantify the degree of fit between voxels and the real ground surface. Vertical occlusion analysis recursively calculates the occlusion effect along the vertical direction, enabling accurate differentiation between ground and non-ground objects (such as vegetation, buildings, and bridges), covering typical complex scenarios such as continuous / discontinuous boundaries, faults, and laser penetration. Subsequently, by removing occluded voxels and retaining only the voxel layer closest to the ground surface, after projection planarization, it combines neighborhood information propagation and bilateral filtering to fill gaps and smooth surface values, effectively suppressing noise and local abrupt changes. This method has better ground point extraction accuracy than traditional filtering algorithms in difficult areas such as steep slopes, under bridges, and vegetated areas.

[0060] Furthermore, this invention employs a virtual hole-filling strategy after marking non-ground points. Using voxel centroids as boundary constraints, a smooth terrain surface is constructed using the thin-plate spline method. Virtual ground points are generated by sampling at equal intervals according to DEM resolution, and finally, boundary fusion and filtering are performed. This mechanism effectively repairs holes caused by occlusion, insufficient reflection, or missing point clouds, resulting in a continuous and complete final ground point cloud. It avoids the problem in traditional methods where hole regions directly lead to DEM distortion or unusable data, providing a reliable data foundation for high-precision generation of digital terrain models. Attached Figure Description

[0061] Figure 1 This is a flowchart of a point cloud ground point extraction method based on voxel morphology in Embodiment 1 of the present invention;

[0062] Figure 2 This refers to the spatial neighborhood of the laser point cloud voxel in Embodiment 1 of the present invention;

[0063] Figure 3 This illustrates the spatial distribution of laser point clouds at ground and non-ground points in Embodiment 1 of the present invention.

[0064] Figure 4 This is a schematic diagram of the surface boundary that can be extracted after plane projection of voxel evaluation values ​​in Embodiment 1 of the present invention;

[0065] Figure 5 This is a schematic diagram of the refinement of the actual surface boundary in Embodiment 1 of the present invention;

[0066] Figure 6 This is a schematic diagram of the electronic device structure in Embodiment 3 of the present invention. Detailed Implementation

[0067] Example 1:

[0068] A method for extracting ground points from massive point clouds based on voxel morphology, such as Figure 1 As shown, it includes:

[0069] Step 1: Divide the massive point cloud into blocks to form several point cloud blocks;

[0070] Step 2: Perform point cloud preprocessing on each point cloud block;

[0071] Step 3: Based on the resolution of the original massive point cloud, voxelize each preprocessed point cloud block to form several three-dimensional voxel meshes.

[0072] Step 4: Calculate the voxel surface value for several 3D voxel meshes by combining vertical occlusion analysis and neighborhood propagation, and remove occluded non-bottom voxels.

[0073] Step 5: Vertically project the effective voxels onto the horizontal plane to generate a two-dimensional raster image, and mark holes as invalid values; fill small holes through neighborhood information propagation and bilateral filtering to smooth abrupt changes in surface values;

[0074] Step 6: Perform attribute constraint optimization on the two-dimensional raster image to fill in missing data, and perform edge refinement operation to extract effective surface boundaries; connect the effective surface boundary edge pixels into a closed vector contour line through an eight-neighbor tracing algorithm and simplify it to obtain the ground contour line;

[0075] Step 7: Based on the ground contour line, traverse each point in the point cloud block and use the ray method to determine whether the point is inside the contour line. Points inside the contour line are marked as non-ground points, and points outside the contour line are marked as ground points.

[0076] Step 8: Based on the marked non-ground points, using the centroid points of the voxels around the holes as boundary constraints, fit a smooth surface using the thin plate spline method, and generate virtual ground points by sampling at DEM resolution; perform boundary fusion and filtering to obtain a complete and continuous ground point cloud.

[0077] Specifically, in this embodiment of the invention, step 1 involves adaptive block partitioning. Directly processing massive point clouds places high demands on hardware performance. Therefore, to improve the adaptability of massive point clouds in general hardware environments, it is necessary to partition them into blocks. To quickly achieve balanced block partitioning of massive point clouds, this invention uses a spatial adaptive block partitioning strategy.

[0078] As should be further explained, the adaptive block segmentation process in this embodiment of the invention includes:

[0079] Step 1.1, Initial Spatial Segmentation. Obtain basic information about the massive point cloud (spatial extent, number of points, file header information, etc.). Based on this information, perform preliminary virtual segmentation of the massive point cloud, dividing the original spatial extent into several virtual spaces, each corresponding to a point cloud block. During the preliminary virtual segmentation process, the size of the segmented space is calculated based on the average density of the massive point cloud, determining the length of the virtual space corresponding to each point cloud block reaching a preset segmentation threshold (e.g., 1 million blocks).

[0080] Step 1.2, Spatial Adaptive Segmentation. During the virtual segmentation process, the number of points in the current point cloud block is counted in real time. After the virtual segmentation ends, ensure that the number of points in each point cloud block reaches the preset segmentation threshold (e.g., 5 million). For point cloud blocks with fewer than 0.5 times the threshold, they are merged according to the spatial neighborhood topology; for point cloud blocks with more than 1.5 times the threshold, a second spatial segmentation is performed to generate smaller point cloud blocks.

[0081] Step 1.3, External Data Storage. Virtual segmentation mainly involves dividing the original massive point cloud space into several virtual spaces (corresponding to various point cloud block files). Then, the massive point cloud is traversed once, and each point in the massive point cloud is allocated and filled into a unique virtual space. To avoid excessive memory consumption due to the large number of points in the virtual space, an I / O (input / output) caching mechanism is set up for each virtual space. When the number of points stored in the virtual space exceeds the cache threshold (e.g., 10,000 points), the data in the cache is written to the corresponding file. Each file write follows the file header information format of the original massive point cloud to avoid information loss and maintain data alignment.

[0082] This invention proposes an adaptive block partitioning strategy that calculates block length based on the average density of point clouds and dynamically adjusts the virtual space segmentation. Through initial segmentation and secondary segmentation, it ensures that the number of points in each point cloud block is close to a preset threshold; at the same time, it adopts an I / O caching mechanism to write in batches according to the file header format during the traversal and filling process, thereby achieving efficient off-core storage and balanced block partitioning of massive point clouds in a general hardware environment.

[0083] This invention employs an adaptive block partitioning strategy, performing initial and secondary spatial segmentation based on the average density of the point cloud and a preset block partitioning threshold. Combined with an external data storage mechanism, it achieves balanced block partitioning and efficient read / write operations for massive point clouds under typical hardware environments. Simultaneously, voxelization transforms irregular point clouds into regular voxel meshes, and centroid attribute settings further compress the data volume. These mechanisms effectively avoid memory overflow caused by loading large-scale point clouds all at once, significantly reducing the dependence of massive point cloud processing on hardware performance, improving computational efficiency and scalability, significantly lowering the hardware threshold, and enhancing the scalability and processing efficiency of the method in large-scale mapping tasks.

[0084] Specifically, in step 2 of this embodiment of the invention, point cloud preprocessing is performed on each point cloud block, including:

[0085] For each point cloud block, noise removal, isolated point removal, and data repair are performed, and adjacent points in overlapping areas are included during processing to ensure continuity between blocks.

[0086] Further explanation is needed regarding the point cloud preprocessing in step 2 of this embodiment. Point cloud preprocessing operations are performed on each point cloud block, mainly including noise removal, isolated point removal, and data repair, ensuring that the point cloud data in each block is clean, continuous, and complete, providing a reliable data foundation for subsequent ground point extraction. Since the point cloud preprocessing process is performed within point cloud blocks, there may be information loss at the boundaries of the point cloud blocks compared to the original massive point cloud. Therefore, the processing of each point cloud block needs to consider the information of adjacent point cloud blocks, designing overlapping areas for each point cloud block. Specifically, this includes:

[0087] Step 2.1 Noise Removal. Noise removal is used to remove discrete noise points caused by sensor errors or environmental interference. For each point cloud block, firstly, the average distance between each point and its k (e.g., 10) neighboring points in its neighborhood is calculated. Assuming that this distance follows a Gaussian distribution, points whose average distance exceeds the mean μ plus λ times (e.g., 3 times) the standard deviation σ are identified as noise and removed. Secondly, a radius filtering algorithm is used to count the number of neighboring points within a given radius of each point. If the number is less than a preset threshold, it is considered an outlier and removed. Since the information at the boundary of the point cloud block may be discontinuous with adjacent blocks, the data of adjacent points in the overlapping area needs to be included in the processing. That is, the extended area of ​​each point cloud block (the part overlapping with adjacent blocks) is also included in the statistical calculation. After removal, only the result of the point cloud within the block is retained, thereby avoiding misjudgment due to missing information at the boundary.

[0088] Step 2.2 Isolated Point Removal. Isolated points refer to small, isolated clusters that are spatially separated from the main point cloud, typically caused by low-reflectivity objects or multipath reflection. Removal is performed using a region-growing clustering method based on Euclidean distance. A distance threshold and a minimum cluster size threshold are set, and points within a point cloud block are clustered. Clusters smaller than the minimum threshold are identified as isolated points and deleted. To prevent incorrectly identifying clusters near the boundaries as isolated points due to point cloud block segmentation, overlapping region information is used to temporarily merge the overlapping areas of the current point cloud block with its neighboring blocks before performing cluster analysis. This ensures that continuous structures located at the block edge but belonging to the main point cloud are not incorrectly removed. After isolated point removal, the original block boundaries are restored to maintain the independence of each point cloud block's data.

[0089] Step 2.3 Data Repair. Data repair is used to fill in voids in the point cloud caused by occlusion or insufficient reflection, such as under bridges or in water areas. Before repair, the boundaries of voids within the current point cloud block are first detected, and then the cause of the voids is determined using additional data (field surveys, image data, etc.). For voids caused by occlusion, interpolation filling is performed by extracting boundary points of the ground (occluded object) to construct closure conditions. For voids caused by insufficient reflection, interpolation filling is performed directly using boundary data. During the repair process, it is necessary to ensure that the density of newly added points is consistent with the average density of the original point cloud and that no additional noise is introduced. Simultaneously, the repaired point cloud should maintain geometric continuity with adjacent point clouds in overlapping areas to avoid seams. After data repair is completed, a smoothing filter is applied to the repaired area to ensure a natural surface transition.

[0090] This invention defines voxel surface values, recursively calculates occlusion effects from top to bottom along the vertical direction, and removes non-surface voxels containing lower voxels from bottom to top, retaining only the bottommost voxels. This method comprehensively considers vertical occlusion relationships and neighborhood information propagation, and can accurately distinguish between ground and non-ground objects in complex scenes such as continuous / discontinuous boundaries, faults, and laser penetration.

[0091] This invention introduces a voxel-based surface value comprehensive quantification method to determine the degree of fit between voxels and the real ground surface. Vertical occlusion analysis recursively calculates the occlusion effect along the vertical direction, accurately distinguishing between ground and non-ground objects (such as vegetation, buildings, and bridges), and covering typical complex scenes such as continuous and discontinuous boundaries, faults, and laser penetration. Subsequently, by removing occluded voxels and retaining only the voxel layer closest to the ground surface, the projection is planarized, and then combined with neighborhood information propagation and bilateral filtering to fill gaps and smooth the surface values, effectively suppressing noise and local abrupt changes. This method outperforms traditional filtering algorithms in ground point extraction accuracy in difficult areas such as steep slopes, under bridges, and vegetated areas. Through voxel morphology and occlusion analysis, the accuracy and robustness of ground point extraction in complex scenes are improved.

[0092] Specifically, in step 3 of this embodiment, voxelization is performed. Based on the resolution of the original massive point cloud, each preprocessed point cloud block is voxelized. Voxelization aims to convert irregularly distributed point cloud data into a regularly arranged three-dimensional voxel mesh, such as... Figure 2 As shown, while effectively reducing the amount of data, it suppresses residual local disturbances and minor fluctuations in the point cloud, providing a stable data structure for subsequent ground point extraction.

[0093] Step 3.1 Voxel Segmentation. The determination of voxel size depends on the resolution of the original point cloud and the scale of the task. The resolution of the original point cloud is usually represented by the point spacing of the acquisition device. In this invention, the voxel side length vl is calculated using formula (1), where min(*) is the minimum value calculation function. If the voxel is too large, it will be over-smoothed, resulting in the loss of terrain details; if the voxel is too small, the downsampling effect will be insignificant and the noise suppression capability will be limited.

[0094] vl = min(2×resolution, 10×scale) (Formula 1)

[0095] Step 3.2 Voxel Attribute Settings. Calculate the centroid (voxel_cg) for each non-empty voxel. This is achieved by using the average coordinates of all points within the voxel as the representative point of that voxel. Compared to the geometric center point, the centroid (voxel_cg) more accurately reflects the spatial distribution characteristics of the point cloud. Especially when the point cloud distribution is uneven, the centroid can effectively reduce positional deviations and has good ground surface representation capabilities. It can be used for subsequent virtual hole filling and rapid generation of the corresponding DEM (Digital Terrain Model) based on ground points. Record the number of original points (voxel_cn) within each voxel. This attribute is used for subsequent evaluation of the surface values ​​within the voxel. Record the vertical coordinate value (voxel_z) (elevation value) of each voxel for subsequent identification of uneven ground surfaces. The combined setting of the centroid and the number of points provides accurate and reliable voxel feature input for voxel morphology-based ground point extraction.

[0096] Step 3.3 Voxel Cleaning. Clean up voxels with a value of zero (0 points), while retaining and recording the voxel spatial topology. The voxel spatial topology records not only spatial adjacency relationships, but also all vertical relationships in the same vertical direction. These vertical relationships are not affected by empty voxels and are mainly used to record vertical occlusion relationships.

[0097] Specifically, in step 4 of this embodiment, the surface value is calculated. The voxel surface value `voxel_v` is used to quantitatively characterize the degree of fit between each voxel and the real surface. The calculation of this value comprehensively considers two factors: First, the voxel occlusion relationship in the vertical direction. Analyzing the voxel sequence from top to bottom in the vertical direction, if the upper layer voxels exist and the lower layer voxels are occluded, the calculation of the surface value of the lower layer voxels is related to the positions of all upper layer voxels (vertical distance from the assigned voxel) and the number of points `voxel_cn`, thus simulating the supporting and occluding effect of the surface on the objects above, and preserving the topological relationship in the vertical direction. Second, the information propagation of neighboring voxels. The surface has spatial continuity, and the surface value of a single voxel needs to be weighted and smoothed through the surface values ​​of its horizontal and diagonal neighboring voxels to avoid abrupt changes in surface value caused by local noise or voids. By combining vertical occlusion analysis with neighborhood propagation, the obtained surface values ​​can accurately reflect the probability that voxels belong to the real surface, providing a reliable basis for subsequent surface morphology calculations.

[0098] Step 4.1 Vertical Occlusion Analysis. Traverse each voxel column from top to bottom along the vertical direction (Z-axis). For voxels in the same column (vertical column), if the upper voxel exists (i.e., is not empty), it will have an occlusion effect on the lower voxel.

[0099] The specific rules for assigning surface values ​​are as follows: when a lower voxel is occluded by an upper voxel, the current voxel's surface value... With all voxel positions above it (vertically). With weight Related to, is the sum of the products of all voxel positions and their weights, voxel positions The distance from the current voxel in the vertical direction. The positively correlated variable with the number of voxel points (voxel_cn) is shown in Formula 2. If there is no occlusion, the surface value of the lower voxels remains at its initial value (e.g., 0). This recursive relationship can simulate the spatial obstruction effect of upper-layer objects (such as vegetation and buildings) on lower-layer surface points in a real terrain, causing voxels actually located on the surface to receive lower surface values, while occluded non-surface voxels receive higher surface values.

[0100] (Formula 2)

[0101] like Figure 3 As shown, the point cloud distribution at the boundary between the ground and non-ground (land features, obstructions) encompasses two main categories: continuous and discontinuous, basically covering all typical scenarios in surveying and production. In the continuous case, Figure 3 Figures a and b illustrate two scenarios, "crossing" and "cutting walls," caused by voxel segmentation, which typically occur when acquiring point clouds using tilted laser scanning. Figure 3 The image in section 'c' shows the fault between the ground point and the obstruction, which is commonly found under the bridge deck. Figure 3 The figure in Figure d illustrates a scenario where insufficient density of the shading material allows the laser point to penetrate and reach the ground, which often occurs in vegetated areas. In either case, the surface value of the voxel can characterize the shading boundary through its large numerical properties.

[0102] Step 4.2 Remove Occlusion Voxels. In the current point cloud block, based on the voxelization results, traverse each column of voxels vertically from bottom to top. For each voxel, if there are no other valid voxels in its vertically downward direction (negative Z-axis direction), it indicates that the voxel is located at the bottom of the column and belongs to the unoccluded potential surface voxel, and is retained; conversely, if there are other voxels below, it indicates that the current voxel is occluded by the voxels below and belongs to the non-surface voxel (such as vegetation, buildings, etc.), and is deleted. This strategy effectively removes all non-bottom-level occluding voxels, retaining only the voxel layer closest to the real surface, significantly reducing interference in subsequent processing.

[0103] Specifically, in step 5 of this embodiment, the projection is planarized. After calculating the surface value, each planar position has only one effective voxel in the vertical direction, and there is no overlapping information between voxels. Then, all effective voxels are projected onto the plane to generate a two-dimensional raster image in planar space. Its resolution is consistent with the voxel side length, and the planar pixel value is the surface value of the voxel. Then, neighborhood information propagation and smoothing optimization are performed.

[0104] Step 5.1: Vertically project all these valid voxels onto a horizontal plane to generate a 2D raster image. Each pixel's value is assigned the corresponding voxel's surface value, reflecting how closely the voxel fits the real surface at that location (lower values ​​indicate closer proximity to the surface). During projection, if a plane location has no valid voxels (e.g., a data hole area), the corresponding pixel is marked as invalid (e.g., -1), and will be interpolated or ignored in subsequent processing. Projection planarization reduces the 3D voxel data to a 2D surface feature map, providing an intuitive and efficient data structure for subsequent surface morphology operations and ground point cloud extraction.

[0105] Step 5.2, Neighborhood Information Propagation and Smoothing. Leveraging the spatial continuity of the land surface, the surface value of each voxel is propagated using a neighborhood-weighted method. First, the hole range (invalid value region) is extracted. When the hole range is less than half of the minimum sensitive area for land features, low-value filling is performed using the boundary values ​​around the hole range. The surface value of the voxel within its spatial neighborhood is then taken, and a bilateral filter is used to update the surface value of the current voxel. This enhances neighborhood edge information, smooths flat areas of similar nature, and considers the coherence of information in the boundary links of different point cloud blocks. This step effectively fills in abrupt changes in surface values ​​caused by local noise or data holes, ensuring a smooth spatial transition of surface values ​​and enhancing robustness to the real land surface morphology.

[0106] Specifically, step 6 of this embodiment extracts non-ground contours. After projection planarization, the generated two-dimensional raster image shows areas corresponding to ground structures (such as buildings, bridges, and elevated facilities) and vegetation with high surface values, forming a significant contrast with surrounding ground areas with low surface values, thus presenting clear contour information. First, attribute constraint optimization is performed on the image to fill in missing data. Then, edge refinement is performed to extract effective surface boundaries. Next, contour line generation is performed, connecting the refined effective surface boundary edge pixels into a closed vector contour line using an eight-neighbor tracing algorithm.

[0107] Step 6.1 Attribute Constraints. Typically, during laser point cloud acquisition, image data is acquired simultaneously. After heterogeneous data processing and alignment, information such as structure boundaries and vegetation boundaries in the image data is used to further compensate for and optimize non-ground contours in the 2D raster image. This prevents discontinuities and lack of prominence in non-ground contours due to occlusion.

[0108] Step 6.2 Edge refinement. For example... Figure 4 As shown, after projection and planarization, the true surface boundary is hidden in areas with higher surface values, due to the existence of (such as...) Figure 3 -a) In cases where the earth crosses a wall, the extreme value at the ground surface is not the ground boundary, and there is a translational phenomenon (e.g., ...). Figure 5As shown in the figure, further refinement is needed in areas with high surface values, using a weighted average of neighborhood spatial location and surface value. However, in order to further eliminate ground errors caused by spatial agitation of laser scanning point cloud, the extracted effective surface boundary should be expanded outward based on the surface boundary to find the first abrupt change point in surface value, or directly expanded to half voxel length.

[0109] Step 6.3 Contour Generation. After obtaining the effective surface boundaries, firstly, unclosed regions (such as steep slopes) are removed. Then, it is determined whether the vertical coordinate value voxel_z within the closed effective surface boundary is greater than the corresponding boundary value. If it is greater, it is retained (convex terrain regions, which include structures, vegetation, etc.); otherwise, the effective surface boundary is deleted (concave terrain regions). The processed effective surface boundary edge pixels are connected into a closed vector contour line using an eight-neighbor tracing algorithm, and then simplified using the Douglas-Peucker algorithm to remove redundant nodes, generating a smooth and compact ground contour line. Finally, all non-ground regions are extracted to provide accurate spatial constraints for subsequent removal or classification of non-ground point clouds.

[0110] Specifically, in step 7 of this embodiment, non-ground points are marked. Based on the extracted non-ground contour lines, each point within the current point cloud block is traversed sequentially. The spatial relationship between the points and the contour lines is determined using the ray method or the wrap-around number method. Points located inside the contour lines are marked as non-ground points, and points located outside the contour lines are marked as ground points.

[0111] The labeled point cloud can be directly used for subsequent ground model generation or non-ground object removal. The labeling results are stored in the point cloud file as attribute fields. During the labeling of non-ground points, the spatial resolution remains consistent with the original point cloud. The voxel edge length is only used to improve efficiency while enhancing the robustness of ground extraction accuracy.

[0112] Specifically, step 8 of this embodiment involves virtual hole filling. After marking non-ground points, a large number of holes remain among the remaining ground points, which is detrimental to the subsequent generation of a DEM (Digital Terrain Model). Therefore, smooth filling is required based on the hole boundary value information (boundary non-ground points). The hole boundary information is a set of boundary points composed of voxel centroids (voxel_cg). These centroids are located in the ground area surrounding the holes and can accurately reflect the local terrain trend. After constructing a fitted surface and sampling filling, a complete and continuous ground point cloud can be obtained, providing a reliable data foundation for high-precision DEM generation.

[0113] Step 8.1 Constructing the Fitted Surface

[0114] Using the set of voxel centroids at non-ground points on the cavity boundary as constraints (control conditions), a smooth terrain surface is fitted using the thin-plate spline method, as shown in Equation 3. Then, samples are taken at equal intervals according to DEM resolution inside the cavity, and the elevation values ​​of each sampling point are calculated to generate virtual ground points.

[0115] (Formula 3)

[0116] in the formula To sample points at equal intervals according to DEM resolution, This is the elevation value. For affine coefficients, The number of boundary points. For index value, These are the coordinates of the control points (coordinates of the centroid of the voxel in the plane). As the base weight, As basis functions, Let be the Euclidean distance between the two points. Basis functions are used to gently propagate the influence of the control points inversely proportional to the distance.

[0117] Step 8.2 Sampling and Filling

[0118] During the surface fitting process, n equations are first constructed using the centroids of n boundary voxels (voxel_cg) to solve for the affine coefficients and basis weights. Then, the corresponding elevation values ​​are calculated by sequentially traversing the sampling points inside the holes.

[0119] During the filling process, the smoothness of the interpolated surface needs to be controlled to avoid abnormal bumps or depressions. Finally, boundary blending and filtering are performed on the filled area to ensure a continuous transition between virtual points and original ground points without any seams. Several processed point cloud blocks are then merged to complete the extraction of massive point cloud ground points.

[0120] This invention proposes a void virtual filling technology that integrates centroid constraints and thin plate splines. The method uses voxel centroids as void boundary constraints, constructs terrain surfaces using thin plate splines, generates virtual ground points by sampling at equal intervals, and then performs boundary fusion and filtering. This effectively repairs data voids caused by occlusion or insufficient reflection, ensures the continuity and integrity of ground point clouds, and supports the reliable generation of high-precision DEMs.

[0121] This invention, after marking non-ground points, further employs a virtual hole-filling strategy to ensure the integrity of ground point cloud data and support high-precision DEM generation. Using voxel centroids as boundary constraints, a smooth terrain surface is constructed using the thin-plate spline method, and virtual ground points are generated by sampling at equal intervals according to DEM resolution. Finally, boundary fusion and filtering are performed. This mechanism effectively repairs holes caused by occlusion, insufficient reflection, or missing point clouds, resulting in a continuous and complete final ground point cloud. It avoids the problem in traditional methods where hole regions directly lead to DEM distortion or unusable data, providing a reliable data foundation for high-precision generation of digital terrain models.

[0122] This invention presents a method for extracting ground points from massive point clouds based on voxel morphology. Addressing the challenges of high hardware dependence, high data noise, and complex occlusion in massive point cloud datasets, the method first achieves efficient block processing of large-scale point clouds through adaptive block segmentation and off-core storage strategies, and introduces an overlapping region mechanism to ensure continuity between blocks. Secondly, noise removal, island removal, and data repair are performed on each point cloud block to improve data quality. Then, voxelization is performed based on the point cloud resolution, and attributes such as voxel centroid and number of points are calculated. Surface values ​​are obtained through vertical occlusion analysis and neighborhood propagation, thereby removing occluded voxels and planarizing the projection, accurately characterizing the degree of fit between voxels and the real surface. Based on this, non-ground areas are extracted through edge refinement and contour generation, completing the marking of non-ground points. Finally, for ground point holes, voxel centroids are used as boundary constraints, and a fitted surface is constructed using the thin-plate spline method and sampled for filling, generating a complete and continuous ground point cloud, providing a reliable foundation for the construction of a high-precision digital terrain model (DEM). This invention combines good hardware adaptability, robustness, and accuracy, making it suitable for extracting ground data from massive point clouds in various complex scenarios. Specifically, it includes the following steps:

[0123] Step 1: Obtain basic point cloud information and perform initial virtual segmentation based on average density; merge neighborhoods for blocks with fewer than the threshold (0.5 times), and perform secondary segmentation for blocks with more than the threshold (1.5 times); use an I / O caching mechanism to allocate points to virtual space and write them to files to maintain data alignment.

[0124] Step 2: Perform noise removal (statistical filtering + radius filtering), island point removal (Euclidean distance clustering), and data repair (interpolation filling based on the cause of holes) on each point cloud block. During processing, adjacent points in overlapping areas are included to ensure continuity between blocks.

[0125] Step 3: Determine the voxel side length based on the point cloud resolution and scale, calculate the centroid, number of points, and elevation value of each non-empty voxel, clean up empty voxels, and record the vertical topological relationships to provide a basis for subsequent surface analysis.

[0126] Step 4: Analyze the occlusion relationship from top to bottom in the vertical direction, assign the surface value of the lower voxel based on the distance and number of points of the upper voxel; then smooth the surface value through neighborhood propagation, and finally remove the occluded non-bottom voxels and retain the potential surface voxels.

[0127] Step 5: Vertically project the effective voxels onto the horizontal plane to generate a two-dimensional raster image (pixel values ​​are surface values), and mark holes as invalid values; fill small holes through neighborhood information propagation and bilateral filtering to smooth abrupt changes in surface values.

[0128] Step 6: In areas with high surface values, perform attribute constraints and edge refinement (expand to abrupt change points or half-voxels), remove unclosed areas and concave terrain, and retain convex terrain boundaries; connect boundary pixels into closed vector contour lines and simplify them to obtain the non-ground range.

[0129] Step 7: Based on the extracted non-ground contour lines, traverse each point within the point cloud block and use the ray casting method to determine whether a point is inside the contour line; internal points are marked as non-ground points, and external points are marked as ground points, and the results are stored in the attribute field.

[0130] Step 8: Using the centroids of the voxels surrounding the cavity as boundary constraints, fit a smooth surface using the thin-plate spline method, and generate virtual ground points by sampling at DEM resolution; perform boundary fusion and filtering to obtain a complete and continuous ground point cloud.

[0131] Application Example: Using the method of this invention, an application was made to a town in a hilly area. The total point cloud data consisted of 120 million points, with each point cloud block containing approximately 5 million points. After initial segmentation, adaptive merging, and secondary segmentation, the number of points in each point cloud block was evenly distributed between 2.5 million and 7.5 million. A resolution of 0.05m (corresponding to an average point spacing of approximately 0.05m for 400 points / square meter) and a scale of 1:500 were used, resulting in a voxel side length of 0.1m. Subsequently, steps such as vertical occlusion analysis, occlusion voxel removal, and projection planarization were performed to extract non-ground contours such as buildings and trees, and non-ground points were marked. Finally, for voids in the ground points (such as areas occluded by tree canopies), the voxel centroids were used as boundary constraints, and a thin-plate spline method was employed to fit the surface and sample and fill it, generating a digital terrain model (DEM) with a resolution of 0.05m. The entire processing was completed on a regular workstation (Intel(R) Core(TM) i9-10980XE CPU @ 3.00GHz, 32GB memory), taking about 3 minutes, which verifies the efficiency and robustness of the invention in medium-density massive point cloud scenarios.

[0132] Example 2:

[0133] Based on the same inventive concept, this invention also provides a massive point cloud ground point extraction system based on voxel morphology, comprising:

[0134] The block segmentation module is used to segment massive point clouds into blocks, forming several point cloud blocks, and to perform point cloud preprocessing on each point cloud block.

[0135] The voxelization module is used to voxelize each preprocessed point cloud block according to the resolution of the original massive point cloud, forming several three-dimensional voxel meshes.

[0136] The surface value calculation module is used to jointly calculate the surface value of several 3D voxel meshes through vertical occlusion analysis and neighborhood propagation, and to remove occluded non-bottom voxels.

[0137] The projection planarization module is used to vertically project effective voxels onto a horizontal plane to generate a two-dimensional raster image, and holes are marked as invalid values; small holes are filled through neighborhood information propagation and bilateral filtering to smooth abrupt changes in surface values;

[0138] The non-ground contour extraction module is used to perform attribute constraint optimization and fill in missing data on the two-dimensional raster image, and to perform edge refinement operations to extract the effective surface boundary; the effective surface boundary edge pixels are connected into a closed vector contour line through an eight-neighbor tracing algorithm and simplified to obtain the ground contour line;

[0139] The non-ground point marking module is used to traverse each point in the point cloud block according to the ground contour line, and use the ray method to determine whether the point is inside the contour line. The internal points are marked as non-ground points, and the external points are marked as ground points.

[0140] The void virtual filling module is used to fill non-ground points based on markers. It uses the centroid points of voxels around the voids as boundary constraints, fits a smooth surface using the thin plate spline method, and generates virtual ground points by sampling at DEM resolution. Boundary fusion and filtering are then performed to obtain a complete and continuous ground point cloud.

[0141] Example 3:

[0142] Based on the same inventive concept, such as Figure 6 As shown, the present invention also provides an electronic device, which may be a computer device, a microcontroller device, a smart mobile device, etc. The electronic device in this embodiment may include a processor, a memory, a transceiver component, etc. The memory, processor, and transceiver component are connected via a bus; the memory can be used to store executable programs, and an exemplary executable program may include instructions; the processor is used to execute the instructions stored in the memory. The memory can also be used to store data, which can be accessed and / or modified when instructions are executed.

[0143] The processor may be a Central Processing Unit (CPU), or it may be other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. It is the computing core and control core of the terminal, and it is suitable for implementing one or more instructions. Specifically, it is suitable for loading and executing one or more instructions in the storage medium to realize the corresponding method flow or corresponding function, so as to realize the steps of the voxel morphology-based massive point cloud ground point extraction method in the above embodiments.

[0144] Example 4:

[0145] Based on the same inventive concept, this invention also provides a readable storage medium, specifically an electronic device readable storage medium (Memory). This readable storage medium is a memory device within an electronic device used to store programs and data. It is understood that the storage medium here can include both built-in storage media within the electronic device and extended storage media supported by the electronic device. The storage medium provides storage space, which stores the terminal's operating system. Furthermore, this storage space also stores one or more instructions suitable for loading and execution by a processor. These instructions can be one or more executable programs (including program code). It should be noted that the storage medium here can be high-speed RAM or non-volatile memory, such as at least one disk storage device. Loading and executing one or more instructions stored in the storage medium by the processor can implement the steps of the voxel morphology-based massive point cloud ground point extraction method described in the above embodiments.

[0146] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0147] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0148] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0149] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0150] The above are merely embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention are included within the scope of the claims of the present invention pending approval.

Claims

1. A method for extracting ground points from massive point clouds based on voxel morphology, characterized in that, include: The massive point cloud is divided into blocks to form several point cloud blocks, and point cloud preprocessing is performed on each point cloud block. Based on the resolution of the original massive point cloud, each preprocessed point cloud block is voxelized. The voxelized point cloud blocks are calculated by combining vertical occlusion analysis and neighborhood propagation to determine the voxel surface value, and the occluded non-bottom voxels are removed. Based on the voxel surface values, the effective voxels are vertically projected onto the horizontal plane to generate a two-dimensional raster image, and the invalid value regions in the two-dimensional raster image are filled through neighborhood information propagation and bilateral filtering. The two-dimensional raster image is optimized by attribute constraints to fill in missing data, and edge thinning is performed to extract effective surface boundaries. The effective surface boundary edge pixels are connected into a closed vector contour line by an eight-neighbor tracing algorithm and simplified to obtain the ground contour line. Based on the ground contour line, traverse each point in all point cloud blocks, and use the ray method to determine whether the point is inside the contour line. Points inside the contour line are marked as non-ground points, and points outside the contour line are marked as ground points. Based on the marked non-ground points, the voxel centroids around the holes are used as boundary constraints. A thin plate spline algorithm is used to fit a smooth surface and virtual ground points are generated by sampling at DEM resolution. Boundary fusion and filtering are performed on all point cloud blocks to obtain a complete and continuous ground point cloud.

2. The method for extracting ground points from massive point clouds based on voxel morphology according to claim 1, characterized in that, The process of dividing the massive point cloud into blocks to form several point cloud blocks includes: Obtain basic information of massive point clouds, and perform preliminary virtual segmentation of massive point clouds based on the basic information, dividing the original spatial range into several virtual spaces, each virtual space corresponding to a point cloud block; During the virtual segmentation process, the number of points in the current point cloud block is counted in real time. After the virtual segmentation is completed, it is ensured that the number of points in each point cloud block reaches the preset segmentation threshold. For point cloud blocks with fewer than 0.5 times the threshold, they are merged according to the spatial neighborhood topology. For point cloud blocks with more than 1.5 times the threshold, a second spatial segmentation is performed to generate smaller point cloud blocks. The point cloud blocks, which are divided into several virtual spaces, are traversed once. Each point in the massive point cloud is allocated and filled into a unique virtual space and stored according to a preset format.

3. The method for extracting ground points from massive point clouds based on voxel morphology according to claim 1, characterized in that, The point cloud preprocessing for each point cloud block includes: For each point cloud block, noise removal, isolated point removal, and data repair are performed, and adjacent points in overlapping areas are included during processing to ensure continuity between point cloud blocks.

4. The method for extracting ground points from massive point clouds based on voxel morphology according to claim 1, characterized in that, The process of voxelizing each preprocessed point cloud block according to the resolution of the original massive point cloud includes: Based on the resolution of the original point cloud and the task scale, the voxel edge length algorithm is used to determine the voxel edge length, and voxel segmentation is performed based on the voxel edge length to form a number of voxels. Calculate the centroid, number of centroids, and elevation value for each non-empty voxel, and clean up voxels with attribute values ​​of zero, while retaining the voxel spatial topological relationships.

5. The method for extracting ground points from massive point clouds based on voxel morphology according to claim 1, characterized in that, The voxelized point cloud blocks are processed by jointly calculating voxel surface values ​​through vertical occlusion analysis and neighborhood propagation, and occluded non-bottom voxels are removed, including: For a number of three-dimensional voxel grids, traverse each voxel column from top to bottom along the vertical direction; For voxels in the same column, the surface value of the lower voxel is assigned based on the distance and number of points of the upper voxel; if the upper voxel exists, it will have an occlusion effect on the lower voxel. The surface values ​​are smoothed by propagating through the neighborhood, and finally the occluded non-bottom voxels are removed, while the potential surface voxels are preserved.

6. The method for extracting ground points from massive point clouds based on voxel morphology according to claim 1, characterized in that, The process of vertically projecting effective voxels onto a horizontal plane based on the voxel surface values ​​to generate a two-dimensional raster image, and filling in invalid value regions in the two-dimensional raster image through neighborhood information propagation and bilateral filtering, includes: All effective voxels are vertically projected onto a horizontal plane to generate a two-dimensional raster image; the value of each pixel is assigned the corresponding voxel's surface value; During projection, if there are no valid voxels at a certain plane position, the corresponding pixel is marked as invalid. Invalid value regions are extracted. When the range of invalid value regions is less than 1 / 2 of the minimum sensitive area of ​​ground features, low values ​​are filled by the boundary values ​​around the void range. The surface values ​​of voxels in the current voxel spatial neighborhood are taken, and the surface values ​​of the current voxels are updated using bilateral filtering to enhance the neighborhood edge information, filter and smooth flat areas of the same nature, and consider the coherence of the information of the boundary link areas of different point cloud blocks.

7. The method for extracting ground points from massive point clouds based on voxel morphology according to claim 1, characterized in that, The process of performing attribute constraint optimization to fill missing data in the two-dimensional raster image, and performing edge refinement to extract effective surface boundaries; connecting the effective surface boundary edge pixels into a closed ground contour line using an eight-neighbor tracing algorithm, includes: For the two-dimensional raster image, attribute constraints are performed based on synchronously acquired image data and after heterogeneous data processing and alignment. For the portion of the land surface value in the attribute-constrained 2D raster image that is higher than a preset threshold, edge refinement is performed by using a weighted average of the neighborhood spatial location and the land surface value to extract the effective surface boundary. Remove the unclosed areas in the effective table boundary, and then determine whether the vertical coordinate value within the closed effective table boundary is greater than the corresponding boundary value. If it is greater, retain it; otherwise, delete the effective table boundary. The processed effective surface boundary edge pixels are connected into a closed vector contour line using an eight-neighbor tracing algorithm, and then simplified using Douglas-Pockmark polyline to remove redundant nodes, generating a smooth and compact ground contour line.

8. The method for extracting ground points from massive point clouds based on voxel morphology according to claim 1, characterized in that, The marker-based non-ground points are constrained by the centroids of the voxels surrounding the holes. A smooth surface is fitted using the thin-plate spline method, and virtual ground points are generated by sampling at DEM resolution. Boundary fusion and filtering are performed on all point cloud blocks to obtain a complete and continuous ground point cloud, including: Using the set of voxel centroid points on the cavity boundary as constraints, the thin plate spline algorithm is used to fit the marked non-ground points to a smooth terrain surface; Inside the cavity, samples are taken at equal intervals according to the DEM resolution, and the elevation value of each sampling point is calculated to generate virtual ground points. Based on the virtual ground points, sampling and filling are performed, and boundary fusion and filtering are applied to the filled area to ensure a continuous transition between the virtual points and the original ground points, thereby obtaining a complete and continuous ground point cloud.

9. The method for extracting ground points from massive point clouds based on voxel morphology according to claim 8, characterized in that, The formula for calculating the elevation value is: in, This is the elevation value. For affine coefficients, The number of boundary points. For index value, To sample points at equal intervals according to DEM resolution, For the coordinates of the control points, As the base weight, As basis functions, The distance between the two points is the Euclidean distance.

10. A system for extracting ground points from massive point clouds based on voxel morphology, characterized in that, include: The block segmentation module is used to segment massive point clouds into blocks, forming several point cloud blocks, and to perform point cloud preprocessing on each point cloud block. The voxelization module is used to voxelize each preprocessed point cloud block according to the resolution of the original massive point cloud. The surface value calculation module is used to jointly calculate the voxel surface value of point cloud blocks after voxelization through vertical occlusion analysis and neighborhood propagation, and remove occluded non-bottom voxels. The projection planarization module is used to vertically project the effective voxels onto the horizontal plane based on the voxel surface values ​​to generate a two-dimensional raster image, and fill the invalid value regions in the two-dimensional raster image through neighborhood information propagation and bilateral filtering. The ground contour extraction module is used to perform attribute constraint optimization and fill in missing data on the two-dimensional raster image, and to perform edge refinement operations to extract effective surface boundaries; the effective surface boundary edge pixels are connected into closed vector contour lines through an eight-neighbor tracing algorithm and simplified to obtain the ground contour lines; The non-ground point marking module is used to traverse each point in all point cloud blocks according to the ground contour line, and use the ray method to determine whether the point is inside the contour line. The internal points are marked as non-ground points, and the external points are marked as ground points. The virtual hole filling module is used to fill non-ground points based on markers. It uses the centroid points of voxels around the holes as boundary constraints, fits smooth surfaces using a thin plate spline algorithm, and generates virtual ground points by sampling at DEM resolution. It then performs boundary fusion and filtering on all point cloud blocks to obtain a complete and continuous ground point cloud.

Images