Building point cloud hierarchical reconstruction method fusing multi-feature constraints

Abstract

The application discloses a building point cloud layered reconstruction method fusing multiple feature constraints and belongs to the technical field of multiple feature constraints. After original point cloud data is acquired, indoor point cloud data sets are obtained by preprocessing, building a point cloud down-sampling data set, posture correction and extraction after screening. Point cloud space indexes are constructed, the curvatures of each point of the indoor point cloud data sets are calculated and sorted, the point with the minimum curvature is taken as an initial seed point, the neighborhood feature extraction is completed after the normal vector consistency and height difference constraints, the clustering and segmentation of the full point cloud are carried out, the point cloud cluster class attributes are obtained, the ground and ceiling attributes of the point cloud are divided by combining the support, level and topology three-factor coupling. Finally, the accurate extraction of the non-regular floor surface is realized through the progressive process of floor interval division and bidirectional cloth simulation. The application significantly improves the integrity and robustness of the layered reconstruction of the non-Manhattan building, and adapts to the layered needs of complex scenes such as curved facades and inclined structures.

Description

Technical Field

[0001] This invention belongs to the field of multi-feature constraint technology, specifically relating to a method for hierarchical reconstruction of building point clouds that integrates multiple feature constraints. Background Technology

[0002] As a core supporting technology in fields such as 3D building reconstruction, spatial information modeling, and smart building management, the building layering method plays an irreplaceable and important role in many key scenarios such as urban planning and design, indoor navigation and positioning, building operation and maintenance, and disaster emergency response.

[0003] In the field of urban planning, accurate building layering results can provide data support for floor area ratio calculation and spatial layout optimization, helping to scientifically formulate planning schemes. In indoor navigation scenarios, layering information can realize inter-floor path planning and spatial positioning, improving the accuracy and practicality of navigation services. In the building operation and maintenance phase, layering data can help managers quickly locate equipment distribution and monitor structural damage, reducing operation and maintenance costs and improving safety assurance levels. In disaster emergency response, a clear building layering structure can provide key spatial references for rescue route planning and search and rescue of trapped personnel, minimizing disaster losses. In addition, this method also provides basic data for emerging technologies such as digital twin buildings and BIM (Building Information Modeling) model construction, promoting the digital and intelligent transformation of the construction industry.

[0004] However, most existing mainstream building stratification methods are designed for Manhattan-style buildings with regular structures, perpendicular walls to the ground, and horizontal floor boundaries. Their core approach relies heavily on the regularity of building height distribution to achieve stratification, such as directly dividing different floor ranges by extracting peak and flat intervals from height distribution maps. The essence of these methods is to utilize the regular structural characteristics of Manhattan architecture—"horizontal stratification and vertical support"—to simplify the stratification logic and calculation process. In practical applications, many non-Manhattan buildings, such as those with curved facades, sloping structures, staggered floors, and irregular outlines, lack regular height distribution patterns. Their floor boundaries are not strictly horizontal, and their supporting structures exhibit non-vertical characteristics, making it difficult for traditional height distribution map-based stratification methods to accurately capture their true floor division logic. Summary of the Invention

[0005] Purpose of the invention: To provide a hierarchical reconstruction method for building point clouds that integrates multiple feature constraints, thereby solving the aforementioned problems existing in the prior art.

[0006] Technical solution: A method for hierarchical reconstruction of building point clouds integrating multiple feature constraints, comprising the following steps:

[0007] The process involves acquiring raw point cloud data, preprocessing the raw point cloud data to obtain point cloud downsampling, building a point cloud downsampling dataset, obtaining the point cloud downsampling dataset, performing attitude correction to align the building facade direction with the coordinate plane, obtaining an attitude-corrected point cloud dataset, and then filtering and extracting from the attitude-corrected point cloud dataset to obtain an indoor point cloud dataset.

[0008] Building point clouds Spatial indexing calculates the initial normal vector for each point in the indoor point cloud dataset, sets a preset angle threshold, and calculates the angle between the initial normal vector and the vertical direction. The non-vertical point cloud regions are filtered out. The curvature of each point in the indoor point cloud dataset is calculated, and the points are sorted according to their curvature. The point with the smallest curvature is used as the initial seed point. For the neighboring points of the initial seed point, the angle between their normal vectors is calculated sequentially. and Axis height difference When the angle between the normal vectors Within the included angle threshold range, the target point is marked as a member of the current cluster and added to the expansion queue to complete the neighborhood feature extraction;

[0009] Then, using the newly added target point as the initial seed point, repeat the above steps until no new points are added to the expanded queue, thus completing the extraction of the current cluster. Then, select the next unlabeled point as the new initial seed point according to the ascending curvature and repeat the above steps to complete the clustering and segmentation of the entire point cloud.

[0010] After completing the clustering and segmentation of the entire point cloud, cluster attributes are then distinguished to obtain point cloud cluster attributes. Elevation normalization is performed based on these attributes to obtain the elevation distribution characteristics of ground and ceiling point clouds. A preset elevation merging threshold is used to merge consecutive elevation intervals with adjacent elevation differences less than the threshold, resulting in several elevation intervals. The original cluster point cloud within each elevation interval is extracted and differential preprocessing is performed to obtain the ground elevation corresponding to the ground point cloud and the ceiling elevation corresponding to the ceiling point cloud. Further preprocessing of the ground and ceiling elevations yields the ground and ceiling point clouds, respectively. A structured trend surface model is then constructed, and the ground and ceiling trend surfaces are extracted based on the point cloud data. Spatial indexing matches the ground trend surface and the ceiling trend surface to obtain the upper and lower limits of the floor elevation. Point clouds whose elevations fall within this range are selected as mezzanine point clouds, thus achieving complete extraction of the target plane.

[0011] The preferred method for constructing a point cloud downsampling dataset is as follows:

[0012] First, the original point cloud data is acquired. Voxel downsampling is then used to divide the original data into several downsampled voxel grids. These downsampled voxel grids are then filtered to obtain non-dropout voxel grids. The centroids of all points within the non-dropout voxel grids are calculated, and these centroids are used to represent the corresponding voxels, thus obtaining the point cloud downsampling. The formula for calculating the centroids of all points within the non-dropout voxel grids is as follows:

[0013] In the formula: Indicates the center of mass of coordinate, Indicates the center of mass of coordinate, Indicates the center of mass of coordinate, This indicates the number of point clouds within the corresponding voxel. Indicates the first voxel in the downsampled voxel grid A number of points are used to calculate the centroid of all point clouds within the non-airborne sampling voxel grid. Indicates the first Points coordinate, Indicates the first Points coordinate, Indicates the first Points coordinate.

[0014] The preferred method for building the pose correction point cloud dataset is as follows:

[0015] pass The algorithm fits a plane, sets a preset threshold range, extracts point clouds of the horizontal plane with deviations within the threshold range, and calculates the horizontal normal vector. and Axis normal vector The angle between The calculation formula is as follows:

[0016] ;

[0017] After completing the horizontal normal vector vector The angle between After calculation, a rotation matrix is ​​constructed based on the Rodriguez rotation formula to complete the point cloud. Orientation and attitude correction, wherein the formula for calculating the rotation matrix is ​​as follows:

