A Fast Progressive TIN Densification Filtering Algorithm for Airborne LiDAR Data using Adjacent Surface Information

Point cloud filtering is a preliminary and essential step in various applications of airborne LiDAR (light detection and ranging) data, with progressive triangulated irregular network (TIN) densification (PTD) being one of the classic methods for filtering LiDAR point clouds. The PTD algorithm densifies ground points through iteration operation based on initial ground seed points. However, the poor performance in steeply sloped areas and time-consuming processing are serious drawbacks for PTD algorithms. In this paper, we propose a fast progressive TIN densification (FPTD) filtering algorithm for airborne LiDAR data using adjacent surface information. After carefully establishing parameters and removing outliers, our improved FPTD uses a sliding window to obtain significantly more initial ground seed points. And we modified some iterative determination criterion , including the definition of maximum relative elevation threshold and the introduction of signed computation, to eliminate avoidable non-ground points. Then adjacent surface information was utilized to iterate each point cloud block, which is the smallest unit that point cloud can be segmented. Additionally, the algorithm is easily run in a multi-threaded environment, further accelerating the filtering process to some extent. Experiments show that our proposed FPTD filtering algorithm is fast and robust. Compared to the PTD, the FPTD algorithm yields better error rates and kappa coefficients in 1/12 of the time required by the PTD.


I. INTRODUCTION
LiDAR (light detection and ranging) has become a popular remote sensing observation technique [1]. LiDAR systems fall into three main categories based on the data acquisition platform [2]: terrestrial laser scanning (TLS), mobile laser scanning (MLS), and airborne laser scanning (ALS). Of these, ALS use has grown significantly due to its ability to capture dense and accurate topographic data at high speed [3]. ALS has been widely employed in various fields [4], including digital terrain model (DTM) generation [5], [6], forest ecosystem investigation [7]- [9], 3D building modeling [10]- [12], and geological disaster surveys [13]- [15]. In most LiDAR applications, separating point clouds into ground and object points is a preliminary but essential step for subsequent processing [16], with various filtering methods proposed for automatically extracting bare earth surface points. These filters fall into seven broad categories according to the underlying technique of mathematical morphology, iterative interpolation, TIN (triangular irregular network), CS (cloth simulation), slope, segmentation, or machine learning.
Several comparisons of filtering algorithms [17]- [19] have shown that TIN-based filters, primarily those using PTD (progressive triangulated irregular network densification) algorithms, offer strong and stable abilities when acquiring bare ground points. Many researchers have improved the PTD algorithm in three large ways. The first improvement concerns the initial ground seed points, for it is widely accepted that the more initial ground seed points, the better filtering results. To identify more initial ground seed points, Zhang and Lin [20] improved the PTD algorithm using a point cloud segmentation method called segmentation using smoothness constraint (SUSC). Specifically, after selecting the lowest points in each grid cell as initial ground seed points, SUSC expands the set of ground seed points as much as possible. Similarly, Xu and Yue [21] firstly segmented point cloud data into several clusters based on plane fitting, then, followed by a coarse spatial clustering process using dual distances to obtain a set of initial ground seed points. Zhao et al. [22] first acquired potential ground seed points using a morphological method, and then used a translation plane fitting method to identify and remove non-ground points from the potential points. By assuming measurements that collide with cloth particles are ground points, the cloth simulation filtering [23] acquired initial ground seed points [4] and subsequently generated a high-quality initial provisional digital terrain model (DTM). These research results show that plentiful and accurate initial ground seed points influence the filtering results to a large extent.
The second PTD improvement concerns the iterative determination criterion. Lin and Zhang [24] modified the PTD algorithm by processing a segment unit rather than a single point. If the number of ground measurements is larger than the number of object measurements in a single segment, all of the points within the current segment are labeled as ground measurements. This segment-based filtering method preserves discontinuities in landscapes and removes the lower parts of large objects attached on the ground surface. To reduce the probability of incorrectly classifying non-ground points on the lower objects as ground points, only one point (the potential ground point with the shortest distance to the corresponding TIN facet) is selected as the ground point. The experimental results show that the revised PTD approach performs better than the PTD method in removing non-ground points, which are from low-to-the-ground objects [25]. And Shi et al. [26] proposed a parameter-free progressive TIN densification filtering algorithm based on the slope estimation.
The third improvement to the PTD is to shorten the processing time. Kang et al. [27] used state-of-the-art multi-core computing facilities to speed up the computationally intensive task by encapsulating the PTD filter into independent computing units with clearly faster results. However, this research ignores filtering accuracy. Other incremental improvements to the PTD were also made in [28], [29].
Most of the above studies point out that the performance of PTD algorithms is poor in steeply sloped areas. Because PTD algorithms rely on two thresholds (maximum angle and maximum distance [20]) and initial ground seed points to a large extent. Complex terrain will prevent PTD algorithms from extracting ground points through these fixed thresholds if there are not enough evenly-distributed initial ground seed points. Additionally, we find the PTD algorithm is very timeconsuming [27]for locating the corresponding triangle facet of a iterative point is time consuming, especially while the amount of ground points increase, causing a large number of triangle networks. However, most triangles are far away the iterative point, and can not help extract ground points but only increase time consumption. Therefore, we propose a fast progressive TIN densification filtering algorithm (FPTD) for airborne LiDAR data using adjacent surface information. In this method, we first use the sliding window algorithm to select large numbers of evenly-distributed initial ground seed points. Next, the point cloud divider breaks the data into a series of minimal processing blocks equal in size to the sliding window. Then, we iterate unclassified points in each block based on adjacent surface information, which greatly accelerates the filtering process. What's more, the improved algorithm is easily run in a multi-threaded environment, which can accelerate the TIN densification process still further. This method notably avoids reductions in filtering accuracy and actually increases it with most datasets. There are three main contributions of our paper: (1) A large number of evenlydistributed initial ground seed points can be readily obtained using the sliding window algorithm. (2) A new filtering strategy based on adjacent surface information is proposed, which can greatly shorten the processing time and improve the filtering accuracy. (3) A method to optimize the threshold selection of PTD filtering algorithms by classifying datasets according to the terrain complexity is proposed.
We organize our paper as follows. We present our improved PTD algorithm in detail in Section II. Section III presents the results and discussion. Finally, Section IV gives our conclusions.

