A method for repairing holes in a point cloud based on a nonlinear radial basis function

Abstract

The application discloses a kind of nonlinear radial basis function hole repair methods based on terrain point cloud, which comprises the following steps: one, obtaining the point cloud with hole to be handled;Two, obtaining point cloud average interval and closed hole boundary point set;Three, the closed hole boundary point set is handled, and true texture roughness and average edge slope angle are obtained;Four, the target excavation depth of hole center is obtained;Five, the reconstruction elevation value of each virtual feature point is obtained;Six, the repair elevation value of each virtual feature point is obtained;Seven, based on the distance between virtual feature point and closed hole boundary point set, boundary suture calculation is carried out and coordinate mapping is carried out, and the point cloud after repair is obtained.The method of the application deduces depth and nonlinear radial basis function interpolation reconstruction by intercept power function depth deduction model, and effectively eliminates the elevation mutation between repair area and original stable area.

Description

Technical Field

[0001] This invention belongs to the field of hole point cloud repair technology, and in particular relates to a nonlinear radial basis function hole repair method based on terrain point clouds. Background Technology

[0002] In recent years, with the rapid development of surveying and remote sensing technologies, three-dimensional information acquisition technologies, represented by photogrammetry, have provided point cloud data as a foundation for high-precision calculations and visualization in areas such as surface deformation research and disaster monitoring. However, due to the complex structure of valleys, severe terrain occlusion, and shadows or light attenuation, the acquired point cloud data often contains holes, which seriously restricts the high-precision representation of three-dimensional model information in valley areas.

[0003] Currently, the repair of holes in 3D point clouds mainly adopts spatial surface interpolation, Poisson surface reconstruction, contextual cloning and copying, and deep learning-based completion algorithms, but these methods have limitations when dealing with eroded terrain with concave morphological features. Repair methods based on spatial surface interpolation construct smooth surfaces using radial basis functions. However, due to their mathematical nature, they tend to generate minimal surfaces, making them prone to planar capping distortion when dealing with concave landforms, and unable to restore the true incision depth. Poisson surface reconstruction algorithms reconstruct closed surfaces by solving the gradient field of the index function. Although they can achieve hole closure, their global smoothing characteristics cause the repaired area to lose fine surface texture, resulting in obvious discontinuities in statistical characteristics such as roughness between the repaired area and the surrounding stable area. Contextual cloning and copying algorithms search and migrate geometrically similar blocks from known areas around the hole for filling. However, in erosion landforms where terrain evolution is unique, it is often difficult to find perfectly matching erosion textures, easily leading to unnatural geometric seams at the repair edges. Deep learning-based completion algorithms use architectures such as generative adversarial networks to learn terrain feature mappings. Although they have predictive capabilities, they lack physical laws and field constraints, so the generated point clouds often do not conform to the real water erosion laws and are limited by the availability of training samples.

[0004] Therefore, a nonlinear radial basis function hole repair method based on topographic point clouds is needed. This method obtains the set of closed hole boundary points from the point cloud with holes, and reconstructs the depth based on the closed hole boundary point set by using a depth extrapolation model with intercept power function and nonlinear radial basis function interpolation. Boundary stitching calculation is then performed, which effectively eliminates the elevation abrupt change between the repair area and the original stable area, achieves a high-precision natural transition between the repair surface and the original landform, and enhances the detail representation of the repair area. Summary of the Invention

[0005] The technical problem to be solved by this invention is to address the shortcomings of the prior art by providing a nonlinear radial basis function hole repair method based on topographic point clouds. The method obtains the set of closed hole boundary points from the point cloud with holes, and reconstructs the depth based on the closed hole boundary point set by using a depth extrapolation model with intercept power function and nonlinear radial basis function interpolation. Boundary stitching calculation is then performed, which effectively eliminates the elevation abrupt change between the repair area and the original stable area, achieves a high-precision natural transition between the repair surface and the original landform, and enhances the detail representation of the repair area.

[0006] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is: a nonlinear radial basis function hole repair method based on terrain point clouds, the method comprising the following steps: Step 1: Acquire and preprocess the terrain point cloud to obtain the point cloud with holes to be processed; Step 2: Traverse and cluster the point cloud with holes to obtain the average spacing ρ of the point cloud and the set of points at the boundary of the closed holes; Step 3: Process the boundary point set of the closed hole to obtain the true texture roughness σ and the average edge slope angle θavg; Step 4: Input the average edge slope angle θavg into the depth extrapolation model with intercept power function to extrapolate the depth radius ratio Rs, and obtain the target excavation depth Ht of the hole center based on the equivalent radius R of the hole; Step 5: Obtain the virtual feature points of the blank area surrounded by the set of boundary points of the closed hole, and subtract the vertical downcut constraint from the benchmark elevation obtained by nonlinear radial basis function interpolation to obtain the reconstructed elevation value of each virtual feature point; Step 6: Perform ground undulation perturbation optimization on the reconstructed elevation value of each virtual feature point to obtain the repaired elevation value of each virtual feature point; Step 7: Calculate the boundary stitching based on the distance between the virtual feature points and the boundary point set of the closed hole, and perform coordinate mapping to obtain the repaired point cloud.