[0018] In the formula: This indicates the point cloud before rotation. This represents the point cloud after rotation. , , Let be the three components of the unit vector representing the axis of rotation; , , The formulas for calculating the three components are as follows: In the formula: Represents the unit vector of the rotation axis;

[0019] After completing the point cloud After axial correction, the image is projected onto the horizontal plane and rotated through the horizontal rotation angles in 1° increments. The optimal rotation angle is calculated using the minimum bounding rectangle area of ​​the projected point cloud as the loss function. After calculating the optimal rotation angle, the horizontal orientation of the point cloud is corrected using the formula for calculating the rotation matrix. The formula for calculating the optimal rotation angle is as follows:

[0020] ; In the formula: Indicates the length of the enclosing rectangle. This represents the width of the enclosing rectangle. This indicates the optimal rotation angle.

[0021] The preferred method for building an indoor point cloud dataset is as follows:

[0022] The process involves acquiring a 3D point cloud after attitude correction, projecting it onto a horizontal plane, setting a preset resolution, dividing the area into several grids, counting the number of 3D points in each grid, setting a preset threshold, and determining whether the threshold is met. If the number of points in each grid cell is greater than the threshold, the cell is assigned a value of 1; if the number of points in each grid cell is less than the threshold, the cell is assigned a value of 0. This generates a binary image of the 3D point cloud distribution. The binary image is then dilated to repair the boundary gaps. The longest closed contour in the binary image is extracted as the horizontal boundary of the building, thus enabling the selection and extraction of the attitude-corrected point cloud dataset and obtaining an indoor point cloud dataset.

[0023] Preferably, the neighborhood feature extraction process is as follows:

[0024] Obtain indoor point clouds from the indoor point cloud dataset as target points. , with target point Construct a spherical neighborhood with a radius equal to twice the size of a voxel, and build a set of neighborhood points. Through neighborhood point set To reflect the local geometric features of the target point, weighted principal component analysis is used to analyze the neighborhood point set. Perform fitting to solve for the target point. Construct the weighted covariance matrix of the optimal tangent plane at the point. First, the neighborhood points are calculated using the following formula. The weight is calculated using the following formula:

[0025] In the formula: Represents the target point, whose value is related to its neighboring points. To the target point European distance Inversely proportional; Indicates the downsampling voxel grid size; after completing the neighborhood points After calculating the weights, substitute them into the following formula to calculate the weighted covariance matrix. :

[0026] In the formula: For all neighboring points The sum of the weights is used to normalize the weights; then the weighted covariance matrix is... Eigenvalue decomposition is performed to obtain decomposed eigenvalues, and the eigenvector corresponding to the smallest eigenvalue is obtained as a point. The normal vector is used to calculate the curvature of a single point in the point cloud based on the decomposed eigenvalues. The calculation formula is as follows:

[0027] In the formula: , , Let them represent the covariance matrix respectively. The decomposed eigenvalues ​​obtained by eigenvalue decomposition are, where, , , Sort from smallest to largest; calculate the curvature and normal vector of a single point in the point cloud in turn, and complete the curvature sorting. The flatness of the location of a single point in the point cloud is determined by the magnitude of the curvature. The smaller the curvature, the flatter the location of the single point in the point cloud. Conversely, the larger the curvature, the more uneven the location of the single point in the point cloud.

[0028] A preset angle threshold is used to calculate the angle between the vector and the vertical direction based on the initial normal vector. The non-vertical point cloud regions are selected, and the curvature of each point in the indoor point cloud dataset is calculated. The points are then sorted according to the magnitude of curvature, and the point with the smallest curvature is selected as the initial seed point. For the neighboring points of the initial seed point, the angle between their normal vectors is calculated sequentially. and Axis height difference When the angle between the normal vectors Within the included angle threshold range, the target point is marked as a member of the current cluster and added to the extended queue to complete the neighborhood feature extraction.

[0029] Preferably, the steps for distinguishing point cloud cluster attributes are as follows:

[0030] For each target cluster, an outer convex envelope is extracted. Vertical point clouds within the outer convex envelope whose elevation difference satisfies a 1m threshold are selected as valid reference samples. Distance intervals are divided with a step size of 0.1m, and Gaussian weights are introduced to calculate a weighted score for each distance interval. The calculation formula is as follows:

[0031] In the formula: Indicates the first The distance from the interval point cloud to the fitting plane; This indicates the number of point clouds within that distance interval; The width of the Gaussian kernel is used to control the weight decay rate; then, the weighted scores of the point clouds above and below the plane are calculated separately. , Obtain the supporting confidence score Complete the analysis of the supporting structure, including the support confidence score. The calculation formula is as follows:

[0032] ;

[0033] When the support confidence score The closer the score is to 1, the more likely it is to be the point cloud above, and the higher the support confidence score. The closer the value is to 0, the more likely it is to be a lower point cloud. The cluster is then mapped to a downsampled voxel grid. Rays are emitted upwards and downwards from the starting voxel, with a preset height threshold used to filter out small-range elevation fluctuations. The attributes of the terminating voxel are then obtained. If the starting and terminating voxel attributes are the same, the starting voxel attributes are modified to satisfy vertical exclusivity. Finally, the number of ceiling point clouds within the cluster is counted. Total number of point clouds Calculate the confidence score of the hierarchical structure. Complete the hierarchical structure analysis, including calculating the hierarchical structure confidence score. The calculation formula is as follows:

[0034] ;

[0035] By comparing the consistency of attributes between the target cluster and neighboring clusters, the analysis of the vertical support structure is more consistent with the attributes of its neighbors, thus increasing... The weighting; when the vertical hierarchy structure validation is more consistent with the neighborhood attributes, the weighting is increased. The weights are calculated by... To complete the topology analysis, the calculation formula is as follows:

[0036] In the formula: This represents the scores for evaluating ceiling and floor attributes. The closer the score is to 1, the higher the ceiling point cloud score is, and the closer the score is to 0, the higher the floor point cloud score is. This represents the support attribute and the adjustment coefficient for the upper and lower cluster structures; the default value is 0.5. This represents the horizontal topology weight adjustment coefficient; it takes a positive value when the support structure analysis and the horizontal topology consistency result are consistent, and a negative value when the hierarchical structure verification and the horizontal topology consistency result are consistent. This is achieved through dynamic adjustment. This ensures that the final decision satisfies the triple constraints of local geometry, vertical structure, and horizontal topology.

[0037] Preferably, the elevation merging threshold setting process is as follows:

[0038] Using a cloth-simulated mesh as the resolution, the 3D point cloud is divided according to the resolution. In planar space, based on the maximum and minimum values ​​of the imported 3D point cloud. The horizontal boundary of the corresponding 3D point cloud is obtained. Based on the horizontal boundary, the number and spatial distribution of the grids in each row and column are calculated. The initial plane of the virtual cloth is set at a fixed distance above the highest point of the target point cloud to ensure that the grid covers the entire point cloud area.

[0039] For each grid after division, the maximum elevation value of the point cloud within the grid is extracted, and the maximum elevation value is used as the grid particle of the corresponding grid. The minimum elevation of the grid particle is used as the constraint threshold. Iterate through each movable grid particle and calculate the elevation of the grid particle under the influence of gravity using the following formula:

