A method and apparatus for converting massive laser point clouds into 3DTiles using LAS.

Abstract

This invention relates to the field of point cloud data processing technology, and provides a method and apparatus for converting massive laser point clouds (LAS) into 3DTiles. The invention determines the spatial range and edge information of the point cloud based on the input point cloud file; based on the point cloud edge information and the point cloud spatial range, it spatially divides the point cloud data in the point cloud file to obtain multi-level data; it calculates the point cloud reading and rendering threshold for each level in the multi-level data to generate a point cloud index file; and it constructs a tile set based on the point cloud index file and the multi-level data, solving the problems of slow slicing or missing slice data in existing technologies for massive laser point clouds (LAS).

Description

Technical Field

[0001] This invention relates to the field of point cloud data processing technology, and in particular to a method and apparatus for converting massive laser point clouds (LAS) into 3DTiles. Background Technology

[0002] Airborne Light Detection and Ranging (LiDAR) emits lasers that can partially penetrate the obstruction of trees, directly acquiring high-precision three-dimensional topographic data of the Earth's surface. This overcomes the shortcomings and deficiencies of traditional remote sensing and geodetic methods in acquiring elevation data and rapid automated processing. LiDAR's high spatial resolution, strong anti-interference capabilities, and high degree of automation enable it not only to rapidly acquire elevation data but also to demonstrate excellent development prospects and strong competitiveness in many fields such as topographic mapping, urban 3D modeling, environmental monitoring, engineering construction, earth sciences, and planetary sciences. In recent years, it has also been widely used in the identification and monitoring of landslide geological hazards, achieving a series of results and new advancements.

[0003] The raw point cloud laser file format (LAS) data of airborne lidar often has tens of millions or hundreds of millions of points. Therefore, if the massive point cloud data is not sliced ​​(i.e. converted to 3DTiles), it is difficult to visualize and render it on the web. 3DTiles refers to the three-dimensional slice format, which is the form of publishing and applying LAS data.

[0004] Existing point cloud slicing tools suffer from problems such as only partially successful slicing or slicing task failure when slicing massive laser point cloud (LAS) data, and slow slicing speed, which seriously affects usability.

[0005] Therefore, overcoming the shortcomings of the existing technology is an urgent problem to be solved in this technical field. Summary of the Invention

[0006] The technical problem to be solved by the present invention is to provide a method and apparatus for converting massive laser point clouds (LAS) into 3DTiles. The purpose is to solve the problems of slow slicing or missing slice data in the prior art.

[0007] The present invention adopts the following technical solution: In a first aspect, the present invention provides a method for converting massive laser point clouds (LAS) into 3DTiles, comprising: Determine the spatial extent and edge information of the point cloud based on the input point cloud file; Based on the point cloud edge information and the point cloud spatial range, the point cloud data of the point cloud file is spatially divided to obtain multi-level data; Calculate the point cloud reading and rendering threshold for each level in the multi-level data to generate a point cloud index file; A tile set is constructed based on the point cloud index file and the multi-level data.

[0008] Furthermore, determining the point cloud spatial range and point cloud edge information based on the input point cloud file includes: Calculate the centroid of the point cloud data; perform mean-reduction processing on each point in the point cloud data according to the centroid to obtain the first data; Calculate the covariance matrix of the first data to obtain the principal direction of the point cloud; correct the principal direction of the point cloud so that it coincides with the coordinate axis direction to obtain the second data. Calculate the point cloud edge information of the second data.

[0009] Furthermore, the point cloud edge information includes the minimum bounding box and the point cloud boundary point set; The point cloud edge information used to calculate the second data includes: Obtain one of the current points of the second data, determine multiple nearest neighbor points of the current point, and obtain a local surface point set; The micro-tangent plane of the local surface point set is obtained by fitting; the angle between the current point and the corresponding projection points of the multiple nearest neighbor points is calculated using the projection of the micro-tangent plane as the reference plane. When the maximum value of the included angles among the plurality of nearest neighbor points is greater than the angle threshold, the current point is determined as a boundary point; Traverse all points in the second data, iteratively taking each point as the current point to obtain multiple boundary points, and use the multiple boundary points to construct a point cloud boundary point set; The minimum bounding box of the second data is solved using the point cloud boundary point set.

[0010] Furthermore, determining the point cloud spatial range and point cloud edge information based on the input point cloud file also includes: The point cloud file is parsed to obtain point cloud attributes; wherein, the point cloud attributes include point cloud bounding boxes; The global bounding box of the second data is calculated using the point cloud bounding box to obtain the spatial extent of the point cloud.

[0011] Furthermore, the spatial partitioning of the point cloud data in the point cloud file based on the point cloud edge information and the point cloud spatial range to obtain multi-level data includes: Using the spatial range of the point cloud as the root tile, calculate the center point and geometric error of the root tile; Obtain the second data corresponding to the point cloud edge information; obtain the target point cloud in the second data, and determine the depth parameter using the geometric error and the number of points contained in the target point cloud; Using the center point as the reference point for octree segmentation, and according to the depth parameters, the point cloud edge information is used to iteratively construct the octree to complete the spatial division along the outer contour of the point cloud and obtain multi-level data.

[0012] Furthermore, the point cloud edge information includes the minimum bounding box and the point cloud boundary point set; The step of iteratively constructing an octree using the point cloud edge information according to the depth parameters includes: The minimum bounding box is uniformly divided into multiple sub-bounding boxes along each directional axis, and the spatial information of all current sub-bounding boxes is determined. The sub-bounding boxes are used as child nodes of the octree, and the point cloud within the coverage of the spatial information is divided into the corresponding sub-bounding boxes to obtain the third data. Project the point cloud boundary point set onto the plane of each directional axis respectively; Determine the projection points of the boundary points in the point cloud boundary point set corresponding to different planes, and determine the maximum and minimum values ​​of the projection points on each directional axis in each child bounding box to obtain the range of the boundary points in each direction. The child bounding boxes corresponding to the minimum value are reorganized from the range domains in each direction. The child bounding boxes corresponding to the minimum value are used as parent nodes to generate new child nodes and child bounding boxes for the next layer. Iteratively reorganize the child bounding boxes corresponding to the minimum value until the depth of the octree satisfies the depth parameter, thus obtaining the octree.