[0007] The aforementioned nonlinear radial basis function hole repair method based on terrain point clouds is further refined in step one, as follows: Step 101: Use the SfM method to perform photogrammetric 3D reconstruction of the valley area to be measured to obtain a topographic point cloud with holes; wherein, the 3D coordinates of the point cloud are in the world coordinate system. Step 102: Use a computer to perform a first-pass filtering on the point cloud with holes in the terrain using the MCC point cloud filtering algorithm to obtain a first-pass filtered point cloud with holes. Step 103: Use a computer and TerraSolid software to perform a second filtering on the point cloud with holes after the first filtering, to obtain a point cloud with holes after the second filtering. Step 104: Input the point cloud with holes after secondary filtering into CloudCompare software for downsampling to obtain the point cloud with holes to be processed.

[0008] The aforementioned nonlinear radial basis function hole repair method based on terrain point clouds further includes step two, which is as follows: Step 201: Use the KDTree construction algorithm to construct the point cloud with holes to be processed, and obtain the point cloud binary tree; Step 202: Iterate through and calculate the Euclidean distance between any two points, sort all Euclidean distances, and take the median as the average spacing ρ of the point cloud. Step 203: Obtain the number of neighborhood search points k according to k = int{-[(kmax-kmin) / (ρmax-ρmin)]×(ρ-ρmin)+kmax}; where kmax is the maximum number of points, kmin is the minimum number of points, ρmax is the maximum spacing, and ρmin is the minimum spacing; int{ } indicates rounding down; Step 204: Use the nearest neighbor search algorithm to search the binary tree of the point cloud to obtain the k nearest neighbor points of point Pi, which constitute the local neighborhood point set of point Pi; Step 205: Perform principal component analysis on the three-dimensional coordinates of the local neighborhood point set of point Pi, extract the eigenvector corresponding to the minimum eigenvalue as the first normal vector, and fit a local micro-tangent plane perpendicular to the first normal vector. Step 206: Project point Pi and its k nearest neighbors orthogonally onto the local micro-tangent plane to obtain the two-dimensional centroid of the projected k nearest neighbors. When the Euclidean distance between the projected point of point Pi and the two-dimensional centroid is greater than the average spacing ρ of the point cloud, then point Pi is marked as the boundary point of the hole. Step 207: Repeat steps 204 to 206 multiple times to complete the traversal and judgment of each point and obtain the set of hole boundary points; Step 208: Use the DBSCAN clustering algorithm to cluster the set of hole boundary points to obtain the set of closed hole boundary points.

[0009] The aforementioned nonlinear radial basis function hole repair method based on terrain point clouds is further refined in step three, which is as follows: Step 301: Perform principal component analysis on the three-dimensional coordinates of the boundary point set of the closed hole, extract the eigenvector corresponding to the minimum eigenvalue as the second normal vector, and fit a local reference plane perpendicular to the second normal vector; wherein, a local coordinate system is established at the center of the local reference plane, the Z-axis of the local coordinate system is along the second normal vector, the X-axis and Y-axis are on the local reference plane and are perpendicular to each other, and both the X-axis and Y-axis are perpendicular to the Z-axis; Step 302: Obtain the distance from each point in the closed hole boundary point set to the local reference plane, and record it as the elevation residual of each point; Step 303: Take the median of the elevation residuals from all points as the true texture roughness σ; Step 304: Obtain the angle between the Z-axis and the second normal vector in the world coordinate system, denoted as the average edge slope angle θavg.

[0010] The aforementioned nonlinear radial basis function hole repair method based on terrain point clouds further includes step four, which is as follows: Step 401: Use the depth extrapolation model with intercept power function to extrapolate the depth-radius ratio Rs, where Rs = c′ + a′(θavg / 45°). b′ Where c′ is the base depth-radius ratio, a′ is the slope growth coefficient, and b′ is the nonlinear growth index; Step 402: Average the three-dimensional coordinates of the closed hole boundary point set to obtain the geometric center of the closed hole boundary point set; Step 403: Obtain the distance from each point in the set of boundary points of the closed hole to the geometric center, and take the median of all distances as the equivalent radius R of the hole; Step 404: Based on Ht=R×Rs, obtain the target excavation depth Ht at the center of the hole.

[0011] The aforementioned nonlinear radial basis function hole repair method based on terrain point clouds, further, step five, is as follows: Step 501: Orthogonally project each point of the closed hole boundary point set onto the local reference plane to obtain each boundary projection point; Step 502: Establish a grid within the blank area enclosed by the projection points of the boundary on the local reference plane, and record the corner points of each grid as virtual feature points, and obtain the two-dimensional coordinates of each virtual feature point; wherein, the length of the grid is the average spacing ρ of the point cloud; Step 503: According to Zdrop=Ht×[arctan(μ×ds / R) / arctan(μ)], obtain the vertical downcut constraint Zdrop for each virtual feature point; where ds is the shortest distance from each virtual feature point to the boundary projection point, and μ is the morphological convergence coefficient; Step 504: Based on the set of boundary points of the closed hole, input the two-dimensional coordinates of each virtual feature point, and use the nonlinear radial basis function of the thin plate spline kernel to perform spatial interpolation on each virtual feature point to obtain the reference elevation of each virtual feature point. Step 505: Subtract the vertical downcut constraint of each virtual feature point from the base elevation of each virtual feature point to obtain the reconstructed elevation value of each virtual feature point.

