HOLE-FILLING ALGORITHM FOR AIRBORNE LIDAR POINT CLOUD DATA

： Due to the influence of the occlusion of objects or the complexity of the measured terrain in the scanning process of airborne lidar, the point cloud data inevitably appears holes after filtering and other processing. The incomplete data will inevitably have an impact on the quality of the reconstructed digital elevation model, so how to repair the incomplete point cloud data has become an urgent problem to be solved. To solve the problem of hole repair in point cloud data, a hole repair algorithm based on improved moving least square method is proposed in this paper by studying existing hole repair algorithms. Firstly, the algorithm extracts the boundary of the point cloud based on the triangular mesh model. Then we use k-nearest neighbor search to obtain the k-nearest neighbor points of the boundary point. Finally, according to the boundary point and its k-nearest neighbor point, the improved moving least squares method is used to fit the hole surface to realize the hole repair. Combined with C++ and MATLAB language, the feasibility of the algorithm is tested by specific application examples. The experimental results show that the algorithm can effectively repair the point cloud data holes, and the repairing precision is high. The filled hole area can be smoothly connected with the boundary.


INTRODUCTION
In recent years, airborne lidar measurement technology has become one of the new technologies in the field of photogrammetry, which has high efficiency of obtaining spatial information, high precision and certain penetration of the emitted laser pulse to vegetation (Zhou G, 2013 and2014). With the development of airborne lidar, in the aspect of digital elevation model (DEM) data acquisition, airborne lidar system can quickly obtain the digital surface model (DSM) of the area (Zhou G, 2018). However, due to the occlusion of ground objects, the limitation of measuring equipment and the complexity of the shape of the object to be measured, the collected point cloud data often contain a variety of areas that can not be measured in the process of collecting data by airborne lidar, and there are a large number of holes. The existence of hole area in point cloud data will directly lead to voids in triangular grid model generated by point cloud. These holes will not only affect the visual effect of 3D model, but also affect a series of operations based on grid model in the later stage. Therefore, point cloud hole repair is very important in the preprocessing process of point cloud data.
In order to solve the problem of hole repair in point cloud model, scholars at home and abroad have made a deep research on it. At present, the hole repair technology at home and abroad is mainly divided into two methods: the first is to repair the discrete point cloud directly, the second is to repair the point cloud based on the grid model, among which there are three main hole repair methods based on the grid model: (1) voxel-based methods; (2) based on triangulation method (3) based on surface fitting method. In the voxel-based method, Joshua and Szymon (2005) split the space into internal and external parts by using the minimum cutting method, and repaired the holes in a globally sensitive manner. Ju Tao (2004) used octree data structure to reconstruct surfaces and repair holes by constructing internal or external volumes and contours. In the triangulation-based method, Chen et al. (2013) used the method of interpolation to calculate the shape and density of the triangular patches around the point cloud cavity. Yongtae Jun (2005) proposed an algorithm that divides a complex hole into simple holes and then repairs them. The method based on surface fitting is the most common in the mesh hole model patching algorithm. Levin (1999) proposed a hole segmentation algorithm based on 1 continuous subdivision surface. G. Casciola (2005) and D. Lazzaro (2002) used radial basis functions to construct curved surfaces for hole filling. J Wang et al (2003) and Zhang (2017) proposed a hole filling algorithm based on moving least squares method, which uses moving least squares to fit a quadric surface for interpolation in the local coordinate system. In the abovementioned point cloud hole repair algorithm, the voxel-based method can effectively repair complex holes, but the repair takes a long time and may generate a wrong topological relationship. The triangulation-based method can quickly fill the holes, but only the boundary points of the holes are considered in the repair, which cannot be smoothly connected with the original mesh model. The hole filling algorithm based on moving least squares can repair the holes more accurately, and the repaired point cloud surface can be more consistent with the original terrain.
However, when using this method to construct the shape function, it is easy to get the solution matrix is ill, even singular, so that it is difficult to solve or get the correct solution.
Our research team won a special focus of Guangxi Innovative Development Grand Grant in 2018. The purpose of this project is to develop a set of airborne laser high-precision threedimensional seabed measuring instruments (Lidar). The schematic structure of the instrument is shown in Figure 1. My task in this project is to repair point cloud holes. In order to solve the above problems, combined with the hole repair method based on moving least square proposed by J.Wang, an improved point cloud hole repair algorithm is proposed by using the weighted orthogonal basis function to construct the shape function.

