TERRAIN-ADAPTIVE GROUND FILTERING OF AIRBORNE LIDAR DATA BASED ON SALIENCY-AWARE THIN PLATE SPLINE

Ground filtering separates the ground and non-ground points from point clouds, which is the essential process for DEM generation, semantic segmentation, model reconstruction and so forth. Considering the topologically complex terrain environments, the segmentation results are prone to be disturbed dealing with steep slopes, buildings, bridges, cliffs, etc. from Airborne LiDAR point clouds. In this paper, a saliency-aware Thin-Plate-Spline (SATPS) interpolation method is proposed including two steps: saliency division and adaptive regularized TPS interpolation with relative variance coefficient. Firstly, the point clouds are indexed in 2D grids and segments are clustered step probing toward 8-adjacent scanning directions. Then, the saliency of each grid is calculated according to the elevation variance of adjacent segments towards each scanning direction. Subsequently, grids of high ground saliency are considered as candidates for seed point selection and then clustered by region growing. The TPS surface is interpolated for each cluster loosely fitting to the seed points involving an adaptive relative variance coefficient which is according to ground saliency and elevation deviation. And finally, the ground points are extracted around the TPS surface. Experimental results indicate that the proposed SATPS algorithm achieves better Type 1 accuracy and total accuracy than the state-of-the-art algorithms in scenes with complex terrain structures, which is practical to generate DEM products.


INTRODUCTION
Ground filtering is one of the most essential and foundational processes for urban mapping applications using point clouds, which plays an important role in DEM generation, model reconstruction, classification, object extraction, etc. Although relevant researches have been carried out for decades, it is obstinate to separate ground and non-ground points from LiDAR (Light Detection and Ranging) point clouds efficiently with high accuracy due to the uneven sampling density and the complicated terrain structures both in the urban scenes and mountainous areas.
Filtering methods for three major categories of LiDAR platforms vary due to the different viewing angle and sampling density (Roberts et al., 2019). Terrestrial Laser Scanning (TLS) and Mobile Laser Scanning (MLS) are ground-based platforms, static and moving respectively, the filtering strategies of which focus on the scan-line analysis and neighborhood analysis (Martí nez Sá nchez et al., Che et al., 2017). The capturing distance of Airborne Laser Scanning (ALS) is relatively consistent and the sampling density is much sparser. Considering limitation of scanning degree and flight elevation, ALS point clouds are commonly processed as 2.5D data for filtering and typical methods mainly include mathematical morphology (Li et al., 2013;Li et al., 2017;Meng et al., 2019), surface interpolation (Hu et al., 2014), segmentation (Chen et al., 2016), etc., which are discussed in this paper.
Since filtering is also a semantic segmentation task, machine learning methods classify the point clouds by trained classifiers * Corresponding author (Wang et al., 2015;Liu et al., 2019), including Information Vector Machine, Relevance Vector Machine, and Random Forest, etc. Deep learning-based methods sprung up as the neural networks towards point clouds processing became mature. Except for point clouds-based networks such as the PointNet (Qi et al., 2017), image-based methods (Yilmaz et al., 2018) are also applied since LiDAR along with multi-view images or photogrammetric point clouds from UAV photos (Gruszczyński et al., 2019;Zeybek and Şanlıoğlu, 2019) became feasible.
Surface interpolation-based methods extract the bare earth by fitting piece-wise successive surfaces for local patches. Triangular Irregular Networks (TIN) are reconstructed and iteratively densified to obtain the ground surface in classic pipelines (Zhao et al., 2016). Considering the uneven density of point clouds, Ma and Li (2019) resort to 3D alpha shapes based on ball-pivoting to generate 3D TINs. Curved surface-based methods perform better to preserve sharp details and conform terrain relief. Least-square, Thin-Plate-Spline (TPS) (Meng et al., 2019), and scanline spline (Martí nez Sá nchez et al., 2019) are commonly used to fit the surface while accurate interpolation points are strictly required for reasonable fitting. Improvements involve seed point selection such as moving-window weighted iterative least-squares fitting (Qin et al., 2017). Cheng et al. (2019) extended the TPS algorithm by improved skewness balancing. The state-of-the-art cloth simulation algorithm (Zhang et al., 2016) covers the surface inverted from points using a rigid cloth, which is easy-to-use and widely applied. However, the lack of explicit structures and the noise of point clouds make it a challenge to fit and accurately determine the segmentation surface for all categories of terrains.
To deal with topographically complex terrain environment, including steep slopes, multi-layer buildings, and cliffs, etc., researchers proposed solutions aimed at adaptive parameter settings for multitudinous terrains (Mongus and Žalik, 2012). Adaptive slope parameters (Susaki, 2012) are frequently optimized. Wan et al. (2018) resorted to the terrain relief index for tuning slope-related parameters. Considering density variance, scholars changed the iterative judgment criteria (Nie et al., 2017) to improve progressive TIN-based methods to decrease type I error (Dong et al., 2018;Zhao et al., 2016). Intrinsic structures of certain terrains are researched particularly and multi-scale filtering is conducted to preserve break lines (Yang et al., 2016). In addition, the ground and non-ground separation can be cast as a binary labeling problem and graph-cuts can be used to obtain a global optimized solution (He et al., 2018). On this basis, Hu et al. (2015) proposed the semi-global filtering method and balanced the energy function by adaptive ground saliency to adapt to discontinuous terrains.
In this paper, a terrain-adaptive ground filtering method is proposed based on saliency division and TPS interpolation. The workflow of the proposed SATPS method is shown in Figure 1. The point clouds are firstly indexed by 2D grids, which are clustered according to terrain saliency analysis, and then the regularized TPS is conducted towards clustered grids to determine the filtering surface. The contributions of this paper mainly consist of: 1. A saliency division method to select candidate ground grids and seed points for TPS interpolation; 2. An adaptive coefficient for the regularized item in TPS interpolation to loosen the seed points constraints.