[0012] The aforementioned nonlinear radial basis function hole repair method based on terrain point clouds, further, step six, is as follows: Step 601: Generate a standard normal distribution using a random number generation algorithm, and take the random variable in the interval [-1,1] of the horizontal axis as Gaussian noise; Step 602: Multiply the Gaussian noise by the actual texture roughness σ to obtain the disturbance component of the simulated ground undulation. ; Step 603: Reconstruct the elevation value of each virtual feature point and the disturbance component of the simulated ground undulation. The values ​​are randomly summed to obtain the repaired elevation value for each virtual feature point.

[0013] The aforementioned nonlinear radial basis function hole repair method based on terrain point clouds, further, step seven, is as follows: Step 701: Based on tj=min(dj / 0.25R,1), obtain the normalization adjustment parameter tj of the j-th virtual feature point; where dj represents the shortest distance from the j-th virtual feature point to the boundary projection point; j is a positive integer, and the value of j ranges from 1 to J, where J is the total number of virtual feature points; Step 702: According to the Smoothstep weight function Wj=3(tj) 2 2 (tj) 3 The fusion weight Wj of the j-th virtual feature point is obtained; Step 703, according to Zfj=Wj×Zxj+(1 We obtain the fused elevation value Zfj of the j-th virtual feature point by multiplying Wj by Zb; where Zxj is the restored elevation value of the j-th virtual feature point, and Zb is the elevation value of the boundary projection point closest to the j-th virtual feature point. Step 704: Combine the two-dimensional coordinates of the j-th virtual feature point with the fused elevation value to obtain the local coordinates of the j-th virtual feature point; Step 705: Obtain the homogeneous transformation matrix of the local coordinate system relative to the world coordinate system; Step 706: Transform the local coordinates of the j-th virtual feature point using a homogeneous transformation matrix to obtain the three-dimensional coordinates of the j-th virtual feature point in the world coordinate system; Step 707: Repeat steps 701 to 706 multiple times to obtain the three-dimensional coordinates of each virtual feature point in the world coordinate system, and then obtain the repaired point cloud.

[0014] Compared with the prior art, the present invention has the following advantages: 1. This invention overcomes the problems of planar capping distortion, texture loss and unnatural geometric seams that are easily generated when the traditional point cloud hole interpolation repair algorithm is used to process complex concave landforms (such as valley areas).

[0015] 2. This invention introduces a local micro-tangent plane orthogonal projection technique based on point cloud average spacing evaluation and principal component analysis. It aims to reduce the dimensionality of local neighborhood point sets in drastically undulating terrain, effectively eliminating the interference caused by drastic elevation changes, greatly reducing misjudgment of hole boundary points in complex terrain, and improving the accuracy of subsequent automatic extraction of hole boundaries based on DBSCAN clustering.

[0016] 3. This invention extracts texture features based on local reference plane fitting, successfully removing the interference of macro slope on the calculation of micro surface roughness, and extracting pure physical roughness features. This allows the subsequently synthesized Gaussian texture to highly reproduce the micro texture of natural soil, avoiding the visual and physical separation between the restored area and the original landform.

[0017] 4. This invention uses a depth extrapolation model with intercept power function to extrapolate the target depth, and introduces the arctan operator to obtain the vertical downcut constraint quantity to construct a funnel geometric constraint model. This model can fit the physical law of water flow downcut in natural landforms (such as water erosion gullies), effectively reconstructing an accurate macroscopic concave trend surface, and avoiding the unnatural flat surface formed above the hole by traditional radial basis function interpolation.

[0018] 5. This invention constructs an initial reconstructed surface based on a nonlinear radial basis function interpolation algorithm, aiming to simulate the complex morphology of natural landforms through high-order continuous mathematical surfaces. Compared with traditional linear interpolation or simple polynomial fitting, this technique can automatically capture and continue the slope changes and terrain orientation of the hole edges, effectively solving the problem of planar capping and texture inconsistency that easily occurs on the reconstructed surface under complex curvature landforms, so that the repaired area achieves a natural and smooth connection with the surrounding original environment in terms of geometric topology.

[0019] 6. This invention addresses the problem of splicing marks between virtual feature points on the repair surface and the boundary point set of closed holes on the edge of the original landform. It uses the Smoothstep weight function, i.e., the nonlinear distance weight function, to perform boundary stitching calculation, effectively eliminating the elevation abrupt change between the repair area and the original stable area, realizing a natural transition between the repair surface and the original landform, and enhancing the detail representation of the repair area.

[0020] In summary, the method of this invention obtains the set of closed hole boundary points from point clouds with holes, and reconstructs the depth based on the closed hole boundary point set by using a depth extrapolation model with intercept power function and nonlinear radial basis function interpolation, and performs boundary stitching calculation. This effectively eliminates the elevation abrupt change between the repair area and the original stable area, achieves a high-precision natural transition between the repair surface and the original landform, and enhances the detail representation of the repair area.

[0021] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0022] Figure 1 This is a flowchart of the method of the present invention.