[0040] In the formula: Indicates the time step. Indicates constant acceleration. This represents the descent distance threshold within a time step. Indicates the mesh particles in the current iteration Elevation value at time, Indicates grid particles in Elevation value at time, Indicates grid particles in Elevation value at time, This represents the instantaneous vertical velocity of the grid particle; after obtaining the elevation value of the grid particle, the elevation value of the grid particle is updated to obtain the grid particle elevation value, and compared with the minimum elevation constraint threshold. The comparison is performed when the grid particle elevation value is less than or equal to the minimum elevation constraint threshold. Then the elevation value of the grid particle will be reset to the minimum elevation constraint threshold. And mark it as immovable;

[0041] For each movable mesh particle, extract the immovable mesh particles within its neighborhood. Based on the stiffness parameter and elevation difference, calculate the displacement adjustment under internal forces. The calculation formula is as follows:

[0042] In the formula: This indicates the amount of displacement adjustment under the action of internal forces. This represents the elevation difference between movable and immovable grid particles. For fabric rigidity parameters;

[0043] When the elevation of a movable grid particle is higher than that of a neighboring immovable grid particle, the corresponding displacement is adjusted along the negative elevation direction; otherwise, it is adjusted along the positive elevation direction.

[0044] Sort all movable grid particles by the distance to the nearest immovable grid particle. Traverse the sorted grid particles and their neighboring movable particles. Adjust the elevation of adjacent movable particles by bidirectional displacement to make them more consistent. This is the elevation merging threshold.

[0045] Preferably, the process of building a structured trend surface model is as follows:

[0046] Traverse the neighborhood of the immovable mesh particle, and when the lowest elevation constraint threshold of the movable particle in the neighborhood is reached, [the following is a process to determine the minimum elevation constraint threshold]. If the difference is less than the elevation merging threshold, the particle is marked as immovable, and its elevation is updated to its lowest elevation constraint threshold. They are then added to an expanded queue, and the iteration is repeated until no new particles are added to the queue.

[0047] Calculate the elevation difference between each original point cloud and its corresponding grid particle, and determine whether the elevation difference meets the elevation merging threshold. If it does, it is retained as the target plane and the target plane point set is constructed. Otherwise, the original point cloud is removed.

[0048] Traverse the neighborhood of each immovable mesh particle and calculate its minimum elevation constraint threshold relative to other immovable particles in the neighborhood. Difference: If the difference exceeds the elevation merging threshold, check if there exists an elevation in the lower grid that matches the minimum elevation constraint threshold of the higher grid. If point clouds with elevation differences within the elevation merging threshold exist, these supplementary point clouds are merged into the target plane point set. Consecutive elevation intervals in the target plane point set whose adjacent elevation differences are less than the elevation merging threshold are merged to obtain several elevation intervals. The original cluster point clouds within each elevation interval are extracted and differential preprocessing is performed to obtain the ground elevation corresponding to the ground-type point cloud and the ceiling elevation corresponding to the ceiling-type point cloud. The ground elevation and ceiling elevation are then preprocessed to obtain the ground point cloud and ceiling point cloud, respectively. The ground point cloud and ceiling point cloud are then extracted separately. , The coordinate range is generated according to the predetermined grid spacing. , Planar meshes, with mesh nodes calculated via linear interpolation. The initial trend surface is constructed using the value, and then the wavefront algorithm is used to fill in the missing values ​​and missing areas in the trend surface. At the same time, the average elevation of each type of initial trend surface is recorded, and a trend surface is constructed that includes at least grid coordinates. A structured trend surface model of the original point cloud values.

[0049] Preferably, the complete extraction process of the target plane is as follows:

[0050] Using pairwise pairings of ground trend surfaces and ceiling trend surfaces as spatial indices, the point cloud between the ground trend surface and the ceiling trend surface represents the point cloud of the target floor. The point cloud for each pair of trend surfaces is then calculated. , The coordinate horizontal intersection is used to filter out the core point cloud points that fall within this intersection range and mark them as unextracted point clouds. Inverse interpolation is then used to obtain the coordinates of these unextracted point clouds on the corresponding ceiling and ground trend surfaces. The value is used to obtain the upper and lower limits of the floor elevation. Point clouds whose elevations fall within this range are selected as mezzanine point clouds, and their masks are marked as extracted point clouds, thus achieving complete extraction of the target plane.

[0051] Beneficial Effects: This invention relates to a hierarchical reconstruction method for building point clouds that integrates multiple feature constraints. First, the original point cloud data is divided and its attitude corrected to align the building facade direction with the coordinate plane, resulting in an attitude-corrected point cloud dataset. Then, the attitude-corrected point cloud dataset is filtered and extracted to obtain an indoor point cloud dataset. Next, features are extracted from the points in the indoor point cloud dataset to complete the clustering and segmentation of the entire point cloud. Finally, a method is constructed... Spatial indexing, combined with neighborhood point weight calculation, normal vector and curvature feature extraction, and grid particle iterative growth strategy, avoids the undersegmentation problem caused by Euclidean clustering's reliance on distance thresholds, and solves the oversegmentation phenomenon caused by the sensitivity of region growth clustering to normal vector fluctuations, thus achieving complete and accurate division of complex structural areas in non-Manhattan buildings.

[0052] Secondly, a three-factor weighted judgment model was constructed, which includes analysis of the supporting structure, analysis of the hierarchical structure, and analysis of the topology. By quantifying the distribution characteristics of the vertical point cloud, verifying the mutual exclusivity of the upper and lower voxel attributes, and dynamically adjusting the feature weight ratio, the interference of auxiliary clusters such as ground protrusions and ceiling droops was effectively suppressed. This solved the problem of attribute misjudgment that easily occurs in the near-ground area when the single-feature or dual-feature model has incomplete viewpoints, and achieved accurate identification of ground and ceiling clusters.

[0053] Finally, a planar extraction scheme combining floor interval division and trend surface was adopted, by setting a minimum elevation constraint threshold for grid particles. By introducing movable particle sorting and cross-grid point cloud retrieval post-processing, it not only avoids the traditional The algorithm (Cloth Simulation Filter) is prone to the problem of excessive sinking of the cloth surface. It also solves the problems of poor plane fitting for discontinuous elevations and fragmentation of the fitting plane in non-vertical structure scenarios. It achieves complete extraction of the target plane, significantly improves the accuracy, completeness and robustness of the layered reconstruction of non-Manhattan buildings, and effectively adapts to the layering requirements of complex building scenarios such as curved facades, inclined structures and staggered floors. Attached Figure Description

[0054] Figure 1 This is a system block diagram of the present invention;

[0055] Figure 2 This is a flowchart of the floor division process of the present invention;

[0056] Figure 3 This is a schematic diagram of point cloud attitude correction according to the present invention;

[0057] Figure 4 This is a flowchart of the point cloud clustering and segmentation process of the present invention;

[0058] Figure 5 This is a comparison diagram of clustering and segmentation in this invention;

[0059] Figure 6 This is a schematic diagram illustrating the three-factor point cloud attribute determination method of the present invention.