SALIENCY DIVISION
The ground saliency (Hu et al., 2015) is an index to reveal the probability of a grid as ground according to the elevation varies for local point neighbors with a range of 0-1. Simply put, The larger the value, the higher the possibility that the point is a ground point. Considering the efficiency, the point clouds are rasterized in a regularized 2D grid to calculate the ground saliency.

Grid Segmentation
The point clouds are firstly indexed by a regularized grid in 2D and for each × grid, we assume the lowest elevation as the elevation of the grid.
Step probing is conducted to divide grids if distinct height differences exist. And adjacent grids with similar elevations are merged. The probing is conducted scanning towards 8 directions for each grid. For example, starting from the grid , along the scanning direction, assume the current grid is and the next grid is . If the elevation difference of and is larger than ℎ , or is empty, cluster grids from the starting grid to the current grid as one segment. And subsequently, starting from the next non-void grid to continue probing. Figure 2 illustrates a segmentation progress of the probing: along the scanning direction, points which are continuous in the horizontal direction and close in the elevation are divided into the same segment, then the grids are segmented into 4 segments caused by the elevation different higer than ℎ and the blank region. To ensure that each grid contains only one segment, grid is partitioned to 4x4 sub-grids and the same step probing process is conducted to detect sub-segments. If a grid contains ( ≥ 1) sub-segments, the grid is duplicated times and sub-segments are assigned to duplicated grids respectively. Mark each grid with the separation as: where is the k th duplicated grid of .

Ground Saliency Calculation
The saliency calculation by single scanning direction is based on the segments. And the saliency of each grid is the superposition of 8 scanning directions.
The ground saliency of the non-void grid is initialized as 1. Scanning towards 8 directions and along each scanning direction, if the difference between the end grid of the current segment and the start grid of the next segment is larger than ℎ , the saliency of grids in the current segment is decreased by 1. Traverse all the directions to obtain the aggregated saliency of each grid: if grids of a segment are higher than surrounding adjacent grids, the ground saliency will be 0 otherwise the value will be 1. The equation of ground saliency calculation is formed as: where ( , ) is saliency of . ( , ) represents the scanning directions to calculate . (⋅) represents the absolute elevation value (in meter) of the segment and is the segment index. In direction , belongs to segment | . and are starting grid and end grid of respectively.
[⋅] equals to 1 if the argument is true or 0 if false. Notably, as shown in Figure  3, among three downwards steps, the 2 nd is skipped for the calculation of 4 since 5 contains only 1 grid, and the 3 rd step condition is passed to 4 .