II. METHOD
Our proposed method contains five key steps: (1) specifying parameters, (2) removing outliers, (3) obtaining initial ground seed points, (4) segmenting the point cloud into blocks, and (5) iterating all unclassified points and densifying the TIN. The detailed procedure of the proposed FPTD algorithm is illustrated in Fig.1. The following subsections provide full details of each step.

A. Specifying parameters
Our proposed algorithm requires seven parameters to be preset: (1) Radius threshold for outlier removal, r. This parameter eliminates outliers below the ground and retains sparse ground points in some areas. r should be appropriately big.
(2) Point number threshold for outlier removal, n. Similar to r, this parameter should be appropriately small.
(3) Maximum relative elevation threshold, e. Different from maximum distance [20], e is the relative elevation threshold from an unclassified point to the highest vertex in the corresponding triangle (seen in Fig. 2). If the elevation of an unclassified point minus that of the highest vertex is larger than e, the point will be labeled as the object point. Otherwise it will be labeled as a potential ground point.
(4) Maximum angle threshold, a. a is the angle threshold between the triangle plane and the line connecting an unclassified point to its corresponding closest triangle vertex [20] (seen in Fig. 2). If the angle is smaller than a, the point will be labeled as a potential ground point. Otherwise, it will be labeled as an object point.
(5) Sliding window size, w. w depends on the size of maximum artificial object (buildings, cars and so on). Specifically, to ensure the initial ground seed points are not obtained from building rooftops or other objects, w must be bigger than the maximum size of artificial objects in the point cloud dataset. . If w is smaller than the size of the biggest artificial objects in the point cloud, the lowest point in a sliding window may be from the objects, not the ground, causing false initial ground seed points. Therefore, to obtain more initial ground points, w should be as small as possible once it is greater than the biggest artificial object. Further, block size is always both equal to w for convenience in this paper.
(6) Sliding stride, s. s is a distance between two adjacent windows. It makes the window slide according to the stride distance.
(7) Neighborhood radius to the center point of a block point cloud, d. The value of d depends on the block size. It determines the initial ground seed points for a block point cloud, with d = 1.5 √ 2w providing a sufficient number in most cases. The sufficient initial ground seed points means that a block of data should contain the initial ground seed points along with those of its 8 neighborhood blocks.