[0060] Figure 7 This is a diagram showing the ablation experiment results for the three-factor attribute determination of this invention;

[0061] Figure 8 This is a schematic diagram illustrating the influence range of the internal forces in this invention;

[0062] Figure 9 This is a schematic diagram for verifying the rationality of the present invention;

[0063] Figure 10 For the purposes of this invention Architectural planar filtering effect;

[0064] Figure 11 This is the floor division result of the sample data in this invention;

[0065] Figure 12 This is a diagram showing the floor division results using the box plot method of the present invention;

[0066] Figure 13 This is a diagram showing the floor division results of the bidirectional hierarchical separation method of the present invention;

[0067] Figure 14 This is a schematic diagram of clustering segmentation in this invention. Detailed Implementation

[0068] like Figures 1 to 13 As shown, the present invention provides a technical solution: a method for hierarchical reconstruction of building point clouds by integrating multiple feature constraints, comprising the following steps:

[0069] The process involves acquiring raw point cloud data, preprocessing it to obtain downsampled point cloud data, and constructing a downsampled point cloud dataset. The construction process is as follows: First, acquire the raw point cloud data. Then, use voxel downsampling to divide the raw data into several downsampled voxel grids. Next, sieve these downsampled voxel grids to obtain non-dropout voxel grids. Finally, calculate the centroid of all points within the non-dropout voxel grid and use this centroid to represent the corresponding voxel, thus obtaining the downsampled point cloud. The formula for calculating the centroid of all points within the non-dropout voxel grid is as follows:

[0070] In the formula: Indicates the center of mass of coordinate; Indicates the center of mass of coordinate; Indicates the center of mass of coordinate; Indicates the number of point clouds within the corresponding voxel; Indicates the first voxel in the downsampled voxel grid A number of points are used to calculate the centroid of all point clouds within the non-airborne sampling voxel grid; Indicates the first Points coordinate; Indicates the first Points coordinate; Indicates the first Points coordinate.

[0071] A point cloud downsampling dataset is obtained, and attitude correction is performed to align the building facade orientation with the coordinate plane, resulting in an attitude-corrected point cloud dataset. The process of building the attitude-corrected point cloud dataset is as follows: [Further details omitted]. The algorithm fits a plane, sets a preset threshold range, extracts point clouds of the horizontal plane with deviations within the threshold range, and calculates the horizontal normal vector. and Axis normal vector The angle between The calculation formula is as follows: ; after completing the horizontal normal vector and Axis normal vector The angle between After calculation, a rotation matrix is ​​constructed based on the Rodriguez rotation formula to complete the point cloud. Axial orientation correction, where the rotation matrix is ​​calculated using the following formula: In the formula: This indicates the point cloud before rotation. This represents the point cloud after rotation. , , Let be the three components of the unit vector representing the axis of rotation; , , The formulas for calculating the three components are as follows: In the formula: Represents the unit vector of the rotation axis;

[0072] After completing the point cloud After axial correction, the image is projected onto the horizontal plane and rotated through the horizontal rotation angles in 1° increments. The optimal rotation angle is calculated using the minimum bounding rectangle area of ​​the projected point cloud as the loss function. After calculating the optimal rotation angle, the horizontal orientation of the point cloud is corrected using the formula for calculating the rotation matrix. The formula for calculating the optimal rotation angle is as follows: ; In the formula: This indicates the length of the enclosing rectangle. This represents the width of the enclosing rectangle. The optimal rotation angle is represented. The original point cloud data is divided and its attitude is corrected to align the building facade with the coordinate plane. The attitude-corrected 3D point cloud is obtained and projected onto the horizontal plane. A preset resolution is used to divide the data into several grids. The number of 3D point clouds in each grid is counted. A preset threshold is used to determine whether the threshold requirement is met. If the number of point clouds in each grid is greater than the threshold, the grid is assigned a value of 1. If the number of point clouds in each grid is less than the threshold, the grid is assigned a value of 0. A binary image of the 3D point cloud distribution is generated. The binary image is dilated to repair the boundary gaps. The longest closed contour in the binary image is extracted as the horizontal boundary of the building to achieve the filtering and extraction of the attitude-corrected point cloud dataset, resulting in an indoor point cloud dataset.

[0073] Building point clouds Spatial indexing involves calculating the normal vector for each point in the indoor point cloud dataset to obtain an initial normal vector. A preset angle threshold is used, and the curvature of each point is calculated and sorted. The point with the smallest curvature is selected as the initial seed point. Normal vector consistency and height difference constraints are applied to complete neighborhood feature extraction. The neighborhood feature extraction process is as follows: Figure 14 As shown, indoor point clouds from the indoor point cloud dataset are obtained as target points. , with target point Construct a spherical neighborhood with a radius equal to twice the size of a voxel, and build a set of neighborhood points. ,like Figure 14 As shown, the neighborhood point set That is Figure 14 All green points inside the sphere, through the set of neighboring points To reflect the local geometric features of the target point, weighted principal component analysis is used to analyze the neighborhood point set. Perform fitting to solve for the target point. Construct the weighted covariance matrix of the optimal tangent plane at the point. First, the neighborhood points are calculated using the following formula. The weight is calculated using the following formula: In the formula: Represents the target point, whose value is related to its neighboring points. To the target point European distance Inversely proportional; Indicates the downsampling voxel grid size; after completing the neighborhood points After calculating the weights, substitute them into the following formula to calculate the weighted covariance matrix. : In the formula: For all neighboring points The sum of the weights is used to normalize the weights; then the weighted covariance matrix is... Eigenvalue decomposition is performed to obtain decomposed eigenvalues, and the eigenvector corresponding to the smallest eigenvalue is taken as the target point. The normal vector is used to calculate the curvature of a single point in the point cloud based on the decomposed eigenvalues. The calculation formula is as follows: In the formula: , , Let them represent the covariance matrix respectively. The decomposed eigenvalues ​​obtained by eigenvalue decomposition are, where, , , Sort the points in ascending order; calculate the curvature and normal vector of each point in the point cloud sequentially, and sort by curvature. The flatness of the point's location is determined by the magnitude of the curvature; the smaller the curvature, the flatter the location of the point; conversely, the larger the curvature, the more uneven the location. A preset angle threshold is used to calculate the angle between the initial normal vector and the vertical direction. The non-vertical point cloud region is selected, which is the point cloud after removing the walls. The curvature of each point in the indoor point cloud dataset is calculated and sorted according to the curvature. The point with the smallest curvature is selected as the initial seed point. For the neighboring points of the initial seed point, the angle between the normal vector of each point in the neighborhood and the initial seed point is calculated. and Axis height difference When the angle between the normal vectors Within the included angle threshold range, the target point is marked as a member of the current cluster and added to the expansion queue to complete neighborhood feature extraction. Then, the newly added target point is used as the initial seed point, and the above steps are repeated until no new points are added to the expansion queue, completing the extraction of the current cluster. Next, the next unmarked point is selected in ascending curvature order as the new initial seed point, and the above steps are repeated to complete the clustering and segmentation of the entire point cloud. In this embodiment, the clustering and segmentation are as follows: Figure 14 As shown, the red points remain unchanged, and the green points have had their normal vectors determined. Axis height difference Then, a new green point is added. After all green points have been judged, the newly added green point is used as a red point for further judgment. This process continues until no point is satisfied. Then, a new point with a smaller curvature is selected as a new red point, and the judgment process is repeated to complete the clustering and segmentation.