[0013] Furthermore, the method also includes: A GPU-based KD-tree spatial index is constructed for the multi-level data; Using the KD tree spatial index, for each current point in the multi-level data, determine the covariance matrix of the neighborhood of each point; use the local covariance matrix to perform eigenvalue decomposition to obtain the degree of change of each directional axis, and sort the eigenvalues ​​in descending order; Initial key points are identified using the feature value ratio condition based on the aforementioned feature values; The specified feature values ​​are used as response intensity indicators, and the initial key points with the strongest response are retained to obtain the target key point set.

[0014] Furthermore, the feature value ratio condition is as follows: the change in the multi-level data in the second principal direction is greater than the change in the multi-level data in the first principal direction relative to the change in the multi-level data in the first principal direction, and the change in the multi-level data in the third principal direction is greater than the change in the second principal direction relative to the change in the second principal direction.

[0015] Secondly, the present invention also provides an apparatus for converting massive laser point clouds (LAS) into 3DTiles, for implementing the method for converting massive laser point clouds (LAS) into 3DTiles as described in the first aspect, the apparatus comprising: At least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor for performing the method for converting massive laser point cloud LAS to 3DTiles as described in the first aspect.

[0016] Thirdly, the present invention also provides a non-volatile computer storage medium storing computer-executable instructions, which are executed by one or more processors to perform the method for converting massive laser point cloud LAS to 3DTiles as described in the first aspect.

[0017] This invention determines the spatial range and edge information of the input point cloud file, and spatially partitions the point cloud data accordingly to obtain multi-level data. Through multi-level spatial partitioning, the point cloud space is divided, ensuring that each point in the point cloud data is processed according to its spatial range and edge information, thus maintaining the structural integrity of the point cloud data. The point cloud reading and rendering threshold for each level in the multi-level data is calculated and used as an index to generate a point cloud index file and construct a tile set. By calculating the corresponding point cloud reading and rendering threshold at each level and generating the index accordingly, data loss in sparse areas or computational redundancy in dense areas is avoided. Furthermore, by calculating the point cloud reading and rendering threshold for each level to establish a spatial index, the relationship between point cloud location and level is recorded, thus supporting rapid retrieval. This solves the problems of slow processing or missing slice data in existing technologies for massive laser point cloud (LAS) point cloud slicing. Attached Figure Description

[0018] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments of the present invention will be briefly described below. Obviously, the drawings described below are merely some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without any creative effort.

[0019] Figure 1This is a flowchart illustrating a method for converting massive laser point clouds (LAS) into 3DTiles according to an embodiment of the present invention. Figure 2 This is a flowchart illustrating step 10 provided in an embodiment of the present invention; Figure 3 This is a flowchart illustrating step 103 provided in an embodiment of the present invention; Figure 4 This is a flowchart illustrating another step 10 provided in an embodiment of the present invention; Figure 5 This is a schematic diagram illustrating a specific example of the spatial range of point cloud data provided in an embodiment of the present invention; Figure 6 This is a flowchart illustrating step 20 provided in an embodiment of the present invention; Figure 7 This is a schematic diagram illustrating a specific example of a box type provided in an embodiment of the present invention; Figure 8 This is a schematic diagram illustrating a specific example of a box-type octree segmentation provided in an embodiment of the present invention; Figure 9 This is a schematic diagram illustrating a specific instance of a region type provided in an embodiment of the present invention; Figure 10 This is a schematic diagram illustrating a specific example of a region-type octree segmentation provided in an embodiment of the present invention; Figure 11 This is a schematic diagram of a process for extracting key points provided in an embodiment of the present invention; Figure 12 This is a schematic diagram illustrating a specific example of LAS point cloud data conversion to 3DTiles tile data preview provided in an embodiment of the present invention; Figure 13 This is a schematic diagram illustrating a specific example of LAS point cloud data conversion to 3DTiles tile data preview provided in an embodiment of the present invention; Figure 14 This is a schematic diagram of the architecture of a device for converting massive laser point clouds into 3DTiles according to an embodiment of the present invention. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0021] In the description of this invention, the terms "inner", "outer", "longitudinal", "lateral", "upper", "lower", "top", "bottom", etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are only for the convenience of describing this invention and do not require that this invention must be constructed and operated in a specific orientation. Therefore, they should not be construed as limiting this invention.

[0022] In this invention, the terms "first," "second," etc., 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. Therefore, a feature defined with "first," "second," etc., may explicitly or implicitly include one or more of that feature. In the description of this application, unless otherwise stated, "a plurality of" means two or more.

[0023] In this application, unless otherwise expressly specified and limited, the term "connection" should be interpreted broadly. For example, "connection" can be a fixed connection, a detachable connection, or an integral part; it can be a direct connection or an indirect connection through an intermediate medium. Furthermore, the term "coupled" can refer to an electrical connection that enables signal transmission.

[0024] Furthermore, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0025] Example 1: Existing commercial tools (such as cesiumLab) have several limitations in their LAS parsing modules: Firstly, the format compatibility is narrow. Although it supports parsing the basic fields of the mainstream LAS 1.2 version, it is not fully compatible with the extended fields added in version 1.4, such as the scanning direction angle, the Global Positioning System (GPS) timestamp, and the category identifier extension. When processing LAS data containing such extended fields, some data blocks cannot be sliced ​​normally due to field parsing misalignment, resulting in only some slices being successful.

[0026] Secondly, parsing the compressed point cloud data format (LAS zip compressed format, abbreviated as LAZ compressed format) depends on third-party plugins. If the plugin version is incompatible or not installed in advance, it will directly trigger a parsing error, causing the entire tiling task to fail.