B. Removing outliers
The PTD filtering algorithm performs well, relying on the hypothesis that there are no outliers with a lower height  than the ground. If there are, the radius outlier removal algorithm can be used to eliminate these low outliers. The radius outlier removal algorithm iterates every point in point cloud, find the points (outliers) whose amount of points in a radius threshold (r) is less than the point number threshold (n), and remove these points out. After using the radius outlier removal algorithm, it is a good idea to check the data again and remove the low outliers manually if any outliers remain.

C. Obtaining initial ground seed points
One of the most critical steps for PTD algorithms is obtaining accurate and appropriate initial ground seed points. We employ a new and simple strategy, the sliding window algorithm, to extract initial ground seed points ( Fig. 3) with the following steps. First, choose a square window of size w and set the sliding stride s. Second, calculate the minimum values of x and y directions to obtain the rectangular boundary of the point cloud, with the origin at the lower left corner. Third, slide the window from left to right within the rectangular boundary using the slide stride increment while finding the lowest point in each window. The point with the lowest elevation in the window is regarded as initial ground seed point. Fourth, slide the window up along the y direction while repeating the third step. Fifth, repeat the previous step until the window has covered the entire point cloud. Compared to most improved PTD algorithms [4], [20]- [22], the sliding window algorithm is much easier to implement, and more importantly, the obtained initial ground seed points are more evenly-distributed.
Note that (1) the coordinates of datasets have been projected before extracting initial ground seed points. (2) If the size of the last evaluated section of the point cloud in the x or y direction is smaller than the slide stride, that last section will be added to the prior window, meaning that the last window size may be bigger than others.

D. Segmenting point cloud into blocks
The PTD algorithm extracts ground points well, but the time it requires to do so is a serious drawback. To reduce the processing time, we segment unclassified point clouds into a series of blocks with overlapping regions (shown in Fig. 4) and use each block as a minimal processing unit. The minimal processing unit means the block can not be segmented again, which can accelerate the filtering process to the greatest extent. Concretely, the block size is equal to that of sliding window w, with the width of overlap equal to a multiple of the sliding stride. In data segmentation, if the width of the last part of the point cloud in the x or y direction is smaller than block size, it will be added to the previous block.
Next, the initial ground seed points (described in Section II-C) within a certain neighborhood distance (d) from the center of a block should be assigned to the corresponding block. These initial ground seed point are the initial block seed points in the next part. Notably, the value of d is extremely important, and its value must make sure the block to be processed is surrounded by adequate initial ground seed points.
, which means the initial ground seed points assigned to the block are from its corresponding 8 adjacent blocks. In this way, sufficient adjacent surface information can be utilized. To ensure that all of the points in each block will be located in the TIN constructed by the initial block seed points, simulated ground seed points (points highlighted in red in Fig. 4) are created for the block on the edge of the rectangle boundary before assigning initial ground seed points to blocks. The elevation of simulated ground seed points is the same as that of its nearest initial ground seed point. Afterwards, these simulated ground seed points are added to the initial ground seed point set, and then the new initial ground seed points are assigned to each block using the criteria above. Once a block and its corresponding initial block seed points are determined, we can employ multiple threads to process a large block point cloud.