[0023] Figure 2 This is a schematic diagram of the topographic point cloud with holes according to the present invention.

[0024] Figure 3 This is a comparison diagram of the repair effects of the method of the present invention and existing methods.

[0025] Figure 4 Comparison images showing the cropped and projected areas of the hole repair area in this invention.

[0026] Figure 5 This is a comparison chart of the roughness repaired by the method of the present invention and existing methods, based on statistical analysis. Detailed Implementation

[0027] like Figure 1 The method for repairing holes based on nonlinear radial basis functions of terrain point clouds, as shown, includes the following steps: Step 1: Acquire and preprocess the terrain point cloud to obtain the point cloud with holes to be processed; Step 2: Traverse and cluster the point cloud with holes to obtain the average spacing ρ of the point cloud and the set of points at the boundary of the closed holes; Step 3: Process the boundary point set of the closed hole to obtain the true texture roughness σ and the average edge slope angle θavg; Step 4: Input the average edge slope angle θavg into the depth extrapolation model with intercept power function to extrapolate the depth radius ratio Rs, and obtain the target excavation depth Ht of the hole center based on the equivalent radius R of the hole; Step 5: Obtain the virtual feature points of the blank area surrounded by the set of boundary points of the closed hole, and subtract the vertical downcut constraint from the benchmark elevation obtained by nonlinear radial basis function interpolation to obtain the reconstructed elevation value of each virtual feature point; Step 6: Perform ground undulation perturbation optimization on the reconstructed elevation value of each virtual feature point to obtain the repaired elevation value of each virtual feature point; Step 7: Calculate the boundary stitching based on the distance between the virtual feature points and the boundary point set of the closed hole, and perform coordinate mapping to obtain the repaired point cloud.

[0028] In this embodiment, step one, the specific process is as follows: Step 101: Use the SfM method to perform photogrammetric 3D reconstruction of the valley area to be measured to obtain a topographic point cloud with holes; wherein, the 3D coordinates of the point cloud are in the world coordinate system. Step 102: Use a computer to perform a first-pass filtering on the point cloud with holes in the terrain using the MCC point cloud filtering algorithm to obtain a first-pass filtered point cloud with holes. Step 103: Use a computer and TerraSolid software to perform a second filtering on the point cloud with holes after the first filtering, to obtain a point cloud with holes after the second filtering. Step 104: Input the point cloud with holes after secondary filtering into CloudCompare software for downsampling to obtain the point cloud with holes to be processed.

[0029] In this embodiment, step two is as follows: Step 201: Use the KDTree construction algorithm to construct the point cloud with holes to be processed, and obtain the point cloud binary tree; Step 202: Iterate through and calculate the Euclidean distance between any two points, sort all Euclidean distances, and take the median as the average spacing ρ of the point cloud. Step 203: Obtain the number of neighborhood search points k according to k = int{-[(kmax-kmin) / (ρmax-ρmin)]×(ρ-ρmin)+kmax}; where kmax is the maximum number of points, kmin is the minimum number of points, ρmax is the maximum spacing, and ρmin is the minimum spacing; int{ } indicates rounding down; Step 204: Use the nearest neighbor search algorithm to search the binary tree of the point cloud to obtain the k nearest neighbor points of point Pi, which constitute the local neighborhood point set of point Pi; Step 205: Perform principal component analysis on the three-dimensional coordinates of the local neighborhood point set of point Pi, extract the eigenvector corresponding to the minimum eigenvalue as the first normal vector, and fit a local micro-tangent plane perpendicular to the first normal vector. Step 206: Project point Pi and its k nearest neighbors orthogonally onto the local micro-tangent plane to obtain the two-dimensional centroid of the projected k nearest neighbors. When the Euclidean distance between the projected point of point Pi and the two-dimensional centroid is greater than the average spacing ρ of the point cloud, then point Pi is marked as the boundary point of the hole. Step 207: Repeat steps 204 to 206 multiple times to complete the traversal and judgment of each point and obtain the set of hole boundary points; Step 208: Use the DBSCAN clustering algorithm to cluster the set of hole boundary points to obtain the set of closed hole boundary points.

[0030] In this embodiment, steps 201 to 208 introduce the local micro-tangent plane orthogonal projection technique of point cloud average spacing evaluation and principal component analysis. This technique aims to reduce the dimensionality of local neighborhood point sets in drastically undulating terrain, effectively eliminating the interference caused by drastic elevation changes, greatly reducing the misjudgment of hole boundary points in complex landforms, and improving the accuracy of subsequent automatic extraction of hole boundaries based on DBSCAN clustering.

[0031] In this embodiment, step three is as follows: Step 301: Perform principal component analysis on the three-dimensional coordinates of the boundary point set of the closed hole, extract the eigenvector corresponding to the minimum eigenvalue as the second normal vector, and fit a local reference plane perpendicular to the second normal vector; wherein, a local coordinate system is established at the center of the local reference plane, the Z-axis of the local coordinate system is along the second normal vector, the X-axis and Y-axis are on the local reference plane and are perpendicular to each other, and both the X-axis and Y-axis are perpendicular to the Z-axis; Step 302: Obtain the distance from each point in the closed hole boundary point set to the local reference plane, and record it as the elevation residual of each point; Step 303: Take the median of the elevation residuals from all points as the true texture roughness σ; Step 304: Obtain the angle between the Z-axis and the second normal vector in the world coordinate system, denoted as the average edge slope angle θavg.