[0027] Third, the massive data parsing logic is rudimentary. It adopts the "overall loading - point-by-point parsing" mode. When processing gigabyte (GB) or terabyte (TB) level LAS files, it is very easy for the slicing process to be interrupted due to memory overflow or input / output read blocking. Either the task fails directly or only a small part of the data fragments are sliced.

[0028] Fourth, the resource monitoring and threshold settings are unreasonable. The CPU and memory usage thresholds of CesiumLab are preset according to the conventional geographic data processing scenario and are not adapted to the characteristics of LAS point cloud tiling of "high memory usage and high input and output throughput". When processing massive amounts of data, it is very easy for the system to be forcibly terminated due to memory usage exceeding the threshold, or for the tiling process to stagnate due to excessive input and output read and write pressure. Ultimately, this manifests as tiling task failure or some data tiling timeouts.

[0029] To solve the above problems, such as Figure 1 As shown, this embodiment of the invention provides a method for converting massive laser point clouds (LAS) into 3DTiles, including: Step 10: Determine the spatial extent and edge information of the point cloud based on the input point cloud file.

[0030] First, the applicable scenarios of the embodiments of the present invention will be explained: After scanning using an airborne LiDAR system, the raw data is acquired and stored in LAS format. Then, the LAS format point cloud data needs to be converted into 3D Tiles format for efficient and smooth 3D visualization on a webpage or client. This invention's method for converting massive laser point cloud data from LAS to 3D Tiles is designed for applications with extremely large amounts of raw data. It converts LAS format point cloud data into 3D Tiles format to obtain a complete point cloud model for subsequent visualization on a webpage or client. LAS format is a standard binary file format for storing point cloud data, recording the 3D coordinates (X, Y, Z) of each point, as well as other attributes such as intensity, echo count, and color information.

[0031] In one embodiment, the laspy library in the Python open-source library can be used to read the LAS file (i.e., the input point cloud file), parse the metadata of the input LAS file (such as boundary information, point density, coordinate system), and count the point cloud attributes. Among them, the point cloud attributes may include: (1) the point cloud bounding box, that is, the minimum and maximum coordinates in the three directions of X, Y and Z axes, which can be: minX, minY, minZ, maxX, maxY, maxZ; (2) the number of points; (3) the point attributes, such as the coordinates, color, classification value, intensity, etc. of the points; (4) the coordinate system, such as the code of the European Petroleum Survey Group (EPSG).

[0032] In an alternative embodiment, point data can also be read in chunks to avoid memory overflow: "with laspy.open('input.las') as f: # Get metadata header = f.header point_format = header.point_format # Read point data block by block for points in f.chunk_iterator(2_000_000): # 2 million points per chunk x = points. y = points.y z = points. # Handling point attributes (such as color) if 'red' in point_format.dimension_names: colors = np.vstack((points.red, points.green, points.blue)).transpose()”.

[0033] Among them, the point cloud spatial range is the spatial range of the entire point cloud data in the point cloud file; the point cloud edge information is the key information describing the shape of the point cloud; specific examples of point cloud spatial range and point cloud edge information will be provided below, and will not be elaborated here.

[0034] Step 20: Based on the point cloud edge information and the point cloud spatial range, spatially divide the point cloud data of the point cloud file to obtain multi-level data.

[0035] Among them, multi-level data refers to point cloud data at different levels. In this embodiment of the invention, the complex and disordered three-dimensional point cloud data is decomposed into subsets with similar characteristics by spatially partitioning the point cloud data, which facilitates subsequent processing and analysis; after segmentation, the point clouds in different regions can be processed independently.

[0036] Step 30: Calculate the point cloud reading and rendering threshold for each level in the multi-level data to generate a point cloud index file.

[0037] In one embodiment, the expression for the point cloud read rendering threshold is: ; in, This refers to the geometric error corresponding to the spatial range of the point cloud; a specific example of geometric error is provided below, which will not be elaborated upon here.

[0038] For each point in the point cloud data, calculate the geometric error of that point, and then calculate the corresponding point cloud reading and rendering threshold for each level in the multi-level data; where n is the level, representing the level of the level.

[0039] In one embodiment, the path of the point cloud folder corresponding to each level of multi-level data and the point cloud rendering reading threshold at each level are read, and a point cloud index file is generated accordingly.

[0040] Step 40: Construct a tile set based on the point cloud index file and the multi-level data.

[0041] The specific method for constructing the tile set based on the point cloud index file and multi-level data shall be determined by those skilled in the art based on the specific application scenario. In one embodiment, the corresponding index in the point cloud index file is used to index the point cloud data in the multi-level data, thereby constructing a complete point cloud model (i.e., a 3DTiles tile set) from the corresponding point cloud data.

[0042] This invention determines the spatial range and edge information of the input point cloud file, and spatially partitions the point cloud data accordingly to obtain multi-level data. Through multi-level spatial partitioning, the point cloud space is divided, ensuring that each point in the point cloud data is processed according to its spatial range and edge information, thus maintaining the structural integrity of the point cloud data. The point cloud reading and rendering threshold for each level in the multi-level data is calculated and used as an index to generate a point cloud index file and construct a tile set. By calculating the corresponding point cloud reading and rendering threshold at each level and generating the index accordingly, data loss in sparse areas or computational redundancy in dense areas is avoided. Furthermore, by calculating the point cloud reading and rendering threshold for each level to establish a spatial index, the relationship between point cloud location and level is recorded, thus supporting rapid retrieval. This solves the problems of slow processing or missing slice data in existing technologies for massive laser point cloud (LAS) point cloud slicing.

[0043] In one embodiment, such as Figure 2 As shown, in order to obtain point cloud edge information, step 10 includes: Step 101a: Calculate the centroid of the point cloud data; perform mean-reduction processing on each point in the point cloud data according to the centroid to obtain the first data.