[0074] After completing the clustering and segmentation of the entire point cloud, cluster attribute differentiation is performed to obtain the point cloud cluster attributes. The steps for differentiating point cloud cluster attributes are as follows: extract the convex envelope of the target cluster, select vertical point clouds within the convex envelope whose elevation difference meets the 1m threshold as valid reference samples, divide the distance interval with a step size of 0.1m, introduce Gaussian weights, and calculate the weighted score of the distance interval. The calculation formula is as follows: In the formula: Indicates the first The distance from the interval point cloud to the fitted plane; This indicates the number of point clouds within that distance interval; The width of the Gaussian kernel is used to control the weight decay rate; then, the weighted scores of the point clouds above and below the plane are calculated separately. , Obtain the support confidence score Complete the analysis of the supporting structure, including the support confidence score. The calculation formula is as follows: When the support confidence score The closer the score is to 1, the more likely it is to be the point cloud above, and the higher the support confidence score. The closer the value is to 0, the more likely it is to be a lower point cloud. The cluster is then mapped to a downsampled voxel grid. Rays are emitted upwards and downwards from the starting voxel, with a preset height threshold used to filter out small-range elevation fluctuations. The attributes of the terminating voxel are then obtained. If the starting and terminating voxel attributes are the same, the starting voxel attributes are modified to satisfy vertical exclusivity. Finally, the number of ceiling point clouds within the cluster is counted. Total number of point clouds Calculate the confidence score of the hierarchical structure. Complete the hierarchical structure analysis, including calculating the hierarchical structure confidence score. The calculation formula is as follows: By comparing the consistency of attributes between the target cluster and neighboring clusters, if the vertical support structure analysis is more consistent with the neighboring attributes, then the calculation is increased. The weights are determined by the vertical support structure analysis and the consistency judgment of neighborhood attributes. First, the neighborhood clusters are determined. The neighborhood clusters are based on the horizontal outline of the target cluster. The adjacent clusters with a spatial distance of ≤1 meter from the target cluster are selected as the neighborhood clusters. The support confidence scores of the target cluster and the neighborhood clusters are then extracted. And binarize it, where, when the support confidence score When the support confidence score is ≥0.5, it is considered the ceiling. If the value is less than 0.5, it is considered ground level. If the vertical support structure analysis of the target cluster and the neighboring clusters is completely identical, it is determined that the vertical support structure analysis is consistent with the neighboring attributes. When the vertical hierarchy structure verification is more consistent with the neighboring attributes, it is increased. The weights are determined by the vertical hierarchy structure verification and the neighborhood attribute consistency judgment. First, the neighborhood clusters are determined. The neighborhood clusters are based on the horizontal outline of the target cluster. The neighboring clusters with a spatial distance of ≤1 meter from the target cluster are selected as the neighborhood clusters. The hierarchy structure confidence scores between the target cluster and the neighborhood clusters are then extracted. And binarize it, where, when the hierarchical confidence score When the confidence score is ≥0.5, it is considered the ceiling. If the value is less than 0.5, it is determined to be ground level. If the vertical hierarchical structure verification of the target cluster and the neighboring clusters is exactly the same, it is determined that the vertical hierarchical structure verification is consistent with the neighborhood attributes; then, through calculation... To complete the topology analysis, the calculation formula is as follows: In the formula: This represents the scores for evaluating ceiling and floor attributes. The closer the score is to 1, the higher the ceiling point cloud score is, and the closer the score is to 0, the higher the floor point cloud score is. This represents the support attribute and the adjustment coefficient for the upper and lower cluster structures; the default value is 0.5. This represents the horizontal topology weight adjustment coefficient; it takes a positive value when the support structure analysis and the horizontal topology consistency result are consistent, and a negative value when the hierarchical structure verification and the horizontal topology consistency result are consistent. This is achieved through dynamic adjustment. This ensures that the final decision satisfies the triple constraints of local geometry, vertical structure, and horizontal topology.

[0075] In this embodiment, as Figure 5 , Figure 6 and Figure 7 As shown, while using only the supporting structure factor for single-factor attribute differentiation can quickly divide the main areas of the ground and ceiling, it can also lead to confusion in attribute identification for auxiliary clusters attached to the ground or ceiling surfaces, such as ground protrusions and ceiling drooping structures, because they are inversely related to the characteristics of the floor and ceiling. Therefore, after introducing a hierarchical structure to construct a two-factor fusion model, the interference of auxiliary clusters is effectively suppressed and attribute misjudgment is significantly reduced. However, in the near-ground area and in areas with incomplete viewpoints, the hierarchical structure features fail due to height threshold settings or viewpoint loss, and attribute misjudgment still exists. In contrast, introducing a topological structure to achieve three-factor joint analysis by modeling the global scene connectivity and spatial distribution pattern can effectively avoid the above limitations and achieve accurate division of the ground and ceiling.