Figure 3. Ground saliency calculation
Saliency for particular terrains is optimized and discussed taking the building as an example. Grids of a building are generally higher than the adjacent grids whereas the buildings with sunken holes may break the rules. To repair the abnormal saliency, the morphological erosion process is conducted to grids whose saliency is higher than . Grids that eliminated by corrosion are marked as Scatter Ground, which will be excluded during the surface fitting but will be still considered as ground grids if they lie on the surrounding surface. Figure 4 shows the point clouds rendered by elevation (the strips in the pictures are caused by higher point clouds density) and the corresponding ground saliency. However, the ground saliency can only be the reference for ground points determination. Taken a building area shown in Figure 5 as example, the ground saliency of the sunkun part is higher but it is actually the inner part of the building. Morphological method can be utilized to revise the saliency to some extent but the substantive filtering issues cannot be resolved. The extraction of the ground surface will be discussed in the next section.
(a) Building structure with sunken parts (b) Ground saliency of a building area Figure 5. Saliency optimization for particular terrains

SURFACE INTERPOLATION
Since TPS is sensitive to the seed points, it is crucial to select reasonable groups of seed points, which should be interrupted by steps and steep slopes. In this section, the selection of seed points is based on the saliency calculated in the last section, and an optimized TPS interpolation is proposed to adapt to complex terrains.

Selection of Ground Seed Point
For grid (not marked as Scatter), if ( , ) > , is considered as the candidate grid for ground seed point selection. Candidate grids are clustered based on region growing which conforms the following conditions: where x, y, z = coordinates of the current grid; , , = coordinates of 8-adjacent grid; , = thresholds for elevation (in meter) and slope ratio.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) For grid that is split into several sub-grids, i.e., # ( , ) > 1, points in are clustered to connected adjacent grids respectively. In Figure 6, the red points have the same horizontal coordinates but with different elevation, and will be split into different sub-grids. This situation is common in multi-layer terrains, for instance, the bridges, roof eaves, cliffs, etc., while the traditional 2.5D ground filtering methods confuse to handle the 3D variance.

Regularized TPS Interpolation
Although grids are clustered for piece-wise TPS fitting, the grid height is roughly taken as the lowest elevation and the uneven sampling density leaves holes, which makes the seed points unreliable. Since the TPS fits seed points rigorously, to avoid over-fitting, an adaptive regularized item with coefficient λ is introduced according to the ground saliency. And ultimately, the ground points are extracted as points around the piece-wise TPS surface.
A seed point is the center of a candidate grid with the lowest elevation of the grid. TPS interpolation fits a smooth surface to interpolates the seed points. The TPS surface is expressed as: where is the distance, = √( − ) 2 + ( − ) 2 . ∑ =1 constrains the skewness of the surface and ensures the minimum curvature. 0 and the weight coefficients are unknown parameters, which can be solved by Equation 8 as given by Mongus and Žalik (2012): where is an N × N matrix and P is an N × 3 matrix, where and ̅ is the standard deviation of elevation and the mean saliency, respectively. Local variance is calculated involving the points in the current grid while the global variance is calculated involving all the points in the whole cluster. The influence of λ is illustrated in Figure 7. Figure 7. Surfaces using different λ fitting the same seed points. 1 < 2 < 3 and 1 = 0.

EXPERIMENTAL RESULTS
The software system for filtering was developed using Microsoft Visual Stdio 2015 under the Windows 7 x64 system with Intel core i5-4590 3.2Ghz CPU and 8G RAM. To verify the effectiveness and efficiency of the proposed SATPS algorithm, 10 test areas from 2 different airborne LiDAR datasets are selected and the characteristics of terrain environments are listed in Table 1. Area 1-9 come from Guangzhou dataset, each of which has 5 million points with a density of 10-20 pts/m 2 . Each area size is 520m×460m, and the elevation is between 38m to 117m, and the average height of buildings is around 24m. The 10th area is from Las Vegas dataset, provided by USGS 1 , which has more than 50 million points with a sparser density of 4-5 pts/ m 2 . The area size is 5380m×5350m, and the elevation is between 279.1m to 326.3m, while the average height of buildings is around 6m.

Area
Terrain characteristics