[0044] Calculate the centroid of point cloud data P In one embodiment, the formula for calculating the centroid is: ; ; ; in, Let P be the size of the point cloud data. Let be the x-coordinates of each point in the point cloud data P. Let be the y-coordinates of each point in the point cloud data P. Let be the z-coordinates of each point in the point cloud data P.

[0045] The specific method for demeaning based on the centroid is determined by those skilled in the art according to the specific application scenario. In one embodiment, demeaning is performed on each point in the point cloud data P along the x, y, and z axes.

[0046] Step 102a: Calculate the covariance matrix of the first data to obtain the principal direction of the point cloud; correct the principal direction of the point cloud so that it coincides with the coordinate axis direction to obtain the second data.

[0047] For the first data obtained after mean reduction, the covariance matrix of the point cloud data P is calculated to obtain the principal direction of the point cloud. The specific method of correcting the principal direction of the point cloud to make it coincide with the coordinate axis direction shall be determined by those skilled in the art based on the specific scenario.

[0048] Step 103a: Calculate the point cloud edge information of the second data.

[0049] In one embodiment, the point cloud edge information includes a minimum bounding box and a set of point cloud boundary points.

[0050] The point cloud boundary point set refers to the set of points located on the outer contour or boundary of a point cloud. In one embodiment, this is achieved by extracting the outermost points of the point cloud surface. The point cloud boundary point set describes the shape boundary of the point cloud, but it is not a regular geometric shape (such as a cuboid). The minimum bounding box (MBB) is a cuboid (in three-dimensional space) with the smallest volume that can completely contain all points of the point cloud. The minimum bounding box is found to minimize its volume by rotation to find the optimal orientation. The minimum bounding box is generally a regular cuboid with definite dimensions and orientation. While the volume of the minimum bounding box is the smallest, its boundary does not necessarily coincide completely with the point cloud boundary point set; it is only an approximate enclosure.

[0051] like Figure 3 As shown, step 103a includes: Step 1031: Obtain one of the current points of the second data, determine multiple nearest neighbor points of the current point, and obtain a local surface point set.

[0052] The second data is obtained through step 102a above. In one embodiment, a point Pi is randomly selected from the second data P, and the k nearest neighbor search algorithm is used to find the k nearest neighbor points of point Pi to form a local surface point set X for forming the local surface.

[0053] Step 1032: Fit the micro-tangent plane of the local surface point set; using the projection of the micro-tangent plane as a reference plane, calculate the angle between the lines connecting the current point and the corresponding projection points of the multiple nearest neighbor points.

[0054] The micro-tangent plane of the local surface point set X is obtained by fitting; in one embodiment, the micro-tangent plane of the point set X is fitted using the least squares method. Taking the projection of the local surface point set X onto the micro-tangent plane as a reference plane, the angle α between the lines connecting point Pi and the corresponding projection points of its k nearest neighbors is calculated.

[0055] Step 1033: When the maximum value of the included angles among the multiple nearest neighbor points is greater than the angle threshold, the current point is determined as the boundary point.

[0056] The angle threshold is determined by those skilled in the art based on the specific application scenario. When the maximum value of the included angle α is greater than the set angle threshold ε, the current point is determined to be a boundary point of the point cloud.

[0057] Step 1034: Traverse all points in the second data, iteratively taking each point as the current point to obtain multiple boundary points, and use the multiple boundary points to construct a point cloud boundary point set.

[0058] In one embodiment, all points in the point cloud data P are traversed (i.e., all points in the point cloud data P are treated as points Pi and the above step 1033 is executed) to extract and record all boundary points, thereby obtaining the cloud boundary point set S.

[0059] Step 1035: Use the point cloud boundary point set to solve for the minimum bounding box of the second data.

[0060] The point cloud boundary point set is the basis for minimum bounding box calculation: when calculating the minimum bounding box, the point cloud boundary point set is extracted first, and the minimum volume cuboid is calculated based on the point cloud boundary point set to obtain the minimum bounding box.

[0061] Based on this, in order to calculate the spatial extent of the point cloud, in one embodiment, such as Figure 4 As shown, step 10 further includes: Step 101b: Parse the point cloud file to obtain point cloud attributes; wherein, the point cloud attributes include point cloud bounding boxes.

[0062] In one embodiment, the point cloud bounding box is: the minimum value minX and the maximum value maxX in the X-axis direction, the minimum value minY and the maximum value maxY in the Y-axis direction, and the minimum value minZ and the maximum value maxZ in the Z-axis direction.

[0063] Step 102b: Calculate the global bounding box of the second data using the point cloud bounding box to obtain the spatial extent of the point cloud.

[0064] Calculate the global bounding box based on the point cloud bounding box, and use the global bounding box as the spatial extent of the point cloud; a specific example of a point cloud spatial extent is shown below. Figure 5 As shown.

[0065] In one embodiment, such as Figure 6 As shown, step 20 includes: Step 201: Using the spatial range of the point cloud as the root tile, calculate the center point and geometric error of the root tile.

[0066] The calculation of the geometric error of the root tile's center point and various points within the root tile, and the specific method for calculating the center point and geometric error, are determined by those skilled in the art. In one embodiment, the formula for calculating the geometric error of the root tile is as follows: geometric_error = dis / height Where dis is the distance from the viewpoint to the point cloud in meters; height is the pixel size of the current viewpoint window.

[0067] Step 202: Obtain the second data corresponding to the point cloud edge information; obtain the target point cloud in the second data, and determine the depth parameter using the geometric error and the number of points contained in the target point cloud.

[0068] This invention uses a global bounding box to define the spatial extent of the octree, uses the center point as the reference point for octree segmentation, and derives the octree depth parameters based on geometric errors. Based on point cloud edge information, the complete point cloud data is spatially divided along the outer contour of the point cloud to form point cloud data of different levels. Specifically, the octree depth parameters are calculated according to the size N of the point cloud data P. octree depth parameters The calculation formula is as follows: ; in, The target point cloud contains the number of points. and These are preset experimental constants, determined by those skilled in the art based on the specific application scenario; in one embodiment, Set to 1000000 Set it to 4.