[0076] Elevation normalization was performed based on point cloud cluster attributes to obtain the elevation distribution characteristics of ground and ceiling point clouds. A preset elevation merging threshold was established, and the threshold setting process is as follows: using the cloth simulation mesh as the resolution, the 3D point cloud was divided according to the resolution. In planar space, based on the maximum and minimum values ​​of the imported 3D point cloud. The horizontal boundary of the corresponding 3D point cloud is obtained, and a mesh is created according to the resolution. Each mesh represents the virtual cloth. By simulating a sinking process, the desired ceiling or floor trend surface can be obtained. Based on the horizontal boundary, the number and spatial distribution of the corresponding meshes in each row and column are calculated. The initial plane of the virtual cloth is set at a fixed distance above the highest point of the target point cloud to ensure that the mesh covers the entire point cloud area. For each mesh, the maximum elevation value of the point cloud inside the mesh is extracted, and the maximum elevation value is used as the mesh particle of the corresponding mesh. The minimum elevation constraint threshold of the mesh particle is used as the constraint threshold. Iterate through each movable grid particle and calculate the elevation of the grid particle under the influence of gravity using the following formula: In the formula: Indicates the time step. Indicates constant acceleration. This represents the descent distance threshold within a time step. Indicates the mesh particles in the current iteration Elevation value at time, Indicates grid particles in Elevation value at time, Indicates grid particles in Elevation value at time, This represents the instantaneous vertical velocity of the grid particle; after obtaining the elevation value of the grid particle, the elevation value of the grid particle is updated to obtain the grid particle elevation value, and compared with the minimum elevation constraint threshold. The comparison is performed when the grid particle elevation value is less than or equal to the minimum elevation constraint threshold. Then the elevation value of the grid particle will be reset to the minimum elevation constraint threshold. The immovable mesh particles are marked as immovable. For each movable mesh particle, the immovable mesh particles in its neighborhood are extracted. Based on the stiffness parameter and elevation difference, the displacement adjustment under the action of internal forces is calculated. The calculation formula is as follows: In the formula: This indicates the amount of displacement adjustment under the action of internal forces. This represents the elevation difference between movable and immovable grid particles. This is the fabric rigidity parameter; when the elevation of a movable grid particle is higher than that of its neighboring immovable grid particles, the corresponding displacement is adjusted along the negative elevation direction, and vice versa; all movable grid particles are sorted according to the distance of the nearest immovable grid particle, and the sorted grid particles and their neighboring movable particles are traversed. The elevations of adjacent movable particles are adjusted in both directions to make them consistent; this is the elevation merging threshold. For example... Figure 9 As shown, regarding the internal force adjustment mechanism of movable particles, the simulation results of circular hole data with and without particle sorting are compared. Figure 9 As shown in a and b, the results demonstrate that after introducing movable particle sorting, the particle distribution is more uniform and better reflects the true laws of physical settlement; secondly, the optimization effect of the post-processing step on cross-elevation fabric bonding is verified, such as... Figure 9 As shown in c and e, the fabric simulation results without post-processing show that when crossing two different elevations, the fabric surface is constrained by the internal forces of immovable particles at different elevations, making it impossible to achieve complete fit to the target plane. Figure 9 As shown in c, a stepped suspended particle layer was formed. The problem was effectively improved after introducing both movable and immovable particles in the post-processing, and the fabric surface basically achieved continuous bonding of high and low planes, as shown in... Figure 9 As shown in d, however, due to the single-value elevation characteristic of the initial grid, a single grid can only extract ground points with the lowest elevation, causing target points with higher elevations to be misclassified as non-ground points and lost. By adding point cloud retrieval post-processing between immovable particle grids, target points on both high and low planes are effectively preserved, such as... Figure 9 As shown in 'e', ​​this effectively improves the completeness and accuracy of the filter, such as... Figure 10As shown, the traditional CSF algorithm without adjusted time step, when filtering point clouds in architectural scenes, is easily constrained by the coexistence of horizontal and steeply descending regions in the scene, leading to excessive sinking of the cloth surface. This results in compromised integrity of the extracted target plane, making it difficult to obtain complete and high-precision architectural target planes. In both complete and non-vertical structure scenes, a large number of target plane points cannot be effectively extracted. Adjusting the time step parameter alleviates the excessive sinking of the cloth surface to some extent, enabling the effective identification of some unextracted target plane points. However, it still struggles to achieve ideal extraction results for inter-plane regions with discontinuous elevations. Furthermore, the continuity of the fitted plane in non-vertical structure scenes is not substantially improved, and the problem of plane fragmentation persists. In contrast, the improved CSF algorithm proposed in this paper, by setting a minimum elevation constraint threshold... By selecting and optimizing the internal force constraints, the problem of excessive sinking of the fabric surface can be effectively avoided, and the simulated fabric surface can be made to accurately fit the spatial morphological characteristics of the target building plane.

[0077] Consecutive elevation intervals with adjacent elevation differences less than the elevation merging threshold are merged to obtain several elevation intervals. The original cluster point cloud within each elevation interval is extracted and differential preprocessing is performed to obtain the ground elevation corresponding to the ground-type point cloud and the ceiling elevation corresponding to the ceiling-type point cloud. The ground elevation and ceiling elevation are then preprocessed to obtain the ground point cloud and ceiling point cloud, respectively. A structured trend surface model is then constructed. The process of constructing the structured trend surface model is as follows: traverse the neighborhood range of immovable grid particles. When the minimum elevation constraint threshold between the movable particles in the neighborhood and the current particle is reached... If the difference is less than the elevation merging threshold, the particle is marked as immovable, and its elevation is updated to its lowest elevation constraint threshold. The original point cloud is then added to the expansion queue, and the iteration is repeated until no new particles are added to the queue. The elevation difference between each original point cloud and its corresponding grid particle is calculated, and it is determined whether the elevation difference meets the elevation merging threshold. If it does, the original point cloud is retained as the target plane and the target plane point set is constructed. Otherwise, the original point cloud is removed.

[0078] Traverse the neighborhood of each immovable mesh particle and calculate its minimum elevation constraint threshold relative to other immovable particles in the neighborhood. Difference: If the difference exceeds the elevation merging threshold, check if there exists an elevation in the lower grid that matches the minimum elevation constraint threshold of the higher grid. If point clouds with elevation differences within the elevation merging threshold exist, these supplementary point clouds are merged into the target plane point set. Consecutive elevation intervals in the target plane point set whose adjacent elevation differences are less than the elevation merging threshold are merged to obtain several elevation intervals. The original cluster point clouds within each elevation interval are extracted and differential preprocessing is performed to obtain the ground elevation corresponding to the ground-type point cloud and the ceiling elevation corresponding to the ceiling-type point cloud. The ground elevation and ceiling elevation are then preprocessed to obtain the ground point cloud and ceiling point cloud, respectively. The ground point cloud and ceiling point cloud are then extracted separately. , The coordinate range is generated according to the predetermined grid spacing. , Planar meshes, with mesh nodes calculated via linear interpolation. The initial trend surface is constructed using values, and then the wavefront algorithm is used to fill in the missing values ​​and areas in the trend surface. The wavefront algorithm involves using existing values ​​from the trend surface... The grid nodes with values ​​are treated as "starting points / wave sources," and like the diffusion of water waves, starting from these grid nodes, the surrounding empty / missing grids are gradually filled. This method uses interpolation to complete the entire trend surface, recording the average elevation of each initial trend surface and constructing a system that includes at least grid coordinates and trend surfaces. A structured trend surface model of the original point cloud values.

[0079] Using pairwise pairings of ground trend surfaces and ceiling trend surfaces as spatial indices, the point cloud between the ground trend surface and the ceiling trend surface represents the point cloud of the target floor. The point cloud for each pair of trend surfaces is then calculated. , The coordinate horizontal intersection is used to filter out the core point cloud points that fall within this intersection range and mark them as unextracted point clouds. Inverse interpolation is then used to obtain the coordinates of these unextracted point clouds on the corresponding ceiling and ground trend surfaces. The value is used to obtain the upper and lower limits of the floor elevation. Point clouds whose elevations fall within this range are selected as mezzanine point clouds, and their masks are marked as "extracted point clouds", thus achieving complete extraction of the target plane.

[0080] like Figures 11 to 13 As shown, the above technical solution is used to verify the floor division in multiple sample scenarios. Multiple data samples are selected as verification scenarios. The floor division of the data samples is observed using the method in this paper. The verification results show that the method can achieve accurate floor division in different sample scenarios, with clear layer boundaries and stable results. It effectively avoids the problems of layer confusion and blurred boundaries that are prone to occur in complex scenarios by traditional methods, and fully verifies the effectiveness and robustness of the technical solution in this paper.