[0032] In this embodiment, step four is as follows: Step 401: Use the depth extrapolation model with intercept power function to extrapolate the depth-radius ratio Rs, where Rs = c′ + a′(θavg / 45°). b′ Where c′ is the base depth-radius ratio, a′ is the slope growth coefficient, and b′ is the nonlinear growth index; Step 402: Average the three-dimensional coordinates of the closed hole boundary point set to obtain the geometric center of the closed hole boundary point set; Step 403: Obtain the distance from each point in the set of boundary points of the closed hole to the geometric center, and take the median of all distances as the equivalent radius R of the hole; Step 404: Based on Ht=R×Rs, obtain the target excavation depth Ht at the center of the hole.

[0033] In this embodiment, step five is as follows: Step 501: Orthogonally project each point of the closed hole boundary point set onto the local reference plane to obtain each boundary projection point; Step 502: Establish a grid within the blank area enclosed by the projection points of the boundary on the local reference plane, and record the corner points of each grid as virtual feature points, and obtain the two-dimensional coordinates of each virtual feature point; wherein, the length of the grid is the average spacing ρ of the point cloud; Step 503: According to Zdrop=Ht×[arctan(μ×ds / R) / arctan(μ)], obtain the vertical downcut constraint Zdrop for each virtual feature point; where ds is the shortest distance from each virtual feature point to the boundary projection point, and μ is the morphological convergence coefficient; Step 504: Based on the set of boundary points of the closed hole, input the two-dimensional coordinates of each virtual feature point, and use the nonlinear radial basis function of the thin plate spline kernel to perform spatial interpolation on each virtual feature point to obtain the reference elevation of each virtual feature point. Step 505: Subtract the vertical downcut constraint of each virtual feature point from the base elevation of each virtual feature point to obtain the reconstructed elevation value of each virtual feature point.

[0034] In this embodiment, the target depth is derived using a depth extrapolation model with intercept power function in step 401. In step 503, the vertical downcut constraint is obtained by introducing the arctan operator, and a funnel geometric constraint model is constructed. This model can conform to the physical law of water flow downcutting in natural landforms (such as water erosion gullies), effectively reconstructing an accurate macroscopic concave trend surface, and avoiding the unnatural flat surface formed above the hole by traditional radial basis function interpolation.

[0035] In this embodiment, step 504 constructs an initial reconstruction surface based on a nonlinear radial basis function interpolation algorithm, aiming to simulate the complex morphology of natural landforms through a high-order continuous mathematical surface. Compared to traditional linear interpolation or simple polynomial fitting, this technique can automatically capture and continue the slope changes and terrain orientation of the hole edges, effectively solving the problem of planar capping and texture inconsistency that easily occurs on the reconstructed surface under complex curvature landforms, so that the repaired area achieves a natural and smooth connection with the surrounding original environment in terms of geometric topology.

[0036] In this embodiment, step six is ​​as follows: Step 601: Generate a standard normal distribution using a random number generation algorithm, and take the random variable in the interval [-1,1] of the horizontal axis as Gaussian noise; Step 602: Multiply the Gaussian noise by the actual texture roughness σ to obtain the disturbance component of the simulated ground undulation. ; Step 603: Reconstruct the elevation value of each virtual feature point and the disturbance component of the simulated ground undulation. The values ​​are randomly summed to obtain the repaired elevation value for each virtual feature point.