[0069] Step 203: Using the center point as the reference point for octree segmentation, and according to the depth parameters, iteratively construct the octree using the point cloud edge information to complete the spatial division along the outer contour of the point cloud and obtain multi-level data.

[0070] Among them, multi-level data refers to point cloud data at different levels after spatial division.

[0071] In one embodiment, the point cloud edge information includes a minimum bounding box and a set of point cloud boundary points; in step 203, iteratively constructing an octree using the point cloud edge information according to the depth parameter includes: Step 2031: Divide the minimum bounding box evenly into multiple child bounding boxes along each directional axis, and determine the spatial information of all current child bounding boxes; use the child bounding boxes as child nodes of the octree, and divide the point cloud within the coverage of the spatial information into the corresponding child bounding boxes to obtain the third data.

[0072] In one embodiment, the minimum bounding box of point cloud data P is uniformly divided into eight sub-bounding boxes along the X, Y and Z axes. The spatial information of all current sub-bounding boxes is recorded. The sub-bounding boxes are used as sub-nodes to divide the point cloud within the coverage area into the corresponding sub-bounding boxes. The data processed by the above steps is the third data.

[0073] Step 2032: Project the point cloud boundary point set onto the plane of each directional axis.

[0074] In one embodiment, the orientation axis may include the X, Y, and Z planes, and the point cloud boundary point set S is projected onto the XOY plane, XOZ plane, and YOZ plane, respectively.

[0075] Step 2033: Determine the projection points of the boundary points in the point cloud boundary point set corresponding to different planes, and determine the maximum and minimum values ​​of the projection points on each directional axis in each child bounding box to obtain the range of the boundary points in each direction.