[0081] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A method for hierarchical reconstruction of building point clouds integrating multiple feature constraints, characterized in that, Includes the following steps: The process involves acquiring raw point cloud data, preprocessing the raw point cloud data to obtain point cloud downsampling, building a point cloud downsampling dataset, obtaining the point cloud downsampling dataset, performing attitude correction to align the building facade direction with the coordinate plane, obtaining an attitude-corrected point cloud dataset, and then filtering and extracting from the attitude-corrected point cloud dataset to obtain an indoor point cloud dataset. Building point clouds Spatial indexing calculates the normal vector for each point in the indoor point cloud dataset to obtain the initial normal vector. A preset angle threshold is set, the curvature of each point is calculated and sorted, and the point with the smallest curvature is used as the initial seed point. Normal vector consistency and height difference constraints are performed to complete the extraction of neighborhood features. Then, select the next unlabeled point as the new initial seed point according to the ascending curvature and repeat the above steps to complete the clustering and segmentation of the entire point cloud; After the clustering and segmentation of the entire point cloud is completed, the cluster attributes are distinguished. The cluster attributes of the point cloud are obtained by combining the three factors of support, hierarchy and topology. The elevation normalization is performed based on the cluster attributes of the point cloud to obtain the elevation distribution characteristics of the ground and ceiling point clouds respectively. The elevation merging threshold is preset. The ground and ceiling point clouds are preprocessed according to the elevation merging threshold to obtain the ground point cloud and ceiling point cloud, and a structured trend surface model is built. Based on the structured trend surface model, the ground trend surface and ceiling trend surface are extracted to achieve complete extraction of the target plane.

2. The method for hierarchical reconstruction of building point clouds integrating multiple feature constraints according to claim 1, characterized in that, The process of building the point cloud downsampling dataset is as follows: First, the original point cloud data is acquired. Voxel downsampling is then used to divide the original data into several downsampled voxel grids. These downsampled voxel grids are then filtered to obtain non-dropout voxel grids. The centroids of all points within the non-dropout voxel grids are calculated, and these centroids are used to represent the corresponding voxels, thus obtaining the point cloud downsampling. The formula for calculating the centroids of all points within the non-dropout voxel grids is as follows: In the formula: Indicates the center of mass of coordinate; Indicates the center of mass of coordinate; Indicates the center of mass of coordinate; Indicates the number of point clouds within the corresponding voxel; Indicates the first voxel in the downsampled voxel grid A number of points are used to calculate the centroid of all point clouds within the non-airborne sampling voxel grid; Indicates the first Points coordinate; Indicates the first Points coordinate; Indicates the first Points coordinate.

3. The method for hierarchical reconstruction of building point clouds integrating multiple feature constraints according to claim 1, characterized in that, The process of building the pose correction point cloud dataset is as follows: pass The algorithm fits a plane, sets a preset threshold range, extracts point clouds of the horizontal plane with deviations within the threshold range, and calculates the horizontal normal vector. and Axis normal vector The angle between The calculation formula is as follows: ； After completing the horizontal normal vector and Axis normal vector The angle between After calculation, a rotation matrix is ​​constructed based on the Rodriguez rotation formula to complete the point cloud. Axial orientation correction, where the rotation matrix is ​​calculated using the following formula: In the formula: This indicates the point cloud before rotation. This represents the point cloud after rotation. , , Let be the three components of the unit vector representing the axis of rotation; , , The formulas for calculating the three components are as follows: In the formula: Represents the unit vector of the rotation axis; After completing the point cloud After axial correction, the image is projected onto the horizontal plane and rotated through the horizontal rotation angles in 1° increments. The optimal rotation angle is calculated using the minimum bounding rectangle area of ​​the projected point cloud as the loss function. After calculating the optimal rotation angle, the horizontal orientation of the point cloud is corrected using the formula for calculating the rotation matrix. The formula for calculating the optimal rotation angle is as follows: ; In the formula: Indicates the length of the enclosing rectangle. This represents the width of the enclosing rectangle. This indicates the optimal rotation angle.

4. The method for hierarchical reconstruction of building point clouds incorporating multiple feature constraints according to claim 2, characterized in that, The process of building an indoor point cloud dataset is as follows: The process involves acquiring a 3D point cloud after attitude correction, projecting it onto a horizontal plane, setting a preset resolution, dividing the area into several grids, counting the number of 3D points in each grid, setting a preset threshold, and determining whether the threshold is met. If the number of points in each grid cell is greater than the threshold, the cell is assigned a value of 1; if the number of points in each grid cell is less than the threshold, the cell is assigned a value of 0. This generates a binary image of the 3D point cloud distribution. The binary image is then dilated to repair the boundary gaps. The longest closed contour in the binary image is extracted as the horizontal boundary of the building, thus enabling the selection and extraction of the attitude-corrected point cloud dataset and obtaining an indoor point cloud dataset.

5. The method for hierarchical reconstruction of building point clouds incorporating multiple feature constraints according to claim 1, characterized in that, The process of neighborhood feature extraction is as follows: Obtain indoor point clouds from the indoor point cloud dataset as target points. , with target point Construct a spherical neighborhood with a radius equal to twice the size of a voxel, and build a set of neighborhood points. Through neighborhood point set To reflect the local geometric features of the target point, weighted principal component analysis is used to analyze the neighborhood point set. Perform fitting to solve for the target point. Construct the weighted covariance matrix of the optimal tangent plane at the point. First, the neighborhood points are calculated using the following formula. The weight is calculated using the following formula: In the formula: Represents the target point, whose value is related to its neighboring points. To the target point European distance Inversely proportional; Indicates the downsampling voxel grid size; after completing the neighborhood points After calculating the weights, substitute them into the following formula to calculate the weighted covariance matrix. : In the formula: For all neighboring points The sum of the weights is used to normalize the weights; then the weighted covariance matrix is... Eigenvalue decomposition is performed to obtain decomposed eigenvalues, and the eigenvector corresponding to the smallest eigenvalue is obtained as the target point. The normal vector is used to calculate the curvature of a single point in the point cloud based on the decomposed eigenvalues. The calculation formula is as follows: In the formula: , , Let them represent the covariance matrix respectively. The decomposed eigenvalues ​​obtained by eigenvalue decomposition are, where, , , Sort the points in ascending order; calculate the curvature and normal vector of each point in the point cloud sequentially, and sort by curvature. The flatness of the point's location is determined by the magnitude of the curvature; the smaller the curvature, the flatter the location of the point; conversely, the larger the curvature, the more uneven the location. A preset angle threshold is used to calculate the angle between the initial normal vector and the vertical direction. The non-vertical point cloud regions are filtered out, the curvature of each point in the indoor point cloud dataset is calculated, and the points are sorted according to the curvature. The point with the smallest curvature is selected as the initial seed point. For the neighboring points of the initial seed point, the angle between the normal vector of each point in the neighborhood and the initial seed point is calculated. and Axis height difference When the angle between the normal vectors Within the included angle threshold range, the target point is marked as a member of the current cluster and added to the extended queue to complete the neighborhood feature extraction.