Quantitative results
The cross-matrices is used to evaluate the filtering results. Type I error the represents rejection of bare-earth points while Type II represents acceptance of points as bare-earth. The total error is the ratio of error classified points and the total points. Three state-of-the-art algorithms are selected as comparison including Progressive Morphological (PM) implemented by PCL 2 , Progressive TIN Densification (PTD) in the commercial software TerraSolid, and Cloth Simulation Filtering (CSF) (Zhang et al., 2016). The quantitative results are listed in Table 2, and the bold fonts represent the lowest value of each error type among all the algorithms. The ground truth of all the points of Guangzhou dataset is annotated manually using TerraSolid software. And for the Las Vegas dataset, the DEM of 1m resolution provided along with the LiDAR dataset is used as reference. The main parameters that affect the filtering results are the saliency threshold ℎ and . ℎ is used to calculate saliency and is set as 1m, which is reasonable for most urban scenes. is used in candidate determination.
is used for Scatter Ground identification, and it can be simply set as same as or slightly higher. The setting of is related to the topography situation of the survey area. In areas with dense and complex ground buildings, such as urban areas, a smaller value, 0.5, is used. In areas with sparse buildings or undulating terrain (rural or mountainous area), a larger value, 0.75, is used. The terrains of two experimental area are flat, thus the value is set to 0.5, and is empirically set to 0.75. However, even if non-ground grids are not predicted correctly at the stage of saliency division, it is still available to filter the non-grounds using TPS interpolation but the time cost may increase sharply. According to the saliency and elevation variance, the variance coefficient is smaller in smooth areas and larger with undulant topographic relief, which adapts the surface to complex terrains.

Qualitative results
The filtering results of site 1 and site 8 are displayed in Figure 8. In the original points, red represents buildings, green the vegetation, yellow the viaduct and pink for the ground. In the filtering results, yellow represents the non-ground points while pink for ground points. The terrain types in site 1 are relatively simple with smooth topological relief, which mainly includes buildings and vegetation. Site 9 is much more complex, which includes viaduct, slopes, buildings, and low vegetation. It can be seen that in most of the cases the filtering results are satisfying except for the low vegetation areas inside the circle. The boundaries of objects and grounds are confusing especially at the entrance of viaducts. Once the points grid is identified as nonground points in saliency division, it is unavailable to be considered as ground and not enough candidate grid can be found 2 PCL, https://www.pointclouds.org/ for the TPS fitting stage, which can be optimized by tuning the saliency thresholds. The DEM is generated from the ground points as shown in Figure  9. Compared to the other three algorithms, it is obvious that the proposed SATPS method performs better to filter building areas and preserve ground points inside the buildings at the same time , as shown in the red circle of Figure 9(d).

Buildings
Low vegetation Non-ground points

Ground points
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition)

Discussion
In most situations, PTD and the proposed SATPS algorithm achieve minimal Type I errors and the error rate is similar. The CSF performs well in flat areas but fails in a complex environment, the parameters are set specifically depending on the terrain structures, which cannot be uniformly set to the whole test area. In flat areas, the SATPS algorithm and CSF have minimal Type II errors while in complex scenes, the results of PTD and SATPS methods stand out. In building areas (Figure 10), PM and CSF consider the roofs as grounds. As for the total error rate, the proposed SATPS algorithm achieves the best results in more than half of the areas. The total error rate increases sharply when dealing with complicated terrain environments. The results of PM are acceptable in flat areas but in the complex scenes, parameters should be precisely tuned for satisfying results, which is not adaptive enough.
Generally, the proposed filtering algorithm performs well in different terrain environments. Especially in complex scenes with hybrid structures, such as multi-layer buildings and bridges, the boundary is preserved and the vegetation can be separated clearly at the same time. Besides, the TPS interpolation is conducted by clusters and reduces the time cost evidently, which can meet the requirement for practical production. As shown in the last line of Table 2, the mean value of each error type calcuated from all the areas indicates that, the proposed SATPS algorithm performs better for Type I errors and total errors among all the four methods.
The SATPS algorithm combines ground saliency calculation and TPS interpolation to filter ground points. Compared with PTD, PM, CSF algorithms, large scale and global operations are avoided, thus the calculation is faster. Tested with point clouds of different point density, different terrain conditions, the proposed algorithm performs good adaptability. When dealing with complex and mixed terrain environment, it can still ensure that the error rate does not increase significantly.

CONCLUSIONS
To provide a feasible and efficient scheme for ground filtering in point clouds with complex or mixed terrains, a saliency-aware TPS interpolation method is proposed in this paper. Ground saliency is calculated for each 2D grid and with a saliency division strategy, potential ground grids are considered as candidate grid for seed point selection. To extract the ground surface, piece-wise TPS interpolation is conducted with a regularized item according to the saliency, which relaxes the constraints that TPS interpolates all the seed points rigorously. The parameters used in the algorithm are mostly adaptive and only few of them need to be tuned manually. The proposed SATPS algorithm achieved better performance in complex terrain structures compared to traditional methods. The grid segmentation and clustering results affect the seed point selection and TPS interpolation process, which can be further improved by supervoxel-based methods under the constraint of structural features and will be researched in the future works.