E. Iterating unclassified points and densifying the TIN
In the PTD algorithm, all of the unclassified points are traversed based on the whole initial ground seed points in the densification process. However, each unclassified point only belongs to one corresponding triangle. Most triangles far from the iterative point have no effect on filtering results but consume a lot of processing time. And the more the ground points extracted, the more prominent this phenomenon becomes. Unlike the PTD, we only iterate the unclassified points in a block based on its corresponding initial block seed points. In this way, we preserve the adjacent surface information of points while reducing unnecessary processing.
The iteration process has five steps: (1) Within a single block, select the unclassified points having an elevation lower than the sum of the highest point in the initial block seed points and the maximum relative threshold. These are called the iterative unclassified points set, or the iterative set for short.  (2) and (3).
In this step, the two variables are slight different from those in the PTD. They are both signed. If their values are negative, it means the point is below the triangle facet and the point can also be identified as a ground point. In this way, it easily solves the problem mentioned in Section 2.3.1 of [25]. (4) Label the unclassified point. If the relative elevation and angle are both smaller than the corresponding threshold, the point will be labeled as a ground point. (5) Iterate all unclassified points in the iterative set and then update the elevation of the simulated ground seed points and the TIN according to the new ground points. Because the initial elevation of the simulated ground seed points are estimated from the corresponding nearest initial ground seed points, the estimated elevation may deviate far from the true adjacent ground surface with more and more ground points being extracted, resulting in that the ground points near the simulated ground seed points may not be extracted efficiently. Steps (1)-(5) are repeated until no unclassified points can be added to the ground point set.
By confining each iteration to a single block of the point cloud, multi-threaded techniques can be utilized to process multiple partitioned datasets simultaneously, shortening the processing time further to some extent.
In these equations, A, B, C, and D are the parameters of the plane equation of a facet. v is the relative elevation between an unclassified point (x u , y u , z u ) and the highest vertex in the corresponding triangle. θ is the angle between a triangle facet and the line connecting an unclassified point and its closest triangle vertex (x 0 , y 0 , z 0 ) in the corresponding triangle. z m is the elevation of the highest vertex in the corresponding facet. d is the signed distance from an unclassified point to its corresponding facet, meaning that the point is under the facet if d is less than zero, and above the facet otherwise. l is the distance between an unclassified point and its closest triangle vertex.

III. RESULTS AND DISCUSSION
To analyze the performance of our proposed algorithm, we used benchmark datasets provided by ISPRS Working Group III/3, including 8 unlabeled datasets from different terrains (4 urban sites and 4 rural sites) with 15 labeled samples from these sites. The eight datasets are named: CSite1, CSite2, CSite3, CSite4, FSite5, FSite6, FSite7, and FSite8. The 15 labeled samples are: s11, s12, s21, s22, s23, s24, s31, s41, s42, s51, s52, s53, s54, s61, and s71. Because these datasets contain different terrains, it is difficult to determine the optimal thresholds of PTD algorithms. Therefore, as shown in Table  I, we tentatively classified these data into four categories according to the terrain complexity (mainly the degree of slope change): flat (D1), flat with gentle slopes (D2), flat with scarp (D3), and steep (D4). For labeled samples, we calculated four indicators to evaluate the filtering results quantitatively as shown in Table II. We implemented our process on a desktop computer with an Intel Core i5-10400F CPU and 16 GB of RAM using Python under the Windows 10 operating system. We used eight threads for the iterative densification process.

A. Parameter analysis
All of the parameters (Section II-A) have different effects on the filtering results, but the maximum relative elevation threshold and maximum angle threshold matter the most. Thus, we paid particular attention to these two parameters.
Referring to other experiments [20], [24], [25], we varied the maximum angle threshold range from 2 to 42 degree in 4 degree increments and the maximum relative elevation threshold from 0.6 to 2.0 m in 0.2 m increments. To evaluate   the filtering results using different parameters, we chose subregions with fewer points as the experimental data for each terrain category: s54 of D1, s71 of D2, s24 of D3, and s52 of D4. We set other parameters based on our experience: r = 5, n = 10, w = 60, and s = 10.
For the first stage, we analyzed the maximum angle threshold while holding the maximum relative elevation threshold at 1.4. As seen in Fig. 5, although the trends in Fig. 5(c) and (d) are not particularly obvious, the kappa coefficient generally increased first and then decreased as the angle threshold increased, and the total error changed in the opposite direction. Additionally, the larger the maximum angle threshold, the smaller the type I error and the larger the type II error, indicating that more ground points were extracted correctly but more object points were misidentified at the same time. The optimal maximum angle thresholds were 18, 22, 30, and 38 for D1, D2, D3, and D4, respectively. From these thresholds, a conclusion can be drawn: the more complicated the terrain, the larger the optimal corresponding maximum angle threshold.
For the second stage, we examined the maximum relative elevation threshold in light of the optimal maximum angle threshold. The kappa coefficient in Fig. 6 fluctuated slightly less than that in Fig. 5, indicating that the FPTD filtering algorithm was more sensitive to the maximum angle threshold. Nonetheless, we were able to determine experimentally the optimal maximum relative elevation thresholds for different terrains: 1.0 for D1, 1.2 for D2, 0.8 for D3, and 1.4 for D4. Except for D3, the maximum relative elevation threshold increased with the terrain complexity. Both the optimal maximum angle and relative elevation thresholds of different terrains are listed in Table I. Our subsequent experiments were based on the these optimal thresholds.