Extract point cloud hole boundaries
The extraction of the boundary of the point cloud hole is a key step in the repair of the point cloud hole. The quality of the hole boundary extraction directly affects the effect of the whole hole repair. In this paper, based on the triangular grid model, the hole extraction is realized by using the topological relationship between the point clouds between the triangular grids. the main steps are as follows: Step1: Building a Delaunay triangulation: Firstly, based on the half-edge data structure, the Delaunay triangulation is generated by the divide-and-conquer algorithm (Xie ,2012 andZhao ,2017).
Step2: Extracting hole boundaries based on Delaunay triangulation: Due to the lack of point data in the point cloud hole, the triangular network in this area is sparse, with fewer triangles and longer side lengths. According to this feature, a certain edge length threshold is set, and the boundary of the hole can be quickly identified by using the topological rules of the triangular mesh.
Step3: Obtaining k nearest neighbor points of hole boundary: In order to make the repaired point cloud smooth and continuous with the original terrain, after extracting the boundary of the point cloud hole, it is necessary to acquire the neighborhood point of the boundary point. This paper uses k-nearest neighbor search to obtain the neighborhood points of the boundary points.

Moving least squares
The Moving Least Squares (MLS) was systematically proposed by P. Lancaster and K. Salkauskas (1981) in the early 1980s, which mainly applied to curve and surface fitting. The MLS method is based on the supported weighting function and the polynomial basis function, and the fitting function suitable for the scatter model is established by the weighted least squares method. This method has greatly improved the least squares method, and the reconstructed curves and surfaces have better precision and smoothness.
The moving least squares method first needs to establish a fitting function, which is composed of the vector coefficient a(x) and its basis function P(x). Suppose U represents a subdomain on a two-dimensional fit region, and the expression of the fitting function is as follows: is basis function. The basis functions are linear, quadratic, ( ) and cubic, and so on, and the most common one is quadratic basis.
By the formula (1), the weighted least squares method is used to form the quadratic form as: Here, n is the number of affected nodes near the point x. is a smooth continuous weight function of the node with tight support properties. Within the tight support domain, >0, at its boundary and external , . is the node value at .

= 0 ( )
In this case, the coefficients that minimize the error are obtained by solving. Deriving the coefficient in formula ( ) where is the following matrix and is a matrix with ( ) diagonal equal to . (4) Substituting equation (5) into equation (1), f(x) can be expressed as:

The improved point cloud hole repair algorithm
The point cloud hole repair algorithm based on moving least squares can effectively repair the point cloud hole, and the repaired point cloud can be compared with the original terrain. However, when the basis function term m is large or the point cloud interval is too small, the matrix A (x) in formula (5) is prone to morbid or even singular, which makes it difficult to find the inverse matrix and thus difficult to obtain the correct solution.
Aiming at this problem, the weighted orthogonal polynomial is chosen as the basis function. The obtained matrix is neither ill nor singular, and the solution of the system of equations can be obtained without solving the matrix A(x) when solving the equations. Based on the improved moving least squares method, an improved hole repairing algorithm is proposed in this paper.

Improved moving least squares
For point set and weight , if a set of { } { }( = 1,2,⋯ ) functions satisfy the condition as 1 ( ), 2 ( ),⋯ ( ) Equation (6) (Cheng, 2003): equations are simplified to: In this way, when finding the coefficient , the inverse of the matrix is avoided, and the solution of the ill-conditioned equations is also avoided.

Point cloud hole repair algorithm based on improved moving least squares
Using the improved moving least squares method, the algorithm is mainly divided into the following steps: Step1: Calculate the feature surface of the hole point and its neighbors, assuming that the hole boundary point and its neighbor points are , ,…, , according to the principle of Here, barycentric coordinates satisfy the relationship: The construction matrix M satisfies the following relationship: Step2: Constructe the local coordinate system of the hole. The three coordinate axes set in the local coordinate system are U, V, S. And the center of gravity O is taken as the origin of the 1 hole coordinate system, and the feature surface is U V plane.

1
Project the boundary point and its neighbor point onto the feature plane. In the projection point, calculate the point farthest from the origin , and use the vector as the U-axis 1 direction, as the v-axis direction and the normal × are used as the S-axis of the coordinate system.
Step3: Establish implicit surfaces. based on improved moving least squares method. In this paper, we choose the quadratic basis as the basis function, and use the = (1, , , 2 , , 2 ) Schmidt orthogonalization method to construct the weighted orthogonality according to equation (8). The base function (1, 1 , the function expression for the fitted surface is: , 2 , 3 , 4 , 5 ) (11) 0 1 1 2 2 3 3 4 4 5 5 Z a a t a t a t a t a t       According to equations (3) and (9), we can solve the coefficient a(x) of the expression.
Step4: Calculate hole fill points. The hole boundary point and its neighbor point are projected into the hole local coordinate system, and then the rectangular bounding box of the hole boundary point on the UV plane is solved, and the rectangular bounding box on the uv plane is composed of , , min max min and . In order to effectively fill the holes, the newly added max sampling points are interpolated in the rectangular area, so the area needs to be meshed, and the grid points are the sampling points. The specific idea is to set a series of straight lines parallel to the U-axis and the V-axis, and set the spacing between the parallel lines to d (d is the average spacing of the original point clouds), then the coordinates of the grid points can be expressed as: The sampling points can be obtained by selecting the grid points in the hole, and then substitute the sample point coordinate into the fitting function (Equation 11) to find the corresponding s coordinate, and execute the process cyclically until the s coordinate values of all the sampling points are calculated. Finally the sampling point in the local coordinate system of the hole is changed to the original coordinate system, thus the sampling and filling of the hole area is completed.