[0037] In this embodiment, step seven is as follows: Step 701: Based on tj=min(dj / 0.25R,1), obtain the normalization adjustment parameter tj of the j-th virtual feature point; where dj represents the shortest distance from the j-th virtual feature point to the boundary projection point; j is a positive integer, and the value of j ranges from 1 to J, where J is the total number of virtual feature points; Step 702: According to the Smoothstep weight function Wj=3(tj) 2 2 (tj) 3 The fusion weight Wj of the j-th virtual feature point is obtained; Step 703, according to Zfj=Wj×Zxj+(1 We obtain the fused elevation value Zfj of the j-th virtual feature point by multiplying Wj by Zb; where Zxj is the restored elevation value of the j-th virtual feature point, and Zb is the elevation value of the boundary projection point closest to the j-th virtual feature point. Step 704: Combine the two-dimensional coordinates of the j-th virtual feature point with the fused elevation value to obtain the local coordinates of the j-th virtual feature point; Step 705: Obtain the homogeneous transformation matrix of the local coordinate system relative to the world coordinate system; Step 706: Transform the local coordinates of the j-th virtual feature point using a homogeneous transformation matrix to obtain the three-dimensional coordinates of the j-th virtual feature point in the world coordinate system; Step 707: Repeat steps 701 to 706 multiple times to obtain the three-dimensional coordinates of each virtual feature point in the world coordinate system, and then obtain the repaired point cloud.

[0038] In this embodiment, in steps 701 to 706, to address the issue of splicing marks between the virtual feature points of the repair surface and the boundary point set of closed holes on the edge of the original landform, the Smoothstep weight function, i.e., the nonlinear distance weight function, is used to perform boundary stitching calculations. This effectively eliminates the elevation abrupt change between the repair area and the original stable area, achieves a natural transition between the repair surface and the original landform, and enhances the detail representation of the repair area.

[0039] In this embodiment, the repaired point cloud and the terrain point cloud with holes will be merged to obtain a terrain point cloud.

[0040] In this embodiment, the elevation value Zb of the boundary projection point closest to the j-th virtual feature point in step 703 is obtained as follows: First, obtain the point corresponding to the boundary projection point closest to the j-th virtual feature point from the set of closed hole boundary points, and then obtain the distance from the point to the local reference plane. This distance is the elevation value Zb.

[0041] In this embodiment, the world coordinate system is the Northeast-Sky coordinate system. The three-dimensional coordinates of the terrain point cloud are based on the Northeast-Sky coordinate system.

[0042] In this embodiment, SfM stands for Structure from Motion, which is a photogrammetric technique that recovers the structure of a three-dimensional scene by analyzing image sequences. It is a conventional method in the field, such as inputting valley photos acquired by camera measurement into Context Capture software to obtain topographic point clouds with holes.

[0043] In this embodiment, it should be noted that the minimum spacing parameter between point clouds is set to 0.005 meters during downsampling in step 104.

[0044] In this embodiment, the maximum number of points kmax is set to 50, the minimum number of points kmin is set to 10, the maximum spacing ρmax is set to 0.1 meters, and the minimum spacing ρmin is set to 0.005 meters.

[0045] In this embodiment, it should be noted that DBSCAN stands for Density-Based Spatial Clustering of Applications with Noise, a density-based spatial clustering algorithm.

[0046] In this embodiment, it should be noted that c′ is the basic depth-radius ratio, which is used to limit the lower physical limit of landform evolution; a′ is the slope growth coefficient, which is used to characterize the increase of the hole incision morphology as the slope increases; b′ is the nonlinear growth index, which controls the degree of accelerated development of the inferred curve as the slope changes. The above parameters can be statistically fitted and calibrated based on the erosion resistance characteristics of different geological soil types (such as loess, red soil, etc.). In a preferred embodiment of the present invention, based on the statistical characteristics of typical samples in the loess hilly valley area, c′=1.002, a′=0.088, b′=6.321 are preferred.

[0047] In this embodiment, it should be noted that μ is the morphological convergence coefficient, which is used to control the steepness or gentleness of the reconstructed funnel profile (such as U-shaped or V-shaped cross-section distribution). For the typical downcutting waterfall morphology of loess gullies, the morphological convergence coefficient μ is preferably set to 3.0.

[0048] In this embodiment, it should be noted that the nonlinear radial basis function of the thin plate spline kernel is used for spatial interpolation, and the smoothing coefficient is set to 0.01.

[0049] To verify the effectiveness of the present invention, this embodiment selects a valley terrain with typical concave features for comparative experiment.

[0050] Experimental data preparation: A simulated water erosion experiment was conducted on the gully slope. Each erosion session lasted 15 minutes. After each 15-minute erosion session, the topography of the gully slope changed, such as... Figure 2 The SfM method was used to obtain the third field of the terrain point cloud with holes in the region. Ground-based 3D laser scanning was used to obtain the terrain point cloud without holes, and the terrain point clouds with and without holes were filtered, denoised, and registered.

[0051] Comparison scheme setup: Ground-based 3D laser scanning is used to acquire hole-free terrain point clouds, denoted as ground-based 3D laser scanning point clouds, such as... Figure 3 As shown in Figure a, the hole in the point cloud with hole in the terrain was repaired using the method of this invention, the traditional RBF interpolation method, and the Poisson surface reconstruction algorithm, respectively. The results are as follows. Figure 3 As shown in b, c, and d, b represents the repair result of the method of the present invention, which accurately restores the concave shape and texture; c represents the repair result of the traditional RBF interpolation method, which shows obvious planar capping phenomenon; and d represents the repair result of the Poisson surface reconstruction algorithm, which shows over-smoothing characteristics.

[0052] Furthermore, using the ground-based 3D laser scan point cloud as the reference ground truth, the quality of the restoration results was verified through the following two comparison methods: 1. Macroscopic geometric topological verification: such as Figure 4 As shown, taking the point cloud of the terrain with holes in the third field as an example, the hole repair area is clipped and projected. The results show that the projection point A of the traditional RBF interpolation method is a straight line spanning the hole, and the projection point B of the Poisson surface reconstruction algorithm, although it has some downcutting, lacks sufficient depth. However, the projection point D of the present invention can more accurately fit the projection point C of the ground 3D laser scanning point cloud, with a maximum local vertical deviation of only 12mm. From a quantitative perspective, this proves that the present invention effectively eliminates the interpolation distortion of concave terrain.

[0053] 2. Micro-statistical characteristic tests: such as Figure 5 As shown, the horizontal axis represents the field number, and the vertical axis represents the average roughness. Roughness statistical analysis was performed on the microscopic geomorphological features of the restored areas in five fields. The results show that, taking field 3 as an example, the average roughness Ra3 of the non-porosity area of ​​the point cloud is 28 μm. The average roughness Ra1 of the restored area using the traditional RBF interpolation method is 67 μm, much higher than the average roughness of the non-porosity area, while the average roughness Ra4 of the restored area using the Poisson surface reconstruction algorithm is 10 μm, lower than the average roughness of the non-porosity area. In contrast, the average roughness Ra2 of the restored area using the method of this invention is 31 μm, close to the average roughness Ra3 of the non-porosity area, maintaining consistency. Similar conclusions were reached for other fields, demonstrating the advantages of this invention in preserving geomorphological statistical features.

[0054] In summary, the method of this invention obtains the set of closed hole boundary points from point clouds with holes, and reconstructs the depth based on the closed hole boundary point set by using a depth extrapolation model with intercept power function and nonlinear radial basis function interpolation, and performs boundary stitching calculation. This effectively eliminates the elevation abrupt change between the repair area and the original stable area, achieves a high-precision natural transition between the repair surface and the original landform, and enhances the detail representation of the repair area.

[0055] The above description is merely a preferred embodiment of the present invention and does not constitute any limitation on the present invention. Any simple modifications, alterations, or equivalent structural changes made to the above embodiments based on the technical essence of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A nonlinear radial basis function hole repair method based on terrain point clouds, characterized in that, The method includes the following steps: Step 1: Acquire and preprocess the terrain point cloud to obtain the point cloud with holes to be processed; Step 2: Traverse and cluster the point cloud with holes to obtain the average spacing ρ of the point cloud and the set of points at the boundary of the closed holes; Step 3: Process the boundary point set of the closed hole to obtain the true texture roughness σ and the average edge slope angle θavg; Step 4: Input the average edge slope angle θavg into the depth extrapolation model with intercept power function to extrapolate the depth radius ratio Rs, and obtain the target excavation depth Ht of the hole center based on the equivalent radius R of the hole; Step 5: Obtain the virtual feature points of the blank area surrounded by the set of boundary points of the closed hole, and subtract the vertical downcut constraint from the benchmark elevation obtained by nonlinear radial basis function interpolation to obtain the reconstructed elevation value of each virtual feature point; Step 6: Perform ground undulation perturbation optimization on the reconstructed elevation value of each virtual feature point to obtain the repaired elevation value of each virtual feature point; Step 7: Calculate the boundary stitching based on the distance between the virtual feature points and the boundary point set of the closed hole, and perform coordinate mapping to obtain the repaired point cloud.

2. The nonlinear radial basis function hole repair method based on terrain point clouds according to claim 1, characterized in that: Step one, the specific process is as follows: Step 101: Use the SfM method to perform photogrammetric 3D reconstruction of the valley area to be measured to obtain a topographic point cloud with holes; wherein, the 3D coordinates of the point cloud are in the world coordinate system. Step 102: Use a computer to perform a first-pass filtering on the point cloud with holes in the terrain using the MCC point cloud filtering algorithm to obtain a first-pass filtered point cloud with holes. Step 103: Use a computer and TerraSolid software to perform a second filtering on the point cloud with holes after the first filtering, to obtain a point cloud with holes after the second filtering. Step 104: Input the point cloud with holes after secondary filtering into CloudCompare software for downsampling to obtain the point cloud with holes to be processed.

3. The nonlinear radial basis function hole repair method based on terrain point clouds according to claim 2, characterized in that: Step two, the specific process is as follows: Step 201: Use the KDTree construction algorithm to construct the point cloud with holes to be processed, and obtain the point cloud binary tree; Step 202: Iterate through and calculate the Euclidean distance between any two points, sort all Euclidean distances, and take the median as the average spacing ρ of the point cloud. Step 203: Obtain the number of neighborhood search points k according to k = int{-[(kmax-kmin) / (ρmax-ρmin)]×(ρ-ρmin)+kmax}; where kmax is the maximum number of points, kmin is the minimum number of points, ρmax is the maximum spacing, and ρmin is the minimum spacing; int{ } indicates rounding down; Step 204: Use the nearest neighbor search algorithm to search the binary tree of the point cloud to obtain the k nearest neighbor points of point Pi, which constitute the local neighborhood point set of point Pi; Step 205: Perform principal component analysis on the three-dimensional coordinates of the local neighborhood point set of point Pi, extract the eigenvector corresponding to the minimum eigenvalue as the first normal vector, and fit a local micro-tangent plane perpendicular to the first normal vector. Step 206: Project point Pi and its k nearest neighbors orthogonally onto the local micro-tangent plane to obtain the two-dimensional centroid of the projected k nearest neighbors. When the Euclidean distance between the projected point of point Pi and the two-dimensional centroid is greater than the average spacing ρ of the point cloud, then point Pi is marked as the boundary point of the hole. Step 207: Repeat steps 204 to 206 multiple times to complete the traversal and judgment of each point and obtain the set of hole boundary points; Step 208: Use the DBSCAN clustering algorithm to cluster the set of hole boundary points to obtain the set of closed hole boundary points.

4. The nonlinear radial basis function hole repair method based on terrain point clouds according to claim 3, characterized in that: Step three, the specific process is as follows: Step 301: Perform principal component analysis on the three-dimensional coordinates of the boundary point set of the closed hole, extract the eigenvector corresponding to the minimum eigenvalue as the second normal vector, and fit a local reference plane perpendicular to the second normal vector; wherein, a local coordinate system is established at the center of the local reference plane, the Z-axis of the local coordinate system is along the second normal vector, the X-axis and Y-axis are on the local reference plane and are perpendicular to each other, and both the X-axis and Y-axis are perpendicular to the Z-axis; Step 302: Obtain the distance from each point in the closed hole boundary point set to the local reference plane, and record it as the elevation residual of each point; Step 303: Take the median of the elevation residuals from all points as the true texture roughness σ; Step 304: Obtain the angle between the Z-axis and the second normal vector in the world coordinate system, denoted as the average edge slope angle θavg.

5. A nonlinear radial basis function hole repair method based on terrain point clouds according to claim 1, characterized in that: Step four, the specific process is as follows: Step 401: Use the depth extrapolation model with intercept power function to extrapolate the depth-radius ratio Rs, where Rs = c′ + a′(θavg / 45°). b′ Where c′ is the base depth-radius ratio, a′ is the slope growth coefficient, and b′ is the nonlinear growth index; Step 402: Average the three-dimensional coordinates of the closed hole boundary point set to obtain the geometric center of the closed hole boundary point set; Step 403: Obtain the distance from each point in the set of boundary points of the closed hole to the geometric center, and take the median of all distances as the equivalent radius R of the hole; Step 404: Based on Ht=R×Rs, obtain the target excavation depth Ht at the center of the hole.

6. A nonlinear radial basis function hole repair method based on terrain point clouds according to claim 4, characterized in that: Step five, the specific process is as follows: Step 501: Orthogonally project each point of the closed hole boundary point set onto the local reference plane to obtain each boundary projection point; Step 502: Establish a grid within the blank area enclosed by the projection points of the boundary on the local reference plane, and record the corner points of each grid as virtual feature points, and obtain the two-dimensional coordinates of each virtual feature point; wherein, the length of the grid is the average spacing ρ of the point cloud; Step 503: According to Zdrop=Ht×[arctan(μ×ds / R) / arctan(μ)], obtain the vertical downcut constraint Zdrop for each virtual feature point; where ds is the shortest distance from each virtual feature point to the boundary projection point, and μ is the morphological convergence coefficient; Step 504: Based on the set of boundary points of the closed hole, input the two-dimensional coordinates of each virtual feature point, and use the nonlinear radial basis function of the thin plate spline kernel to perform spatial interpolation on each virtual feature point to obtain the reference elevation of each virtual feature point. Step 505: Subtract the vertical downcut constraint of each virtual feature point from the base elevation of each virtual feature point to obtain the reconstructed elevation value of each virtual feature point.

7. A nonlinear radial basis function hole repair method based on terrain point clouds according to claim 1, characterized in that: Step six, the specific process is as follows: Step 601: Generate a standard normal distribution using a random number generation algorithm, and take the random variable in the interval [-1,1] of the horizontal axis as Gaussian noise; Step 602: Multiply the Gaussian noise by the actual texture roughness σ to obtain the disturbance component of the simulated ground undulation. ; Step 603: Reconstruct the elevation value of each virtual feature point and the disturbance component of the simulated ground undulation. The values ​​are randomly summed to obtain the repaired elevation value for each virtual feature point.

8. A nonlinear radial basis function hole repair method based on terrain point clouds according to claim 6, characterized in that: Step seven, the specific process is as follows: Step 701: Based on tj=min(dj / 0.25R,1), obtain the normalization adjustment parameter tj of the j-th virtual feature point; where dj represents the shortest distance from the j-th virtual feature point to the boundary projection point; j is a positive integer, and the value of j ranges from 1 to J, where J is the total number of virtual feature points; Step 702: According to the Smoothstep weight function Wj=3(tj) 2 2 (tj) 3 The fusion weight Wj of the j-th virtual feature point is obtained; Step 703, according to Zfj=Wj×Zxj+(1 We obtain the fused elevation value Zfj of the j-th virtual feature point by multiplying Wj by Zb; where Zxj is the restored elevation value of the j-th virtual feature point, and Zb is the elevation value of the boundary projection point closest to the j-th virtual feature point. Step 704: Combine the two-dimensional coordinates of the j-th virtual feature point with the fused elevation value to obtain the local coordinates of the j-th virtual feature point; Step 705: Obtain the homogeneous transformation matrix of the local coordinate system relative to the world coordinate system; Step 706: Transform the local coordinates of the j-th virtual feature point using a homogeneous transformation matrix to obtain the three-dimensional coordinates of the j-th virtual feature point in the world coordinate system; Step 707: Repeat steps 701 to 706 multiple times to obtain the three-dimensional coordinates of each virtual feature point in the world coordinate system, and then obtain the repaired point cloud.

Images