B. Results
After specifying the optimal thresholds, we applied these values to all of the unlabeled datasets and obtained the filtering results shown in Table III. Although the number of extracted ground points did not show much useful information, the time cost was significantly reduced by the improved filtering algorithm. We have provided a number-time diagram (Fig. 7) to show the relationship between the total number of points and the processing time, using a linear function to fit these two variables. The processing time varied linearly with the number of points, indicating the high and robust processing efficiency.
High efficiency is usually accompanied by poor extraction results. Thus, to evaluate the performance of our algorithm, we selected the filtering results for FSite5, a relatively complex terrain, for visualization. Fig. 8 shows the overall high performance in our experiments. The buildings in the left area were eliminated robustly, and the vegetation in the lower right corner was also noticeably filtered out. In Fig. 8, (a-1) and (b-1) are vertical cross-sections of (a) and (b) in red highlight. In the same figure, (a-2) and (b-2) are horizontal cross-sections of (a) and (b), respectively, showing more specifically that the filtering results are consistent with the actual terrain.

C. Performance evaluation between the classic and improved PTD methods
We have worked to improve the PTD algorithm in three aspects: (1) the method of selecting initial ground seed points, (2) the iterative judgement criterion, and (3) the processing time. The filtering results in Section III-B indicate the improved PTD algorithm was effective and efficient. However, the differences in the results from the PTD and our improved PTD methods were obscure, so we tested the filtering performance of these two methods on labeled samples. The PTD was implemented according to [20]. Four different criteria were calculated, as given in Table IV. The data in the table shows that the improved PTD filter earned a lower T I score, but a  higher T II score than the PTD, indicating that our method extracted the ground points more reliably but generated more Total number: Total number of points in dataset (points) * These datasets were used to analyse the optimal parameters object points at the same time. Therefore, it is difficult to determine which filtering algorithm was better in terms of T I and T II. From the overall indicators T E and K, our improved algorithm performed better in most cases. More specifically, there were 12 datasets where the filtering results using the proposed algorithm were better than those using the PTD, while three datasets (s12, s22, and s71) showed slightly worse performance from our improved algorithm. The average total error of our FPTD was 6.12, while that of the PTD is nearly 2 times higher. The average kappa coefficient of the FPTD was almost 10% higher than that of the PTD. It is also worth mentioning the much shorter processing time required by our algorithm. For all of the labeled samples, the processing time is significantly reduced. In addition, the average processing time is nearly shortened by 12 times, and even nearly 30 times in some datasets (s42, s23).
Although the performance of our proposed algorithm was generally robust and efficient, there were some areas with poor filtering results. To explore the limitations in our method, we selected s12, s61, s24, and s53 for detailed analysis. The filtering results for these four regions were poorer than others, and each comes from a different type of terrain: s12 from D1, s61 from D2, s24 from D3, and s53 from D4. As shown in Figs. 9(b), 10(b), 11(b), and 12(b), the initial ground seed points obtained by two methods were markedly different. There were more initial ground seed points using our FPTD algorithm than the PTD. Most of the initial ground seed points acquired by the PTD were included as points acquired by the FPTD, apart from several simulated ground seed points. Thus, the FPTD readily obtained more dense initial ground seed points.
For s12 (Fig. 9), both the PTD and the FPTD performed terribly, especially in the white highlighted area. We classified this dataset into D1 because most areas of it are flat. However, there is a noticeable change in the topography (highlighted in white) that caused the optimal threshold to be unsuitable for the area. We tested this area using the optimal thresholds from D3, a more complex terrain, and obtained much better   filtering results. The scores of T I, T II, T E, and K using our proposed algorithm were 1.58%, 5.67%, 3.58%, and 92.83%, respectively. K was 6.83% higher than the original one. Notably, it is not that the larger the thresholds, the better the results. The scores of T I, T II, T E and K using the optimal thresholds from D4 were 0.61%, 12.44%, 6.38%, and 87.2%, indicating a drop in accuracy. For region s61 (Fig.  10), the same phenomenon occurred and the scores of T I, T II, T E, and K using the optimal thresholds of D4 were 0.96%, 29.02%, 1.93%, and 70.73%, respectively, indicating a better result. From these two datasets, we conclude that the datasets for mixed terrains should be classified based on the most complex area.
As shown in Fig. 11, most points in the right side (highlighted in white) of the dataset were filtered out using both the PTD and the FPTD, causing a high value of T I. There are two reasons for this. The first is that the dataset was too small, with a 121 m by 72 m rectangular boundary, resulting in incomplete continuity of the terrain and very few initial ground seed points. The block-wise operation magnified the problem. The second is that this dataset has noticeably abrupt changes in terrain. According to elevation, this dataset could be divided into three parts ( Fig. 11(a)): the lower left corner (area A), the top part (area B), and the right part (area C). The elevation of area C is the highest, but the initial ground seed points are concentrated in area A, making it hard to identify the ground points in areas B and C. However, the FPTD produced results that were still superior to the PTD's due to the larger number of initial ground seed points.
An example of datasets with steep terrain is shown in Fig.  12. To some extent, both the PTD and the FPTD algorithms failed to extract ground points in steep areas, and this drawback is common to most PTD-based algorithms. However, our proposed filtering performed much better than the PTD due to the larger number of initial ground seed points ( Fig. 12 (b)).
In addition, because the sum of the maximum terrain angle [20] plus maximum angle [20] was greater than 90 degrees, some object points with high elevation were identified as ground points using the PTD. The use of the maximum relative elevation threshold can help solve this problem.