6. The method for hierarchical reconstruction of building point clouds incorporating multiple feature constraints according to claim 2, characterized in that, The steps for distinguishing point cloud cluster attributes are as follows: For each target cluster, an outer convex envelope is extracted. Vertical point clouds within the outer convex envelope whose elevation difference satisfies a 1m threshold are selected as valid reference samples. Distance intervals are divided with a step size of 0.1m, and Gaussian weights are introduced to calculate a weighted score for each distance interval. The calculation formula is as follows: In the formula: Indicates the first The distance from the interval point cloud to the fitting plane; This indicates the number of point clouds within that distance interval; This represents the Gaussian kernel width and controls the weight decay rate. Then calculate the weighted scores of the point clouds above and below the plane respectively. , Obtain the supporting confidence score Complete the analysis of the supporting structure, including the support confidence score. The calculation formula is as follows: ； When the support confidence score The closer the score is to 1, the more likely it is to be the point cloud above, and the higher the support confidence score. The closer the value is to 0, the more likely it is to be a lower point cloud. The cluster is then mapped to a downsampled voxel grid. Rays are emitted upwards and downwards from the starting voxel, with a preset height threshold used to filter out small-range elevation fluctuations. The attributes of the terminating voxel are then obtained. If the starting and terminating voxel attributes are the same, the starting voxel attributes are modified to satisfy vertical exclusivity. Finally, the number of ceiling point clouds within the cluster is counted. Total number of point clouds Calculate the confidence score of the hierarchical structure. Complete the hierarchical structure analysis, including calculating the hierarchical structure confidence score. The calculation formula is as follows: ； By comparing the consistency of attributes between the target cluster and neighboring clusters, if the vertical support structure analysis is consistent with the neighboring attributes, then... The weight is increased when the vertical hierarchy validation matches the neighborhood attributes. The weights are calculated by... To complete the topology analysis, the calculation formula is as follows: In the formula: This represents the scores for evaluating ceiling and floor attributes. The closer the score is to 1, the higher the ceiling point cloud score is, and the closer the score is to 0, the higher the floor point cloud score is. This represents the support attribute and the adjustment coefficient for the upper and lower cluster structures; the default value is 0.

5. This represents the horizontal topology weight adjustment coefficient; it takes a positive value when the support structure analysis and the horizontal topology consistency result are consistent, and a negative value when the hierarchical structure verification and the horizontal topology consistency result are consistent. This is achieved through dynamic adjustment. This ensures that the final decision satisfies the triple constraints of local geometry, vertical structure, and horizontal topology.

7. The method for hierarchical reconstruction of building point clouds incorporating multiple feature constraints according to claim 4, characterized in that, The process for setting the elevation merging threshold is as follows: Using a cloth-simulated mesh as the resolution, the 3D point cloud is divided according to the resolution. In planar space, based on the maximum and minimum values ​​of the imported 3D point cloud. The horizontal boundary of the corresponding 3D point cloud is obtained. Based on the horizontal boundary, the number and spatial distribution of the grids in each row and column are calculated. The initial plane of the virtual cloth is set at a fixed distance above the highest point of the target point cloud to ensure that the grid covers the entire point cloud area. For each grid after division, the maximum elevation value of the point cloud within the grid is extracted, and the maximum elevation value is used as the grid particle of the corresponding grid. The minimum elevation of the grid particle is used as the constraint threshold. Iterate through each movable grid particle and calculate the elevation of the grid particle under the influence of gravity using the following formula: In the formula: Indicates the time step. Indicates constant acceleration. This represents the descent distance threshold within a time step. Indicates the mesh particles in the current iteration Elevation value at time, Indicates grid particles in Elevation value at time, Indicates grid particles in Elevation value at time, This represents the instantaneous vertical velocity of the grid particle; after obtaining the elevation value of the grid particle, the elevation value of the grid particle is updated to obtain the grid particle elevation value, and compared with the minimum elevation constraint threshold. The comparison is performed when the grid particle elevation value is less than or equal to the minimum elevation constraint threshold. Then the elevation value of the grid particle will be reset to the minimum elevation constraint threshold. And mark it as immovable; For each movable mesh particle, extract the immovable mesh particles within its neighborhood. Based on the stiffness parameter and elevation difference, calculate the displacement adjustment under internal forces. The calculation formula is as follows: In the formula: This indicates the amount of displacement adjustment under the action of internal forces. This represents the elevation difference between movable and immovable grid particles. For fabric rigidity parameters; When the elevation of a movable grid particle is higher than that of a neighboring immovable grid particle, the corresponding displacement is adjusted along the negative elevation direction; otherwise, it is adjusted along the positive elevation direction. Sort all movable grid particles by the distance to the nearest immovable grid particle. Traverse the sorted grid particles and their neighboring movable particles. Adjust the elevation of adjacent movable particles by bidirectional displacement to make them more consistent. This is the elevation merging threshold.

8. The method for hierarchical reconstruction of building point clouds incorporating multiple feature constraints according to claim 7, characterized in that, The process of building a structured trend surface model is as follows: Traverse the neighborhood of the immovable mesh particle, and when the lowest elevation constraint threshold of the movable particle in the neighborhood is reached, [the following is a process to determine the minimum elevation constraint threshold]. If the difference is less than the elevation merging threshold, the particle is marked as immovable, and its elevation is updated to its lowest elevation constraint threshold. They are then added to an expanded queue, and the iteration is repeated until no new particles are added to the queue. Calculate the elevation difference between each original point cloud and its corresponding grid particle, and determine whether the elevation difference meets the elevation merging threshold. If it does, it is retained as the target plane and the target plane point set is constructed. Otherwise, the original point cloud is removed. Traverse the neighborhood of each immovable mesh particle and calculate its minimum elevation constraint threshold relative to other immovable particles in the neighborhood. Difference: If the difference exceeds the elevation merging threshold, check if there exists an elevation in the lower grid that matches the minimum elevation constraint threshold of the higher grid. If point clouds with elevation differences within the elevation merging threshold exist, these supplementary point clouds are merged into the target plane point set. Consecutive elevation intervals in the target plane point set whose adjacent elevation differences are less than the elevation merging threshold are merged to obtain several elevation intervals. The original cluster point clouds within each elevation interval are extracted and differential preprocessing is performed to obtain the ground elevation corresponding to the ground-type point cloud and the ceiling elevation corresponding to the ceiling-type point cloud. The ground elevation and ceiling elevation are then preprocessed to obtain the ground point cloud and ceiling point cloud, respectively. The ground point cloud and ceiling point cloud are then extracted separately. , The coordinate range is generated according to the predetermined grid spacing. , Planar meshes, with mesh nodes calculated via linear interpolation. The initial trend surface is constructed using the value, and then the wavefront algorithm is used to fill in the missing values ​​and missing areas in the trend surface. At the same time, the average elevation of each type of initial trend surface is recorded, and a trend surface is constructed that includes at least grid coordinates. A structured trend surface model of the original point cloud values.

9. The method for hierarchical reconstruction of building point clouds incorporating multiple feature constraints according to claim 4, characterized in that, The complete process of extracting the target plane is as follows: Using pairwise pairings of ground trend surfaces and ceiling trend surfaces as spatial indices, the point cloud between the ground trend surface and the ceiling trend surface represents the point cloud of the target floor. The point cloud for each pair of trend surfaces is then calculated. , The coordinate horizontal intersection is used to filter out the core point cloud points that fall within this intersection range and mark them as unextracted point clouds. Inverse interpolation is then used to obtain the coordinates of these unextracted point clouds on the corresponding ceiling and ground trend surfaces. The value is used to obtain the upper and lower limits of the floor elevation. Point clouds whose elevations fall within this range are selected as mezzanine point clouds, and their masks are marked as extracted point clouds, thus achieving complete extraction of the target plane.

Images