EXPERIMENT ANALYSIS
In order to verify the effectiveness of the proposed method, We have implemented the hole-filling algorithm described using C++ and MATLAB. The experimental data of this paper is selected from the Airborne Lidar point cloud data on the Open Topography website. The website is a portal that provides highresolution terrain data and operating tools. It can download LiDAR data from the United States, Canada, Australia and other places for free.

Mountain area repair
The data one downloaded in this article is located in Marburg, New Zealand. The original point cloud data format is las format, and the point cloud density is 10.93pts/m2. After conversion to xyz format, a total of 4922397 data points are read, the minimum point cloud coordinate value is: (1683925.227579, 5362768.967441), and the maximum point cloud coordinate value is: (1685911.815391, 5364537.986929). The original point cloud is shown in Figure 2. The main features in the mountainous area are some bushes and some tall trees.

Figure 2. Original point cloud
In order to improve the efficiency of processing, this paper intercepts one of the areas to conduct experiments, a total of 346658 data points are read, as shown in Figure 3(a). Due to vegetation obstruction or terrain complexity, after filtering the original point cloud, 159135 non-ground points are filtered out and 187523 ground points are extracted. As shown in Figure 3(b), a point cloud hole appears. Using the method of this paper, the boundary extraction, knearest neighbor search, and hole repair are performed on the holes respectively. The results of the processing in each step are shown in Figure 4.
In order to visually see the effect of the repair, Figure 5 shows the side view before and after the hole repair.

Urban area repair
As shown in Figure 6, the data volume of airborne LiDAR point cloud in urban area is 122416, the data range is 160m × 255m, and the density of point cloud is 2.94 pts/m2. Due to the occlusion of ground object, a large number of holes appear in the original point cloud, of which hole 1, 2, 3, 4 is caused by the occlusion of buildings, hole 5 is caused by the occlusion of trees, and hole 6 may be caused by the absorption of water.

Comparison of various methods and accuracy analysis
In order to detect the repair accuracy of this method, the experiment also selects one of the regions to make a hole artificially, as shown in Figure 8. The point cloud data points that are artificially removed are taken as the known values, and the elevation values repaired by the algorithm are used as the predicted values, and the residual values and root mean square errors are calculated respectively. In addition, the inverse distance interpolation and nearest neighbor interpolation are used to repair, and compared with the repair accuracy of the algorithm. Using the proposed algorithm and two other different interpolation methods, the elevation values of the point cloud in the hole area are interpolated separately. The result of fitting the surface of the hole area is shown in Figure 9.
The predicted value is compared with the known value. And the residual value, the mean error and the root mean square errors are respectively calculated. The precision analysis is performed and the results are as follows:  Figure 9, it can be concluded that the effect of repairing the hole by using the inverse distance interpolation and the nearest neighbor interpolation is the worst.
The elevation values inserted in the two methods are prone to sudden changes. The interpolation and the terrain at the boundary of the hole is also quite consistent, but in the middle part of the hole, the interpolated point cloud surface is not continuous. The effect of repairing with this algorithm is the most ideal. The surface of the repaired point cloud is smooth and continuous, and the precision is high.

CONCLUSIONS
Airborne lidar can effectively obtain high-precision spatial data, which is a hot topic in current research. In the process of point cloud data scanning and measurement, due to the influence of the external environment, there will be certain errors and noises in the data, and even large blank areas will appear, which will have a negative impact on the subsequent generation of digital elevation models. In this paper, a hole repair algorithm based on moving least square is proposed. This method also improves the problem that there may be singularity or ill-condition in solving the matrix in the moving least square method. Without inversion of the matrix, the unknown coefficients of the implicit surface can be solved. The experimental results show that the proposed method can effectively repair point cloud holes with high accuracy, but there are still some problems in this method that need to be further studied, such as: for complex and small holes, the extracted boundary is not accurate enough. At the hole boundary, the accuracy of point cloud repair is not enough.