IV. CONCLUSION
Due to the poor performance in steeply sloped areas and the time-consuming processing of classic PTD, we improved the PTD algorithm in three perspectives: the method of selecting initial ground seed points, the iterative judgement criterion, and the processing time. Experiments show our proposed FPTD filtering algorithm is fast and robust. Its average processing time was 1/12 of the PTD, and its average kappa coefficient of the FPTD was almost 10% higher than that of the PTD. Additionally, we draw three additional conclusions: (1) Using the sliding window algorithm, we can readily obtains a large number of evenly-distributed initial ground seed points. Decreasing the window size and stride size produces more initial ground seed points.
(2) Segmenting datasets into blocks and assigning the initial ground seed points in a certain neighborhood to blocks guarantees the surface spatial relationship between blocks and greatly shortens the processing time.
(3) The classification of datasets according to the terrain complexity is very helpful for determining the thresholds of PTD filtering algorithms.
Hongfu Li received the B.Eng. degree in space science and technology from the Chengdu University of Technology, China, in 2019, where he is currently pursuing the master's degree in Earth exploration and information technology. His main research interests comprise geo-hazard remote sensing applications, point cloud processing, hyperspectral image processing and deep learning.
Chengming Ye received the Ph.D. degree in Earth exploration and information technology from the Chengdu University of Technology, Chengdu, China, in 2011. He is currently an Associate Professor with the College of Geophysics, Chengdu University of Technology. He was a Visiting Scholar with the Department of Geography and Environmental Management, University of Waterloo, Canada, from 2016 to 2017. He has published more than 20 articles in peer-reviewed journals including the IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING (JSTARS). His main research interests include geo-hazard remote sensing applications, ecological remote sensing, and LiDAR data processing.
Zixuan Guo received the B.Eng. degree in space science and technology from the Chengdu University of Technology, China, in 2021, where she is currently pursuing the master's degree in Earth exploration and information technology. Her main research interests comprise geo-hazard remote sensing applications, and point cloud processing.
Ruilong Wei received the B.Eng. degree in Saptial information and digital technology from the Chengdu Univerisity of Technology, Chengdu, China, where he is currently working toward the MA.Eng. degree in Earth exploration and information technology. His research interests include Geo-hazard risk assessment, machine learning, and remote sensing applications.
Lixuan Wang received the B.Sc. degree in geoscience from Hubei University of Arts and Science, China, 2019. She is currently pursuing the master's degree in Earth exploration and information technology in Chengdu University of Technology. Her main research interests comprise Ecological and geohazard remote sensing application.