[0076] In one embodiment, the projection points of boundary point Si corresponding to different planes are set as Si. The maximum and minimum values ​​of the corresponding projection points in the X, Y, and Z axes are found in each child-level bounding box to obtain the range of point cloud boundary points in the three main directions. , and ,in, This indicates finding the maximum value of the corresponding projection point along the X-axis in each child bounding box. This indicates finding the minimum value along the X-axis for the corresponding projection point within each child bounding box. This indicates finding the maximum value of the corresponding projection point along the Y-axis in each child bounding box. This indicates finding the minimum value along the Y-axis for the corresponding projection point within each child bounding box. This indicates finding the maximum value of the corresponding projection point along the Z-axis in each child bounding box. This indicates finding the minimum value of the corresponding projection point along the Z-axis in each child bounding box; based on the above three ranges (i.e., , and Reorganize the minimum bounding box of this sub-level.

[0077] Step 2034: Reorganize the child bounding boxes corresponding to the minimum values ​​of the range domains in each direction, and use the child bounding boxes corresponding to the minimum values ​​as parent nodes to generate new child nodes and child bounding boxes for the next layer.

[0078] In one embodiment, following step 2034, the child minimum bounding boxes are iteratively reconstructed based on the point cloud boundary information until the depth parameters of the octree are constructed. Complete the spatial division along the outer contour of the point cloud.

[0079] Step 2035: Iteratively reorganize the child bounding boxes corresponding to the minimum value until the depth of the octree satisfies the depth parameter, thus obtaining the octree.

[0080] In one embodiment, following step 2034, the child minimum bounding boxes are iteratively reconstructed based on the point cloud boundary information until the depth parameters of the octree are constructed. Complete the spatial division along the outer contour of the point cloud.

[0081] This invention uses an octree to hierarchically refine the point cloud data in a point cloud file: the structure is subdivided level by level based on the octree, thereby achieving the effect of prioritizing the processing of high-density areas and reducing computational complexity.

[0082] In practical applications, 3DTiles supports several spatial range cuboid segmentation methods: box type, region type, and sphere type. For the box type, in one embodiment, a specific implementation of the original box is... Figure 7 As shown, a specific implementation of octree partitioning is as follows: Figure 8 As shown. For the region type, a specific implementation of the original region is as follows: Figure 9 As shown; a specific implementation of octree partitioning is shown below. Figure 10 As shown. This embodiment of the invention does not currently support the sphere type because this type cannot be recursively partitioned into quadtrees or octrees.

[0083] For TB-level laser point cloud data, this invention provides an improved Intrinsic Shape Signatures (ISS) algorithm for keypoint extraction. Specifically, it utilizes a K-Dimensional (KD) tree on a Graphics Processing Unit (GPU) for acceleration, while maintaining the geometric representativeness of the keypoints to the original data. In one embodiment, such as... Figure 11 As shown, the method further includes: Step 501: Construct a GPU-based KD tree spatial index for the multi-level data.

[0084] The specific method for constructing the KD-tree spatial index of multi-level data implemented on the GPU shall be determined by those skilled in the art based on the specific use case.

[0085] Step 502: Using the KD tree spatial index, determine the covariance matrix of the neighborhood of each current point in the multi-level data; use the local covariance matrix to perform eigenvalue decomposition to obtain the degree of change of each directional axis, and arrange the eigenvalues ​​in descending order.

[0086] The core of extracting key points from laser point clouds in this invention lies in identifying points with anisotropic geometric features through eigenvalue analysis of the local covariance matrix.

[0087] Step 503: Identify the initial key points using the feature value ratio condition according to the feature values.

[0088] In one embodiment, the feature value ratio condition is: the change in the multi-level data in the second principal direction relative to the change in the first principal direction is greater than a first threshold, and the change in the multi-level data in the third principal direction relative to the change in the second principal direction is greater than a second threshold. In one embodiment, when the axes are X, Y, and Z, the first principal direction is the X-axis principal direction, the second principal direction is the Y-axis principal direction, and the third principal direction is the Z-axis principal direction; wherein, the first and second thresholds are determined by those skilled in the art based on the specific application scenario.

[0089] Step 504: Use the specified feature values ​​as response intensity indicators, retain the initial key points with the strongest response, and obtain the target key point set.

[0090] The improved ISS algorithm of this invention can adaptively extract key points reflecting the local shape features of point clouds, balancing geometric saliency and spatial uniformity. In practical applications, parameters need to be adjusted based on data density and scene requirements to balance the number of key points and feature integrity. This invention can efficiently process TB-level laser point cloud data; by adjusting the salient_radius parameter of the ISS algorithm and performing secondary sampling on low-density regions, it solves the problem of uneven key point distribution and maintains the representativeness of the geometric features of the key points to the original data.

[0091] A specific example of the above steps is as follows: import numpy as np from sklearn.neighbors import KDTree def compute_covariance_matrix(points, neighbors_indices): """ Calculate the covariance matrix of the point cloud :param points: Input point cloud :param neighbors_indices: Neighborhood point indices :return: Covariance matrix """ # Calculate the neighborhood point set neighbors = points[neighbors_indices] # Calculate the centroid centroid = np.mean(neighbors, axis=0) # Decentralization centered_neighbors = neighbors - centroid # Calculate the covariance matrix covariance_matrix = np.dot(centered_neighbors.T, centered_neighbors) / len(neighbors) return covariance_matrix def iss_keypoint_detection(points, salient_radius=0.1, non_max_radius=0.2, gamma_21=0.7, gamma_32=0.7, min_neighbors=5): """ Keypoint detection using the ISS algorithm :param points: Input point cloud, Nx3 NumPy array :param salient_radius: The neighborhood radius used to calculate the covariance matrix. :param non_max_radius: The radius of non-maximum suppression. :param gamma_21: Eigenvalue ratio threshold γ21 :param gamma_32: Eigenvalue ratio threshold γ32 :param min_neighbors: Minimum number of neighboring nodes :return: Key point index list """ # Constructing a KD-tree for fast neighborhood search tree = KDTree(points) # Initialize variables n_points = len(points) `eigenvalues ​​= np.zeros((n_points, 3))` # Stores the three eigenvalues ​​for each point (i.e., in the X, Y, and Z axes). keypoint_indices = [] # Step 1: Calculate the covariance matrix and its eigenvalues ​​for each point. for i in range(n_points): # Query Neighborhood neighbors_indices = tree.query_radius([points[i]], r=salient_radius)[0] # Ensure there are enough neighboring points if len(neighbors_indices) <min_neighbors: continue # Calculate the covariance matrix cov_matrix = compute_covariance_matrix(points, neighbors_indices) # Calculate and sort the eigenvalues evals, _ = np.linalg.eigh(cov_matrix) evals = np.sort(evals)[::-1] # Sort from largest to smallest eigenvalues[i] = evals # Step 2: Apply eigenvalue ratio criteria to filter key points valid_indices = np.where( (eigenvalues[:, 1]>gamma_21eigenvalues[:, 0])& (eigenvalues[:, 2]>gamma_32eigenvalues[:, 1]) [0] # Step 3: Calculate the response value (using the third characteristic value as the response value) responses = eigenvalues[valid_indices, 2] # Step 4: Non-maximum suppression if len(responses) > 0: # Create a mapping between response values ​​and indexes response_indices = np.array(list(zip(responses, valid_indices)),dtype=[('response', float), ('index', int)]) # Sort by response value in descending order response_indices = np.sort(response_indices, order='response')[::-1] # Nonmaximum suppression selected_indices = [] for resp, idx in response_indices: # Check if it is within the non-maximum radius of the selected key point. is_max = True for s_idx in selected_indices: dist = np.linalg.norm(points[idx] - points[s_idx]) if dist <non_max_radius: is_max = False break if is_max: selected_indices.append(idx) keypoint_indices = selected_indices return keypoint_indices # Usage Example if __name__ == "__main__": # Generate example point cloud np.random.seed(42) points = np.random.rand(1000, 3) # 1000 random points # Extracting key points keypoint_indices = iss_keypoint_detection( points, salient_radius=0.1, non_max_radius=0.2 )” Then, in one embodiment, in step 30, the present invention employs voxel downsampling to thin out point cloud data at different levels in multi-level data, calculates the point cloud reading and rendering thresholds for different levels, and generates a point cloud index file. First, based on the average spacing of the point cloud data P... The maximum dimension D is determined using the following formula: ; in, This indicates finding the maximum value of the corresponding projection point along the X-axis in each child bounding box. This indicates finding the minimum value along the X-axis for the corresponding projection point within each child bounding box. This indicates finding the maximum value of the corresponding projection point along the Y-axis in each child bounding box. This indicates finding the minimum value along the Y-axis for the corresponding projection point within each child bounding box. This indicates finding the maximum value of the corresponding projection point along the Z-axis in each child bounding box. This indicates finding the minimum value of the corresponding projection point along the Z-axis in each child bounding box.

[0092] Calculate the voxel grid scale required for voxel downsampling at each level of point cloud The calculation formula is as follows: ; Where n is the current point cloud level, n = 2, 3, 4, ... The value of n corresponding to the root node is Its child node has an n value of " -1”, and so on; among them, the lowest level point cloud does not need to be sampled.

[0093] Next, coordinate quantization is performed, converting floating-point coordinates into 16- or 32-bit integers to reduce storage space. A specific example code is as follows: “def quantize_coordinates(points, range_min, range_max, quantization_bits=16): # Calculate scaling factor and offset scale = (range_max - range_min) / (2quantization_bits - 1) offset = range_min # Quantized coordinates quantized = ((points - offset) / scale).astype(np.uint16) return quantized, scale, offset".

[0094] Finally, following step 40, a 3DTiles tile set is constructed using the laser point cloud data. This tile set includes .pnts and .json files.

[0095] A common example of the generated directory structure is as follows: “output / 3dtiles / ├── tileset.json # Main configuration file ├── tiles / # Tile Data Directory │├── root.pnts# Root tiles │├── 1-0-0.pnts# First tile of the first layer │└── ... └── textures / # texture file (if it exists) Specifically, folders are created with the point cloud level number as the filename, and each folder stores the tile files of each level of point cloud. Each level of point cloud generates a corresponding pnts format point cloud tile file based on the divided voxel grid, and the point cloud tile file is named with the current voxel number. A text file named "tileset.json" is created in each level of point cloud folder to store the transformation matrix of each voxel relative to the center of the point cloud.

[0096] Among them, pnts is a file type in the 3D Tiles data format specifically used for point cloud data. Point cloud data is a dataset composed of a large number of small points, which is usually used to represent the three-dimensional shape and features of the real world. Each point can contain a variety of information, such as position, color, intensity, transparency, normal vector, temperature, reflectivity, etc. The pnts format PNTS file consists of two parts: (1) Header: contains metadata such as file type ("pnts"), version, and byte length; (2) Feature Table: stores information such as the number of points, bounding boxes, and coordinates. In one embodiment, it may also include a batch table for storing additional attributes (such as classification, building identifiers).

[0097] In the 3DTiles data format, JSON files are typically used as Tileset JSON files, which describe the structure, hierarchy, and attributes of the entire 3DTiles dataset. Its main features include: dataset structure description; Levels of Detail (LOD) information; spatial partitioning and indexing; metadata and attributes. Tiles are organized in a tree structure using a Hierarchical Level of Detail (HLOD) model to optimize spatial data rendering. Each tile has a bounding volume property and a content object defining the bounding scope of its control. The index tree has spatial continuity, with child tiles completely enclosed within the outer box of their parent tiles. The type and specific content of the generated 3DTiles tile set file are determined by those skilled in the art based on the specific use case and are not limited here.

[0098] Subsequently, the tile set obtained in step 40 is integrated into a 3D visualization platform (such as Cesium, Unity, and Unreal), and the point cloud tiles are dynamically loaded through the Application Programming Interface (API) for visualization verification.

[0099] Below is a specific visualization example: "01A Province Key Geological Hazard Area LiDAR Data Acquisition and Processing (2020) Results - Area B", the original point cloud data is 1.47TB (of which, there are 6000 files), and the generated 3D Tiles data size is 845.81GB; "02A Province Key Geological Hazard Area LiDAR Data Acquisition and Processing (2020) Results - Area B", the original point cloud data is 1.26TB (of which, there are 6177 files), and the generated 3D Tiles data size is 720.29GB; "A Province Key Geological Hazard Area LiDAR Data Acquisition and Processing (2020) Results - Area C", the original point cloud data is 1.13TB (of which, there are 8436 files), and the generated 3D Tiles data size is 684GB.

[0100] Figure 12 and Figure 13 The image shown is a preview of 3DTiles tile data generated using the method for converting massive laser point clouds to 3DTiles according to an embodiment of the present invention, after being published as a 3DTiles data service and loaded onto a 3D map. Figure 12 As can be seen, the method for converting massive laser point cloud LAS data into 3D tiles using the embodiments of the present invention can support the generation of 3D tiles from massive laser point cloud data; such as Figure 13 The visualization verification of the generated 3D tiles data preview shows that the 3D tiles data generated using the method of the present invention can be displayed normally and the integrity of the data can be guaranteed.

[0101] The method for converting massive laser point cloud LAS to 3DTiles according to embodiments of the present invention is designed for end-to-end point cloud data processing. It has a built-in native and complete LAS and / or LAZ format parsing engine, which can deeply adapt to the header information, extended fields and LAZ compression format decoding of all versions of LAS 1.0-1.4. It also optimizes the memory reading logic for massive data parsing. In particular, it automatically adopts a phased slicing and parallel processing strategy for ultra-large LAS files. The core processing flow is block and parallel processing, efficient key point extraction, generation of ".pnts" format files and tile set construction.

[0102] The method for converting massive laser point cloud LAS to 3DTiles using embodiments of the present invention can expand the format adaptation range of existing technologies, support the processing of LAS data containing extended fields, and avoid the phenomenon that some data blocks cannot be sliced ​​normally due to field parsing misalignment, resulting in only partial slicing. The parsing of LAZ compression format also does not rely on third-party plugins, avoiding the failure of the entire slicing task due to plugin version incompatibility or failure to install in advance. For scenarios requiring the processing of massive amounts of data, a set of processing logic is provided, introducing optimization techniques such as memory mapping and data block preloading. This avoids interruption of the slicing process due to memory overflow or input / output read blocking during processing, and avoids forced termination by the system due to memory usage exceeding the threshold, or slicing process stagnation due to excessive input / output read / write pressure. Therefore, slicing task failures or partial data slicing timeouts are avoided.

[0103] Example 2 like Figure 14 The diagram shown is a schematic representation of the architecture of a device for converting massive laser point clouds (LAS) to 3DTiles according to an embodiment of the present invention. The device for converting massive laser point clouds (LAS) to 3DTiles according to this embodiment includes one or more processors 21 and a memory 22. Figure 14 Take a processor 21 as an example.

[0104] Processor 21 and memory 22 can be connected via a bus or other means. Figure 14 Taking the example of a connection between China and Israel via a bus.

[0105] The memory 22, as a non-volatile computer-readable storage medium, can be used to store non-volatile software programs and non-volatile computer-executable programs, such as the method for converting massive laser point clouds to 3DTiles in Embodiment 1. The processor 21 executes the method for converting massive laser point clouds to 3DTiles by running the non-volatile software programs and instructions stored in the memory 22.

[0106] Memory 22 may include high-speed random access memory, and may also include non-volatile memory, such as at least one disk storage device, flash memory device, or other non-volatile solid-state storage device. In some embodiments, memory 22 may optionally include memory remotely located relative to processor 21, which can be connected to processor 21 via a network. Examples of such networks include, but are not limited to, the Internet, intranets, local area networks, mobile communication networks, and combinations thereof.

[0107] The program instructions / modules are stored in the memory 22. When executed by one or more processors 21, they perform the method for converting massive laser point clouds LAS to 3DTiles as described in Embodiment 1 above. For example, they perform each step of the method for converting massive laser point clouds LAS to 3DTiles as described above.

[0108] It is worth noting that the information interaction and execution process between the modules and units in the above-mentioned device and system are based on the same concept as the processing method embodiment of the present invention. For details, please refer to the description in the method embodiment of the present invention, and will not be repeated here.

[0109] Those skilled in the art will understand that all or part of the steps in the various methods of the embodiments can be implemented by a program instructing related hardware. The program can be stored in a computer-readable storage medium, which may include: read-only memory (ROM), random access memory (RAM), magnetic disk or optical disk, etc.

[0110] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for converting massive laser point clouds (LAS) into 3DTiles, characterized in that, include: Determine the spatial extent and edge information of the point cloud based on the input point cloud file; Based on the point cloud edge information and the point cloud spatial range, the point cloud data of the point cloud file is spatially divided to obtain multi-level data; Calculate the point cloud reading and rendering threshold for each level in the multi-level data to generate a point cloud index file; A tile set is constructed based on the point cloud index file and the multi-level data.

2. The method of converting massive laser point cloud LAS to 3D Tiles according to claim 1, wherein, The method includes: Calculate the centroid of the point cloud data; perform mean-reduction processing on each point in the point cloud data according to the centroid to obtain the first data; Calculate the covariance matrix of the first data to obtain the principal direction of the point cloud; correct the principal direction of the point cloud so that it coincides with the coordinate axis direction to obtain the second data. Calculate the point cloud edge information of the second data.

3. The method of converting massive laser point cloud LAS to 3D Tiles according to claim 2, wherein, The point cloud edge information includes the minimum bounding box and the point cloud boundary point set; The method includes: Obtain one of the current points of the second data, determine multiple nearest neighbor points of the current point, and obtain a local surface point set; The micro-tangent plane of the local surface point set is obtained by fitting; the angle between the current point and the corresponding projection points of the multiple nearest neighbor points is calculated using the projection of the micro-tangent plane as the reference plane. When the maximum value of the included angles among the plurality of nearest neighbor points is greater than the angle threshold, the current point is determined as a boundary point; Traverse all points in the second data, iteratively taking each point as the current point to obtain multiple boundary points, and use the multiple boundary points to construct a point cloud boundary point set; The minimum bounding box of the second data is solved using the point cloud boundary point set.

4. The method of converting massive laser point cloud LAS to 3D Tiles according to claim 2, wherein, The method further includes: The point cloud file is parsed to obtain point cloud attributes; wherein, the point cloud attributes include point cloud bounding boxes; The global bounding box of the second data is calculated using the point cloud bounding box to obtain the spatial extent of the point cloud.

5. The method of converting massive laser point cloud LAS to 3D Tiles according to claim 1, wherein, The method includes: Using the spatial range of the point cloud as the root tile, calculate the center point and geometric error of the root tile; Obtain the second data corresponding to the point cloud edge information; obtain the target point cloud in the second data, and determine the depth parameter using the geometric error and the number of points contained in the target point cloud; Using the center point as the reference point for octree segmentation, and according to the depth parameters, the point cloud edge information is used to iteratively construct the octree to complete the spatial division along the outer contour of the point cloud and obtain multi-level data.

6. The method for converting massive laser point clouds into 3DTiles according to claim 5, characterized in that, The point cloud edge information includes the minimum bounding box and the point cloud boundary point set; The method includes: The minimum bounding box is uniformly divided into multiple sub-bounding boxes along each directional axis, and the spatial information of all current sub-bounding boxes is determined. The sub-bounding boxes are used as child nodes of the octree, and the point cloud within the coverage of the spatial information is divided into the corresponding sub-bounding boxes to obtain the third data. Project the point cloud boundary point set onto the plane of each directional axis respectively; Determine the projection points of the boundary points in the point cloud boundary point set corresponding to different planes, and determine the maximum and minimum values ​​of the projection points on each directional axis in each child bounding box to obtain the range of the boundary points in each direction. The child bounding boxes corresponding to the minimum value are reorganized from the range domains in each direction. The child bounding boxes corresponding to the minimum value are used as parent nodes to generate new child nodes and child bounding boxes for the next layer. Iteratively reorganize the child bounding boxes corresponding to the minimum value until the depth of the octree satisfies the depth parameter, thus obtaining the octree.

7. The method for converting massive laser point clouds into 3DTiles according to claim 1, characterized in that, The method further includes: A GPU-based KD-tree spatial index is constructed for the multi-level data; Using the KD tree spatial index, for each current point in the multi-level data, determine the covariance matrix of the neighborhood of each point; use the local covariance matrix to perform eigenvalue decomposition to obtain the degree of change of each directional axis, and sort the eigenvalues ​​in descending order; Initial key points are identified using the feature value ratio condition based on the aforementioned feature values; The specified feature values ​​are used as response intensity indicators, and the initial key points with the strongest response are retained to obtain the target key point set.

8. The method for converting massive laser point clouds into 3DTiles according to claim 7, characterized in that, The feature value ratio condition is as follows: the change in the multi-level data in the second principal direction is greater than the change in the multi-level data in the first principal direction relative to the change in the multi-level data in the first principal direction, and the change in the multi-level data in the third principal direction is greater than the change in the second principal direction relative to the change in the second principal direction.

9. A non-volatile computer storage medium, characterized in that, The computer storage medium stores computer-executable instructions, which are executed by one or more processors to perform the method for converting massive laser point clouds into 3DTiles as described in any one of claims 1-8.

10. A device for converting massive laser point clouds into 3DTiles using LAS, characterized in that, include: At least one processor; And a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the processor for performing the method for converting massive laser point cloud LAS to 3DTiles as described in any one of claims 1-8.

Images