Filtering Airborne LiDAR Data Through Complementary Cloth Simulation and Progressive TIN Densiﬁcation Filters

: Separating point clouds into ground and non-ground points is a preliminary and essential step in various applications of airborne light detection and ranging (LiDAR) data, and many ﬁltering algorithms have been proposed to automatically ﬁlter ground points. Among them, the progressive triangulated irregular network (TIN) densiﬁcation ﬁltering (PTDF) algorithm is widely employed due to its robustness and e ﬀ ectiveness. However, the performance of this algorithm usually depends on the detailed initial terrain and the cautious tuning of parameters to cope with various terrains. Consequently, many approaches have been proposed to provide as much detailed initial terrain as possible. However, most of them require many user-deﬁned parameters. Moreover, these parameters are di ﬃ cult to determine for users. Recently, the cloth simulation ﬁltering (CSF) algorithm has gradually drawn attention because its parameters are few and easy-to-set. CSF can obtain a ﬁne initial terrain, which simultaneously provides a good foundation for parameter threshold estimation of progressive TIN densiﬁcation (PTD). However, it easily causes misclassiﬁcation when further reﬁning the initial terrain. To achieve the complementary advantages of CSF and PTDF, a novel ﬁltering algorithm that combines cloth simulation (CS) and PTD is proposed in this study. In the proposed algorithm, a high-quality initial provisional digital terrain model (DTM) is obtained by CS, and the parameter thresholds of PTD are estimated from the initial provisional DTM based on statistical analysis theory. Finally, PTD with adaptive parameter thresholds is used to reﬁne the initial provisional DTM. These contributions of the implementation details achieve accuracy enhancement and resilience to parameter tuning. The experimental results indicate that the proposed algorithm improves performance over their direct predecessors. Furthermore, compared with the publicized improved PTDF algorithms, our algorithm is not only superior in accuracy but also practicality. The fact that the proposed algorithm is of high accuracy and easy-to-use is desirable for users.

improvements have been proven successful, most algorithms require many user-defined parameters when densifying ground seed points. Moreover, these parameters are sensitive to particular scenes. Recently, a new surface-based method called cloth simulation filtering (CSF) algorithm has begun receiving attention because its user-defined parameters are few and easy to set [44,48,49]. CSF approximates the initial terrain with simulated cloth (CS) and then adds ground points from the unfiltered points to the initial terrain by using elevation information only, based on the assumption that the actual terrain should be close to a horizontal plane within a certain local area [17,34,39]. This algorithm, which only needs an integer parameter, can obtain an initial terrain that is near the actual terrain. However, unfiltered points are easily misclassified during the refinement of the initial terrain for steep mountainous areas, because the feature differences for ground and non-ground objects on steep surfaces are significantly different from flat surfaces. Sithole and Vosselman [2] suggested using additional contexts to solve this problem. In this sense, PTDF is one of the best implementations [45].
In addition, the challenge with PTDF is that the cautious tuning of parameters is necessary for various landforms. To overcome this limitation, several solutions for parameter threshold adaption have been proposed to improve the automation of PTDF. Axelsson [36] first proposed a statistical analysis-based method for adaptive determination of parameter thresholds. The method used characteristic values of statistics from ground points as the parameter thresholds. Experimental results show that this method can make PTDF achieve promising filtering performance in most types of landscape. However, the details of this method have not been reported. The possible reason is that the academic practitioner tends to maintain work proprietary [37]. Shi et al. [43] used kriging method to predict terrain slope and proposed a parameter threshold estimation formula by analyzing the connection between the terrain slope and the parameters of progressive TIN densification (PTD). However, the method may fail to provide satisfactory results in dense forest areas, because the kriging method cannot accurately predict the terrain slope when the ground seed point is sparse [5]. As mentioned above, estimated parameter thresholds are usually derived from the ground points. In addition, in order to obtain accurate thresholds, it is necessary to acquire as many ground points as possible before PTD.
To summarize, PTDF can acquire a filtering result with high precision, even in steep mountainous areas. The problem is that it usually relies on the detailed initial terrain and suitable parameter thresholds. CSF can obtain a high-quality initial terrain, which also provides a good foundation for parameter threshold estimation of PTD. However, misclassification can occur easily during the refinement of the initial terrain with respect to steep mountainous areas. To make use of the advantages of PTDF and CSF, this study proposes a novel algorithm to differentiate between ground and non-ground points from airborne LiDAR point clouds in a high-precision and easy-to-use manner by combining cloth simulation (CS) and progressive TIN densification (PTD). Recently, a new surface-based method called cloth simulation filtering (CSF) algorithm has begun receiving attention because its user-defined parameters are few and easy to set [44,48,49]. CSF approximates the initial terrain with simulated cloth (CS) and then adds ground points from the unfiltered points to the initial terrain by using elevation information only, based on the assumption that the actual terrain should be close to a horizontal plane within a certain local area [17,34,39]. This algorithm, which only needs an integer parameter, can obtain an initial terrain that is near the actual terrain. However, unfiltered points are easily misclassified during the refinement of the initial terrain for steep mountainous areas, because the feature differences for ground and non-ground objects on steep surfaces are significantly different from flat surfaces. Sithole and Vosselman [2] suggested using additional contexts to solve this problem. In this sense, PTDF is one of the best implementations [45].
In addition, the challenge with PTDF is that the cautious tuning of parameters is necessary for various landforms. To overcome this limitation, several solutions for parameter threshold adaption have been proposed to improve the automation of PTDF. Axelsson [36] first proposed a statistical analysis-based method for adaptive determination of parameter thresholds. The method used characteristic values of statistics from ground points as the parameter thresholds. Experimental results show that this method can make PTDF achieve promising filtering performance in most types of landscape. However, the details of this method have not been reported. The possible reason is that the academic practitioner tends to maintain work proprietary [37]. Shi et al. [43] used kriging method to predict terrain slope and proposed a parameter threshold estimation formula by analyzing the connection between the terrain slope and the parameters of progressive TIN densification (PTD). However, the method may fail to provide satisfactory results in dense forest areas, because the kriging method cannot accurately predict the terrain slope when the ground seed point is sparse [5]. As mentioned above, estimated parameter thresholds are usually derived from the ground points. In addition, in order to obtain accurate thresholds, it is necessary to acquire as many ground points as possible before PTD.
To summarize, PTDF can acquire a filtering result with high precision, even in steep mountainous areas. The problem is that it usually relies on the detailed initial terrain and suitable parameter thresholds. CSF can obtain a high-quality initial terrain, which also provides a good foundation Remote Sens. 2019, 11, 1037 4 of 23 for parameter threshold estimation of PTD. However, misclassification can occur easily during the refinement of the initial terrain with respect to steep mountainous areas. To make use of the advantages of PTDF and CSF, this study proposes a novel algorithm to differentiate between ground and non-ground points from airborne LiDAR point clouds in a high-precision and easy-to-use manner by combining cloth simulation (CS) and progressive TIN densification (PTD).

Test Data
Two airborne LiDAR point clouds namely ISPRS (International Society for Photogrammetry and Remote Sensing) data and dense point cloud were used to test the performance of the proposed algorithm. ISPRS data, which were collected with an Optech ALTM scanner, include eight sites (named as Site 1-8) consisting of different environments: four urban sites and four rural sites. Fifteen reference samples with different terrain features and land-cover types were selected from Site 1-7 to evaluate the filtering accuracy, as listed in Table 1. Site 8 was excluded due to lack of reference data. The point density is 0.4-1 points/m 2 for urban sites and 0.08-0.25 points/m 2 for rural sites, respectively. For each sample, the reference data was manually generated with prior knowledge and available aerial imagery [2]. Since the ISPRS data were collected almost two decades ago, the average point density was relatively low ranging from 0.08-1 points/m 2 . Recent advances in LiDAR technology allow continuous improvement of point density. To further test the performance of the proposed algorithm on high density of point cloud, data with an average point density of 31.12 points/m 2 were used. The data were obtained in 2018 using an unmanned aerial vehicle laser scanning system (SZT-R250) flying about 50 m above the ground level. The data were obtained over 35 ha, 10,912,743 points in total (including 1,749,516 ground points and 9,163,227 non-ground points), and elevations ranging from 0 m to 398.82 m, including the steep slopes, vegetation on hillside and discontinuities which are all typical features of the terrain in mountainous regions, as shown in Figure 2. We produced the reference data using a widely adopted method, which classifies point clouds into ground points and non-ground points by automatic classification tool in Terrasolid, and then the filtering results were further improved manually, through a repeated process in which the results were checked and rectified to ensure accuracy [50][51][52]. Although the classification results of the reference data are not completely accurate, the misclassification in reference data is much lower than that of the automatic filtering Remote Sens. 2019, 11, 1037 5 of 23 results [53]. Thus, in this study, we assume that the reference data is a correct classification of the LiDAR data.

Methods
A series of critical improvements are implemented to increase the accuracy and practicality of CSF and PTDF in this paper. First, we use CS to select the ground seed points instead of using the lowest points in a user-defined grid cell. This algorithm provides more ground seed points that are almost evenly distributed in general, thereby generating a high-quality initial provisional DTM. Second, parameter thresholds (i.e., the maximum angle (θ) and the maximal terrain slope (s)) for PTD are derived from the initial provisional DTM based on the statistical analysis theory to improve the automation of the filter. Third, the initial provisional DTM is refined by PTD with adaptive parameter thresholds. The entire workflow is composed of three parts, i.e., generation of initial provisional DTM, estimation of parameter thresholds and refinement of the initial provisional DTM, as shown in Figure 3.

Methods
A series of critical improvements are implemented to increase the accuracy and practicality of CSF and PTDF in this paper. First, we use CS to select the ground seed points instead of using the lowest points in a user-defined grid cell. This algorithm provides more ground seed points that are almost evenly distributed in general, thereby generating a high-quality initial provisional DTM. Second, parameter thresholds (i.e., the maximum angle (θ) and the maximal terrain slope (s)) for PTD are derived from the initial provisional DTM based on the statistical analysis theory to improve the automation of the filter. Third, the initial provisional DTM is refined by PTD with adaptive parameter thresholds. The entire workflow is composed of three parts, i.e., generation of initial provisional DTM, estimation of parameter thresholds and refinement of the initial provisional DTM, as shown in Figure 3.

Methods
A series of critical improvements are implemented to increase the accuracy and practicality of CSF and PTDF in this paper. First, we use CS to select the ground seed points instead of using the lowest points in a user-defined grid cell. This algorithm provides more ground seed points that are almost evenly distributed in general, thereby generating a high-quality initial provisional DTM. Second, parameter thresholds (i.e., the maximum angle (θ) and the maximal terrain slope (s)) for PTD are derived from the initial provisional DTM based on the statistical analysis theory to improve the automation of the filter. Third, the initial provisional DTM is refined by PTD with adaptive parameter thresholds. The entire workflow is composed of three parts, i.e., generation of initial provisional DTM, estimation of parameter thresholds and refinement of the initial provisional DTM, as shown in Figure 3.

Cloth Simulation
This method simulates the physical process of cloth-touching objects. After a cloth drops on an inverted (upside-down) landscape, the final shape of the cloth can be confirmed and regarded as an approximated terrain, as shown in Figure 4.

Cloth Simulation
This method simulates the physical process of cloth-touching objects. After a cloth drops on an inverted (upside-down) landscape, the final shape of the cloth can be confirmed and regarded as an approximated terrain, as shown in Figure 4.  A cloth is placed above the inverted measurements. (b) Each particle drops to the inverted measurements under the external force, and their displacement is computed based on the second law of Newton. If the cloth particles collide with the measurements, these particles will stop falling and be labeled unmovable particles. (c) Except for unmovable particles, each particle is moved according to the internal force produced by neighboring particles. Steps (b) and (c) are repeated until the maximum height variation of all particles is sufficiently small.
The process consists of an external force operation followed by an internal force operation. According to the second law of Newton, the relationship between position and forces is decided by where m represents the mass of the cloth particle, X is the position of a cloth particle at time t, F ext (X, t) is the external force, that is, gravity, and F int (X, t) is the internal force that is produced by interconnections between particles. When the cloth drops to the inverted measurement surface, the displacement of each cloth particle only under the external force is computed by Equation (2). The process consists of an external force operation followed by an internal force operation. According to the second law of Newton, the relationship between position and forces is decided by where m represents the mass of the cloth particle, X is the position of a cloth particle at time t, F ext (X, t) is the external force, that is, gravity, and F int (X, t) is the internal force that is produced by interconnections between particles. When the cloth drops to the inverted measurement surface, the displacement of each cloth particle only under the external force is computed by Equation (2). where ∆t represents the time step, and g is the acceleration. If the cloth particles collide with the inverted surface in the direction of its movement, these particles will stop moving and be labeled unmovable particles. After cloth particles are moved by an external force, an internal force is exerted to constrain the particles' displacement in the void areas of the inverted surface. The displacement (vector) of each particle is calculated as follows: where → d is the displacement vector of a particle, b is zero when the particle is unmovable, otherwise it equals to one, → p 0 represents the position of the current particle, → p i represents the position of the neighboring particle that connects with p 0 , and → n represents a normalized vector that points to the vertical direction, → n = (0, 0, 1) T . RT is the repeated time that controls the rigidness of the cloth. As shown in Figure 5, when RT is set to one, the movable particle will move 1/2 vertical distance (VD) between the two particles. When RT is set to two, the movable particle will move 3/4VD. When RT is set to three, the movable particle will move 7/8VD. When the terrain is extremely steep, RT is specified as a small value (RT = 1). When the terrain is relatively steep, RT is set to a medium value (RT = 2). When handling a flat terrain, a large value (RT = 3) must be set. where ∆t represents the time step, and g is the acceleration. If the cloth particles collide with the inverted surface in the direction of its movement, these particles will stop moving and be labeled unmovable particles. After cloth particles are moved by an external force, an internal force is exerted to constrain the particles' displacement in the void areas of the inverted surface. The displacement (vector) of each particle is calculated as follows: where d �⃗ is the displacement vector of a particle, b is zero when the particle is unmovable, otherwise it equals to one, p �⃗ 0 represents the position of the current particle, p �⃗ i represents the position of the neighboring particle that connects with p 0 , and n �⃗ represents a normalized vector that points to the vertical direction, n �⃗ = (0,0,1) T . RT is the repeated time that controls the rigidness of the cloth. As shown in Figure 5, when RT is set to one, the movable particle will move 1/2 vertical distance (VD) between the two particles. When RT is set to two, the movable particle will move 3/4VD. When RT is set to three, the movable particle will move 7/8VD. When the terrain is extremely steep, RT is specified as a small value (RT = 1). When the terrain is relatively steep, RT is set to a medium value (RT = 2). When handling a flat terrain, a large value (RT = 3) must be set.

Ground Seed Point Acquisition
When the maximum height variation of all cloth particles is sufficiently small, the cloth is regarded as stationary, and the CS (cloth simulation) process is terminated. The measurements that collide with cloth particles are labeled as ground seed points. The process is presented in Figure 6. The nearest measurement for each cloth particle in the X-Y plane is initially searched. If the heights of the cloth particle and its nearest measurement are equal, then this measurement is considered to collide with the corresponding cloth particle and labeled as a ground seed point. Figure 7 illustrates the labeled results by CS.

Ground Seed Point Acquisition
When the maximum height variation of all cloth particles is sufficiently small, the cloth is regarded as stationary, and the CS (cloth simulation) process is terminated. The measurements that collide with cloth particles are labeled as ground seed points. The process is presented in Figure 6. The nearest measurement for each cloth particle in the X-Y plane is initially searched. If the heights of the cloth particle and its nearest measurement are equal, then this measurement is considered to collide with the corresponding cloth particle and labeled as a ground seed point. Figure 7 illustrates the labeled results by CS. Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 23 Cloth particle Ground seed point Measurements Figure 6. Extraction of ground seed points by collision detection.

Initial Provisional DTM Construction
After ground seed points are acquired, the initial provisional DTM is constructed through Delaunay triangulating, as shown in Figure 8. Note that to ensure that each unfiltered point is located in this DTM, the four corners of the study areas are simulated as ground seed points whose heights are equal to the heights of their nearest ground seed points on the X-Y plane [37,38].

Parameter Threshold Estimation Based on Statistical Analysis
The performance of PTD greatly relies on three key parameters, including: the maximum angle (θ), the maximum distance (d) and the maximal terrain slope (s), as shown in Figure 9. The maximum angle (θ) is the maximum angle between the TIN facet and the line that links the candidate point to the facet's closest vertex. When the threshold of the maximum angle (θ) is set to

Initial Provisional DTM Construction
After ground seed points are acquired, the initial provisional DTM is constructed through Delaunay triangulating, as shown in Figure 8. Note that to ensure that each unfiltered point is located in this DTM, the four corners of the study areas are simulated as ground seed points whose heights are equal to the heights of their nearest ground seed points on the X-Y plane [37,38].

Parameter Threshold Estimation Based on Statistical Analysis
The performance of PTD greatly relies on three key parameters, including: the maximum angle (θ), the maximum distance (d) and the maximal terrain slope (s), as shown in Figure 9. The maximum angle (θ) is the maximum angle between the TIN facet and the line that links the candidate point to the facet's closest vertex. When the threshold of the maximum angle (θ) is set to

Initial Provisional DTM Construction
After ground seed points are acquired, the initial provisional DTM is constructed through Delaunay triangulating, as shown in Figure 8. Note that to ensure that each unfiltered point is located in this DTM, the four corners of the study areas are simulated as ground seed points whose heights are equal to the heights of their nearest ground seed points on the X-Y plane [37,38].
Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 23 Cloth particle Ground seed point Measurements Figure 6. Extraction of ground seed points by collision detection.

Initial Provisional DTM Construction
After ground seed points are acquired, the initial provisional DTM is constructed through Delaunay triangulating, as shown in Figure 8. Note that to ensure that each unfiltered point is located in this DTM, the four corners of the study areas are simulated as ground seed points whose heights are equal to the heights of their nearest ground seed points on the X-Y plane [37,38].

Parameter Threshold Estimation Based on Statistical Analysis
The performance of PTD greatly relies on three key parameters, including: the maximum angle (θ), the maximum distance (d) and the maximal terrain slope (s), as shown in Figure 9. The maximum angle (θ) is the maximum angle between the TIN facet and the line that links the candidate point to the facet's closest vertex. When the threshold of the maximum angle (θ) is set to

Parameter Threshold Estimation Based on Statistical Analysis
The performance of PTD greatly relies on three key parameters, including: the maximum angle (θ), the maximum distance (d) and the maximal terrain slope (s), as shown in Figure 9. The maximum angle (θ) is the maximum angle between the TIN facet and the line that links the candidate point to the facet's closest vertex. When the threshold of the maximum angle (θ) is set to a relatively large value, these non-ground points on the lower objects may be mistakenly classified as ground points. If the threshold of the maximum angle (θ) is small, many ground points related to key terrain characteristic may be mistakenly classified as non-ground points. The maximum distance (d) is the maximum distance from the candidate point to the corresponding TIN facet. The function of the maximum distance (d) is to guarantee that when the edge length of TIN facet is relatively large, the non-ground points on the low-large areas objects are not mistakenly classified as ground points [37,42]. The maximum terrain slope (s) is the maximum slope between the horizontal plane and the line that links the candidate point to the facet's lowest vertex. The value of the maximum terrain slope (s) decides whether to perform the mirrored operation, which makes the mirrored position of the candidate point is located symmetrically on the other side of the vertex in the corresponding TIN facet, relative to the original position, on the candidate point to handle the discontinuous surface [36][37][38]42]. More details are provided in Section 3.3.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 23 a relatively large value, these non-ground points on the lower objects may be mistakenly classified as ground points. If the threshold of the maximum angle (θ) is small, many ground points related to key terrain characteristic may be mistakenly classified as non-ground points. The maximum distance (d) is the maximum distance from the candidate point to the corresponding TIN facet. The function of the maximum distance (d) is to guarantee that when the edge length of TIN facet is relatively large, the non-ground points on the low-large areas objects are not mistakenly classified as ground points [37,42]. The maximum terrain slope (s) is the maximum slope between the horizontal plane and the line that links the candidate point to the facet's lowest vertex. The value of the maximum terrain slope (s) decides whether to perform the mirrored operation, which makes the mirrored position of the candidate point is located symmetrically on the other side of the vertex in the corresponding TIN facet, relative to the original position, on the candidate point to handle the discontinuous surface [36][37][38]42]. More details are provided in Section 3.3. In this study, the maximum angle (θ) is determined by statistical analysis, which is inspired by the work by Axelsson [36]. Statistic from initial provisional DTM is collected in the form of a discrete histogram of the terrain slopes, as shown in Figure 10. Finally, θ based on the median value, which is the value separating the higher half from the lower half of a data sample (a probability distribution), is determined from this histogram. For the maximum distance (d), it is set to the largest elevation difference in the study region. As shown by the example in Figure 11, when the horizontal distance between the non-ground point and the ground seed point is large, d should be set to a suitable threshold to prevent non-ground point from being incorrectly classified as ground point by the maximum angle threshold (θ). In this In this study, the maximum angle (θ) is determined by statistical analysis, which is inspired by the work by Axelsson [36]. Statistic from initial provisional DTM is collected in the form of a discrete histogram of the terrain slopes, as shown in Figure 10. Finally, θ based on the median value, which is the value separating the higher half from the lower half of a data sample (a probability distribution), is determined from this histogram.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 23 a relatively large value, these non-ground points on the lower objects may be mistakenly classified as ground points. If the threshold of the maximum angle (θ) is small, many ground points related to key terrain characteristic may be mistakenly classified as non-ground points. The maximum distance (d) is the maximum distance from the candidate point to the corresponding TIN facet. The function of the maximum distance (d) is to guarantee that when the edge length of TIN facet is relatively large, the non-ground points on the low-large areas objects are not mistakenly classified as ground points [37,42]. The maximum terrain slope (s) is the maximum slope between the horizontal plane and the line that links the candidate point to the facet's lowest vertex. The value of the maximum terrain slope (s) decides whether to perform the mirrored operation, which makes the mirrored position of the candidate point is located symmetrically on the other side of the vertex in the corresponding TIN facet, relative to the original position, on the candidate point to handle the discontinuous surface [36][37][38]42]. More details are provided in Section 3.3. In this study, the maximum angle (θ) is determined by statistical analysis, which is inspired by the work by Axelsson [36]. Statistic from initial provisional DTM is collected in the form of a discrete histogram of the terrain slopes, as shown in Figure 10. Finally, θ based on the median value, which is the value separating the higher half from the lower half of a data sample (a probability distribution), is determined from this histogram. For the maximum distance (d), it is set to the largest elevation difference in the study region. As shown by the example in Figure 11, when the horizontal distance between the non-ground point and the ground seed point is large, d should be set to a suitable threshold to prevent non-ground point from being incorrectly classified as ground point by the maximum angle threshold (θ). In this For the maximum distance (d), it is set to the largest elevation difference in the study region. As shown by the example in Figure 11, when the horizontal distance between the non-ground point and the ground seed point is large, d should be set to a suitable threshold to prevent non-ground point from being incorrectly classified as ground point by the maximum angle threshold (θ). In this study, dense ground seed points generate a relatively short edge length of TIN facet. Thus, d does not work and can be ignored. It usually is set to a rather large value, and we use the largest elevation difference in the study region as a threshold.  Figure 11. Simple demonstration of the maximum distance threshold setting.
In theory, the maximal terrain slope (s) should be the maximum slope derived from the final DTM. This threshold is usually not available prior to acquirement because ground points are not yet identified [20]. Fortunately, considering the initial provisional DTM is rather close to the actual terrain, we use this provisional DTM to obtain this threshold.

Refinement of Initial Provisional DTM Based on Progressive TIN Densification
After the initial provisional DTM and key parameters in PTD (progressive TIN densification) are determined, unfiltered points (including ground and non-ground points) are iteratively classified using PTD, and the detected ground points are used to update the initial provisional DTM. The procedure contains six core steps as follows.
Step1: Find the corresponding TIN facet of an unfiltered point.
Step2: Calculate the slope of this TIN facet. If the slope is greater than the maximal terrain slope (s), the unfiltered point will be replaced by the mirrored point in Step 3 to preserve discontinuous terrain. Then, go to Step4. Otherwise, the unfiltered point will be judged directly in Step 4. Step3: The mirrored operation is illustrated in Figure 12. Find the vertex P vertex (x v , y v , z v ) with the highest elevation value in the corresponding TIN facet of an unfiltered point P unfiltered (x u , y u , z u ). The mirrored point P mirror (x m , y m , z m ) of P unfiltered (x u , y u , z u ) is calculated as follows: Step4: Calculate the distance to the TIN facet and the angle to the vertex. If both computed values are lower than the maximum distance (d) and the maximum angle (θ), then the unfiltered point is labeled as a ground point. In theory, the maximal terrain slope (s) should be the maximum slope derived from the final DTM. This threshold is usually not available prior to acquirement because ground points are not yet identified [20]. Fortunately, considering the initial provisional DTM is rather close to the actual terrain, we use this provisional DTM to obtain this threshold.

Refinement of Initial Provisional DTM Based on Progressive TIN Densification
After the initial provisional DTM and key parameters in PTD (progressive TIN densification) are determined, unfiltered points (including ground and non-ground points) are iteratively classified using PTD, and the detected ground points are used to update the initial provisional DTM. The procedure contains six core steps as follows.
Step1: Find the corresponding TIN facet of an unfiltered point.
Step2: Calculate the slope of this TIN facet. If the slope is greater than the maximal terrain slope (s), the unfiltered point will be replaced by the mirrored point in Step 3 to preserve discontinuous terrain. Then, go to Step 4. Otherwise, the unfiltered point will be judged directly in Step 4.
Step3: The mirrored operation is illustrated in Figure 12. Find the vertex P vertex x v , y v , z v with the highest elevation value in the corresponding TIN facet of an unfiltered point P unfiltered x u , y u , z u .
The mirrored point P mirror x m , y m , z m of P unfiltered x u , y u , z u is calculated as follows:  Figure 11. Simple demonstration of the maximum distance threshold setting.
In theory, the maximal terrain slope (s) should be the maximum slope derived from the final DTM. This threshold is usually not available prior to acquirement because ground points are not yet identified [20]. Fortunately, considering the initial provisional DTM is rather close to the actual terrain, we use this provisional DTM to obtain this threshold.

Refinement of Initial Provisional DTM Based on Progressive TIN Densification
After the initial provisional DTM and key parameters in PTD (progressive TIN densification) are determined, unfiltered points (including ground and non-ground points) are iteratively classified using PTD, and the detected ground points are used to update the initial provisional DTM. The procedure contains six core steps as follows.
Step1: Find the corresponding TIN facet of an unfiltered point.
Step2: Calculate the slope of this TIN facet. If the slope is greater than the maximal terrain slope (s), the unfiltered point will be replaced by the mirrored point in Step 3 to preserve discontinuous terrain. Then, go to Step4. Otherwise, the unfiltered point will be judged directly in Step 4. Step3: The mirrored operation is illustrated in Figure 12. Find the vertex P vertex (x v , y v , z v ) with the highest elevation value in the corresponding TIN facet of an unfiltered point P unfiltered (x u , y u , z u ). The mirrored point P mirror (x m , y m , z m ) of P unfiltered (x u , y u , z u ) is calculated as follows: Step4: Calculate the distance to the TIN facet and the angle to the vertex. If both computed values are lower than the maximum distance (d) and the maximum angle (θ), then the unfiltered point is labeled as a ground point.
Step5: Calculate edge length of the TIN facet in the X-Y plane. If the ratio between the maximum edge length and the minimum edge length is smaller than the length ratio threshold (it is set to four), then the ground point is added to update the initial provisional DTM. Note that the introduction of length ratio threshold helps to avoid long narrow triangles, which is low credibility in terrain expression [40].
Step6: Repeat Steps 1-5 until all unfiltered points have been judged. Figure 13 illustrates the final DTM. Step5: Calculate edge length of the TIN facet in the X-Y plane. If the ratio between the maximum edge length and the minimum edge length is smaller than the length ratio threshold (it is set to four), then the ground point is added to update the initial provisional DTM. Note that the introduction of length ratio threshold helps to avoid long narrow triangles, which is low credibility in terrain expression [40].
Step6: Repeat Steps 1-5 until all unfiltered points have been judged. Figure 13 illustrates the final DTM.

Accuracy Indexes
Three accuracy indexes (i.e., type I, type II, and total errors) were used to assess the accuracy of the proposed algorithm quantitatively. The type I error is the percentage of ground points incorrectly classified as non-ground points. The type II error is the percentage of non-ground points incorrectly classified as ground points. The total error is the percentage of all mistakenly classified points. These indexes are calculated as follows [2]: where a is the number of reference ground points incorrectly classified as non-ground points, b is the number of reference non-ground points incorrectly classified as ground points, and c and d are the numbers of reference ground and non-ground points, respectively.

Testing with ISPRS Dataset
The distance between neighboring cloth particles, which is defined as cloth resolution (CR), has a strong relationship with the number of ground seed points. Generally, the smaller the CR value, the more ground seed points are obtained, and vice versa. Meanwhile, the parameter affects the time of cloth simulation because it decides how many cloth particles are generated and modified for the dataset. In CSF [44], CR is set to 0.5 m based on the principle of optimal accuracy. To objectively compare the proposed algorithm and CSF, the proposed algorithm under CR = 0.5 m was first tested. Then, the proposed algorithm under CR = 1 m was tested based on the balance between processing speed and filtering accuracy, which is discussed in Section 5.2. In addition, optimally tuned thresholds were used to compare our results with optimally tuned results. Specifically, we used a back-calibration optimization process based on the ground reference data provided by ISPRS to find the lowest total error of each sample. We tested the maximum angle threshold with 5° increments from 5° to 30°. The result with the lowest total error identifies the provisional optimized threshold (POT). POT was further refined with 1° increments from POT-2°

Accuracy Indexes
Three accuracy indexes (i.e., type I, type II, and total errors) were used to assess the accuracy of the proposed algorithm quantitatively. The type I error is the percentage of ground points incorrectly classified as non-ground points. The type II error is the percentage of non-ground points incorrectly classified as ground points. The total error is the percentage of all mistakenly classified points. These indexes are calculated as follows [2]: where a is the number of reference ground points incorrectly classified as non-ground points, b is the number of reference non-ground points incorrectly classified as ground points, and c and d are the numbers of reference ground and non-ground points, respectively.

Testing with ISPRS Dataset
The distance between neighboring cloth particles, which is defined as cloth resolution (CR), has a strong relationship with the number of ground seed points. Generally, the smaller the CR value, the more ground seed points are obtained, and vice versa. Meanwhile, the parameter affects the time of cloth simulation because it decides how many cloth particles are generated and modified for the dataset. In CSF [44], CR is set to 0.5 m based on the principle of optimal accuracy. To objectively compare the proposed algorithm and CSF, the proposed algorithm under CR = 0.5 m was first tested. Then, the proposed algorithm under CR = 1 m was tested based on the balance between processing speed and filtering accuracy, which is discussed in Section 5.2. In addition, optimally tuned thresholds were used to compare our results with optimally tuned results. Specifically, we used a back-calibration optimization process based on the ground reference data provided by ISPRS to find the lowest total error of each sample. We tested the maximum angle threshold with 5 • increments from 5 • to 30 • . The result with the lowest total error identifies the provisional optimized threshold (POT). POT was further refined with 1 • increments from POT-2 • to POT+2 • for obtaining the final optimized threshold (FOT). For example, if POT equals to 15 • , FOT will be selected from 13 • , 14 • , 15 • , 16 • and 17 • , and the result with the lowest total error identifies the FOT. Similarly, we tested the maximum distance threshold with 0.5 m increments from 0.5 m to 1.5 m. POT was further refined with 0.1 m increments from POT−0.2 m to POT + 0.2 m. We also tested the maximal terrain slope threshold with 10 • increments from 60 • to 90 • . POT was further refined with 1 • increments from POT-2 • to POT+2 • . Figure 14 shows to POT+2° for obtaining the final optimized threshold (FOT). For example, if POT equals to 15°, FOT will be selected from 13°, 14°, 15°, 16° and 17°, and the result with the lowest total error identifies the FOT. Similarly, we tested the maximum distance threshold with 0.5m increments from 0.5 m to 1.5 m. POT was further refined with 0.1 m increments from POT−0.2 m to POT + 0.2 m. We also tested the maximal terrain slope threshold with 10° increments from 60° to 90°. POT was further refined with 1° increments from POT-2° to POT+2°. Figure 14 shows   As shown in Figure 15, the average total error of the proposed algorithm under CR = 1 m is in the middle of all algorithms. However, the time consumption of the proposed algorithm under CR = 1 m is only about 0.25 times against CR = 0.5 m. We compared the proposed algorithm under CR = 1 m with all publicized improved PTDF algorithms, which completely used the 15 reference samples of the ISPRS dataset to assess their performances in terms of average total error. All parameters in the first three publicized improved PTDF algorithms were determined by the developers' experienced judgment for a particular terrain type, which is somewhat subjective. However, in the proposed algorithm, these parameters were automatically extracted from the initial provisional DTM. The detailed values of parameters in these algorithms are listed in Table 2. Note that the two algorithms proposed by Zhang and Lin [37], Lin and Zhang [38] share the same parameter values. 25  As shown in Figure 15, the average total error of the proposed algorithm under CR = 1 m is in the middle of all algorithms. However, the time consumption of the proposed algorithm under CR = 1 m is only about 0.25 times against CR = 0.5 m. We compared the proposed algorithm under CR = 1 m with all publicized improved PTDF algorithms, which completely used the 15 reference samples of the ISPRS dataset to assess their performances in terms of average total error. All parameters in the first three publicized improved PTDF algorithms were determined by the developers' experienced judgment for a particular terrain type, which is somewhat subjective. However, in the proposed algorithm, these parameters were automatically extracted from the initial provisional DTM. The detailed values of parameters in these algorithms are listed in Table 2. Note that the two algorithms proposed by Zhang and Lin [37], Lin and Zhang [38] share the same parameter values. We also tested the maximal terrain slope threshold with 10° increments from 60° to 90°. POT was further refined with 1° increments from POT-2° to POT+2°. Figure 14 shows       Table 2. The values of parameters in all filters for each sample. Note that "c", "θ", "d" and "s" represent "grid cell size", "maximum angle", "maximum distance" and "maximal terrain slope", respectively.  Table 3a shows the total error of the proposed algorithm against four other improved PTDF algorithms. The average total error of the proposed algorithm is 6.95%. Compared to the results of other improved PTDF algorithms, the proposed algorithm achieves the highest precision for 7/15 samples in terms of the total error and is also the best for average values. Seven samples achieving the highest precision are S22, S23, S52, S53, S54, S61, and S71. These samples contain a variety of complex terrain features, including steep slopes, terraced terrain, and discontinuity. The statistical results in Table 3b suggest the average type I and II errors of the proposed algorithm are 4.6% and 11.42%. Compared to the results of other improved PTDF algorithms, the proposed algorithm achieves the highest precision for 11/15 samples in terms of the type I error and is also the best for average values. However, the proposed algorithm is poor in avoiding the type II error. The main reason is more ground seed points will increase the risk of the type II error, while they are controlled to an acceptable level, because the accuracy of ground seed points obtained from CS is extremely high (the average overall precision is 98.39%). Moreover, the inclination to type II error may not be a fatal flaw for the filter strategy, taking into consideration that type II error can be more easily handled by human editing than type I error [2]. Table 3a. Performance comparison between the proposed algorithm and the publicized improved PTDF algorithms in terms of the total error (%). The best results across all the algorithms are marked in red.

Samples
Zhang and Lin Lin Figure 16 shows our results and optimally tuned results under CR = 0.5 m and CR = 1 m in terms of average total error. We can observe two characteristics. First, it appears that our results are relatively close to the optimally tuned results. The fact that the proposed algorithm is effective to parameter estimation demonstrates the wide applicability of this algorithm. Second, the difference between our result and optimally tuned result under CR = 0.5 m is lower than the difference under CR = 1 m, and the difference between optimally tuned results under CR = 0.5 m and CR = 1 m is small. It indicates that the detail degree of the initial terrain and the accuracy of parameter estimation have a significant positive correlation relationship.  I  II  I  II  I  II  I  II  I  II  S11 25      Figure 17 shows the DTMs and the space distribution of type I and II errors (SDEs) generated by the proposed algorithm under CR = 0.5 m and CR = 1 m in several representative samples, which are named as DTM CR=0.5 , DTM CR=1 , SDE CR=0.5 and SDE CR=1 . S11 has a considerable slope change, and many mixtures of buildings and vegetation are on the hillside. The visualizations of DTM CR=0.5 and DTM CR=1 exhibit similar appearances. Exploration of further details for some recognizable differences between them shows DTM CR=0.5 is closer to the reference DTM than DTM CR=1 in the terrain discontinuities, as shown in the rectangle regions in Figure 17b-d. However, compared to DTM CR=1 , there are more non-ground objects in DTM CR=0.5 , as shown in the circle regions in Figure 17b-d. Compared to SDE CR=1 , SDE CR=0.5 shows lower type I error but higher type II error. The reason can be explained by exploring the procedure of CS. When cloth gradually approaches terrain with more cloth particles, a more detailed initial terrain is constructed, thereby assisting PTD to correctly classify the ground point in complex terrain. However, the ground seed points are not 100% correctly obtained by CS. Thus more seed points will increase the risk of type II error. For S52, the ground surface has many abrupt changes. In comparison with the reference DTM, DTM CR=0.5 and DTM CR=1 can effectively preserve the steep and terraced slopes. and DTM CR=1 exhibit similar appearances. Exploration of further details for some recognizable differences between them shows DTM CR=0.5 is closer to the reference DTM than DTM CR=1 in the terrain discontinuities, as shown in the rectangle regions in Figure 17b-d. However, compared to DTM CR=1 , there are more non-ground objects in DTM CR=0.5 , as shown in the circle regions in Figure  17b-d. Compared to SDE CR=1 , SDE CR=0.5 shows lower type I error but higher type II error. The reason can be explained by exploring the procedure of CS. When cloth gradually approaches terrain with more cloth particles, a more detailed initial terrain is constructed, thereby assisting PTD to correctly classify the ground point in complex terrain. However, the ground seed points are not 100% correctly obtained by CS. Thus more seed points will increase the risk of type II error. For S52, the ground surface has many abrupt changes. In comparison with the reference DTM, DTM CR=0.5 and DTM CR=1 can effectively preserve the steep and terraced slopes.
S11 (a) S11 (b) S11 (c) S11 (d) S11 (e) S11 (f) Ground Non-ground Type I error Type II error Figure 17. Results of each group (S11 and S52 are selected as representatives). First column: DSMs, second column: DTMs generated from the references, third column: DTMs produced by the proposed algorithm under CR = 0.5 m, fourth column: DTMs produced by the proposed algorithm under CR = 1 m, fifth column: spatial distributions of the type I and type II errors produced by the proposed algorithm under CR = 0.5 m, sixth column: spatial distributions of the type I and type II errors produced by the proposed algorithm under CR = 1 m.

Testing with Dense Point Cloud
To quantitatively assess the advantages of the proposed algorithm under CR = 1 m, which achieves the complementary advantages of CSF and PTDF, we compared the three types of errors of the proposed algorithm with CSF and PTDF, where CR value in CSF was set to 1 m. The results are shown in Figure 18. We can observe three characteristics. First, the type I error is obviously larger than the type II error for PTDF, CSF, and the proposed algorithm. The main reason is that the number of non-ground points in this data is more than the ground points [34]. Specifically, non-ground points comprise 83.97% of all the points. Thus, only a few ground points being misclassified as non-ground points will yield a large type I error. Second, in terms of the type I and total errors, the proposed algorithm is better than PTDF. It occurs because CS is capable of providing enough ground seed points evenly distributed in general, which can cover terrain characteristics and retains finer details, especially for hilltops and steep slopes. The initial provisional DTM based on the ground seed points is close to the actual terrain, which is important because the subsequent PTD is based on the initial provisional DTM [40]. Third, compared with CSF, the proposed algorithm performs better in terms of three types of errors. It indicates that, when refining the initial Figure 17. Results of each group (S11 and S52 are selected as representatives). First column: DSMs, second column: DTMs generated from the references, third column: DTMs produced by the proposed algorithm under CR = 0.5 m, fourth column: DTMs produced by the proposed algorithm under CR = 1 m, fifth column: spatial distributions of the type I and type II errors produced by the proposed algorithm under CR = 0.5 m, sixth column: spatial distributions of the type I and type II errors produced by the proposed algorithm under CR = 1 m.

Testing with Dense Point Cloud
To quantitatively assess the advantages of the proposed algorithm under CR = 1 m, which achieves the complementary advantages of CSF and PTDF, we compared the three types of errors of the proposed algorithm with CSF and PTDF, where CR value in CSF was set to 1 m. The results are shown in Figure 18. We can observe three characteristics. First, the type I error is obviously larger than the type II error for PTDF, CSF, and the proposed algorithm. The main reason is that the number of non-ground points in this data is more than the ground points [34]. Specifically, non-ground points comprise 83.97% of all the points. Thus, only a few ground points being misclassified as non-ground points will yield a large type I error. Second, in terms of the type I and total errors, the proposed algorithm is better than PTDF. It occurs because CS is capable of providing enough ground seed points evenly distributed in general, which can cover terrain characteristics and retains finer details, especially for hilltops and steep slopes. The initial provisional DTM based on the ground seed points is close to the actual terrain, which is important because the subsequent PTD is based on the initial provisional DTM [40]. Third, compared with CSF, the proposed algorithm performs better in terms of three types of errors. It indicates that, when refining the initial provisional DTM, PTD outperforms the method that uses elevation information only. This finding is also shared by previous studies, which have illustrated that additional contexts can improve filtering results [2,47].
The generated DTM and the cross-section of point clouds for quality evaluation are shown in Figure 19. In general, the DTM generated from the proposed algorithm shows the closer appearances with the reference DTM, i.e., most of non-ground points are removed and ground points remain. Details of the comparison show that the proposed algorithm can preserve the ground points on hill tops and steep slopes, whereas PTDF and CSF hardly preserve terrain characteristics in areas with large terrain variations.
Remote Sens. 2018, 10, x FOR PEER REVIEW  16 of 23 provisional DTM, PTD outperforms the method that uses elevation information only. This finding is also shared by previous studies, which have illustrated that additional contexts can improve filtering results [2,47]. The generated DTM and the cross-section of point clouds for quality evaluation are shown in Figure 19. In general, the DTM generated from the proposed algorithm shows the closer appearances with the reference DTM, i.e., most of non-ground points are removed and ground points remain. Details of the comparison show that the proposed algorithm can preserve the ground points on hill tops and steep slopes, whereas PTDF and CSF hardly preserve terrain characteristics in areas with large terrain variations. Figure 18. Three types of errors on data with high point density for PTDF, CSF, and our algorithm. An awareness of the relationship between filtering algorithms and point density helps to balance the processing speed and filtering accuracy [53]. To analyze the effect of point density on the proposed algorithm, the low-density versions of original LiDAR data were generated using a simple point thinning algorithm, keeping only every nth (for example, the 10th) point along the time line from the original LiDAR data [62]. In this study, we used every 2nd, 5th, 10th and 15th points from original LiDAR data and generated a series of low-density versions with average point density of 15.56, 6.22, 3.11and 2.07 points/m 2 (named as LD 2nd, LD 5th, LD 10th and LD 15th). As shown in Figure 20, it can be seen that point density has relatively less effect on the filtering results. Specifically, the root mean square errors (RMSEs) of the type I, type II and total errors are 2.02%, provisional DTM, PTD outperforms the method that uses elevation information only. This finding is also shared by previous studies, which have illustrated that additional contexts can improve filtering results [2,47]. The generated DTM and the cross-section of point clouds for quality evaluation are shown in Figure 19. In general, the DTM generated from the proposed algorithm shows the closer appearances with the reference DTM, i.e., most of non-ground points are removed and ground points remain. Details of the comparison show that the proposed algorithm can preserve the ground points on hill tops and steep slopes, whereas PTDF and CSF hardly preserve terrain characteristics in areas with large terrain variations. Figure 18. Three types of errors on data with high point density for PTDF, CSF, and our algorithm. An awareness of the relationship between filtering algorithms and point density helps to balance the processing speed and filtering accuracy [53]. To analyze the effect of point density on the proposed algorithm, the low-density versions of original LiDAR data were generated using a simple point thinning algorithm, keeping only every nth (for example, the 10th) point along the time line from the original LiDAR data [62]. In this study, we used every 2nd, 5th, 10th and 15th points from original LiDAR data and generated a series of low-density versions with average point density of 15.56, 6.22, 3.11and 2.07 points/m 2 (named as LD 2nd, LD 5th, LD 10th and LD 15th). As shown in Figure 20, it can be seen that point density has relatively less effect on the filtering results. Specifically, the root mean square errors (RMSEs) of the type I, type II and total errors are 2.02%, An awareness of the relationship between filtering algorithms and point density helps to balance the processing speed and filtering accuracy [53]. To analyze the effect of point density on the proposed algorithm, the low-density versions of original LiDAR data were generated using a simple point thinning algorithm, keeping only every nth (for example, the 10th) point along the time line from the original LiDAR data [62]. In this study, we used every 2nd, 5th, 10th and 15th points from original LiDAR data and generated a series of low-density versions with average point density of 15.56, 6.22, 3.11and 2.07 points/m 2 (named as LD 2nd, LD 5th, LD 10th and LD 15th). As shown in Figure 20, it can be seen that point density has relatively less effect on the filtering results. Specifically, the root mean square errors (RMSEs) of the type I, type II and total errors are 2.02%, 0.81%, and 0.45%, respectively. Additionally, there is a general tendency for the type II error to be slightly negative. The possible reason is that the effect of the type II error from high-density data is limited because non-ground points are surrounded closely by ground points [63]. In terms of the type I error, the proposed algorithm does not show a recognizable tendency.
Remote Sens. 2018, 10, x FOR PEER REVIEW 17 of 23 0.81%, and 0.45%, respectively. Additionally, there is a general tendency for the type II error to be slightly negative. The possible reason is that the effect of the type II error from high-density data is limited because non-ground points are surrounded closely by ground points [63]. In terms of the type I error, the proposed algorithm does not show a recognizable tendency.

Accuracy of Ground Seed Points
To explore the effects of ground seed point extraction, the numbers of ground seed points for 15 samples obtained by PTDF and the proposed algorithm under CR = 1 m are listed in Table 4. The results show that, on average, approximately 370 ground seed points for 15 samples are extracted by PTDF. However, the proposed algorithm can provide approximately 10368 ground seed points, which indicates that it can significantly increase the number of ground seed points. In addition, overall precision (OP) is employed to evaluate the quality of ground seed points. OP is the rate of correctly classified ground points in all extracted ground seed points and is calculated as: where a is the number of ground points correctly classified, and c is the number of non-ground points incorrectly classified as ground points. The statistics in Table 4 indicate that the accuracy of ground seed points acquired by the proposed algorithm is extremely high, and the average OP is 98.39%. To assess the quality of ground seed point extraction intuitively, S11 is selected as the representative case for demonstration. As shown in Figure 21a,b, the ground seed points provided by the proposed algorithm can cover the terrain features and retain fine details, especially for steep slopes. The detailed changes in the terrain can be presented relatively well, especially for the regions of significant variations in the terrain. Figure 21d shows that the initial provisional DTM constructed from the proposed algorithm is rather near the actual terrain, which is important because the subsequent PTD is based on the initial terrain.

Accuracy of Ground Seed Points
To explore the effects of ground seed point extraction, the numbers of ground seed points for 15 samples obtained by PTDF and the proposed algorithm under CR = 1 m are listed in Table 4. The results show that, on average, approximately 370 ground seed points for 15 samples are extracted by PTDF. However, the proposed algorithm can provide approximately 10368 ground seed points, which indicates that it can significantly increase the number of ground seed points. In addition, overall precision (OP) is employed to evaluate the quality of ground seed points. OP is the rate of correctly classified ground points in all extracted ground seed points and is calculated as: where a is the number of ground points correctly classified, and c is the number of non-ground points incorrectly classified as ground points. The statistics in Table 4 indicate that the accuracy of ground seed points acquired by the proposed algorithm is extremely high, and the average OP is 98.39%. To assess the quality of ground seed point extraction intuitively, S11 is selected as the representative case for demonstration. As shown in Figure 21a,b, the ground seed points provided by the proposed algorithm can cover the terrain features and retain fine details, especially for steep slopes. The detailed changes in the terrain can be presented relatively well, especially for the regions of significant variations in the terrain. Figure 21d shows that the initial provisional DTM constructed from the proposed algorithm is rather near the actual terrain, which is important because the subsequent PTD is based on the initial terrain.  Figure 21. Results of representative sample: (a) Overlay result of reference DTM and ground seed points (red points) produced by PTDF, (b) overlay result of reference DTM and ground seed points (red points) produced by the proposed algorithm, (c) initial provisional DTM generated from ground seed points of (a), (d) initial provisional DTM generated from ground seed points of (b).

Parameter Analysis
For ground seed points of the 15 reference samples, the simulation time and the total error at different cloth resolution (CR) values (ranging from 0.3 m to 1.5 m) are shown in Figure 22. The results show that a small CR value usually has large time consumption, and the time cost is relatively stable after the CR value exceeds 0.8 m. The accuracies are relatively stable, and the average total error is less than 10% when CR ranges from 0.3 m to 1.1 m, and the accuracies are the highest around 0.6 m. In view of the balance of efficiency and accuracy, the value of this parameter is usually set from 0.8 m to 1.1 m. In terms of the principle of optimal accuracy, the value of this parameter is usually set from 0.5 m to 0.7m. Total error Figure 21. Results of representative sample: (a) Overlay result of reference DTM and ground seed points (red points) produced by PTDF, (b) overlay result of reference DTM and ground seed points (red points) produced by the proposed algorithm, (c) initial provisional DTM generated from ground seed points of (a), (d) initial provisional DTM generated from ground seed points of (b).

Parameter Analysis
For ground seed points of the 15 reference samples, the simulation time and the total error at different cloth resolution (CR) values (ranging from 0.3 m to 1.5 m) are shown in Figure 22. The results show that a small CR value usually has large time consumption, and the time cost is relatively stable after the CR value exceeds 0.8 m. The accuracies are relatively stable, and the average total error is less than 10% when CR ranges from 0.3 m to 1.1 m, and the accuracies are the highest around 0.6 m. In view of the balance of efficiency and accuracy, the value of this parameter is usually set from 0.8 m to 1.1 m. In terms of the principle of optimal accuracy, the value of this parameter is usually set from 0.5 m to 0.7 m.
After ground seed points are obtained, the remaining ground points are extracted using PTD. To analyze the sensitivity of the parameters (θ and d) on the improvement of filtering accuracy, S11 was selected as a representative in this section. We adjusted the two parameters as follows: θ ∈ 10  Figure 23. It shows that the total error of the proposed algorithm is lower than PTDF under all parameter configurations, which indicates obtaining more ground seed points is an important way to improve the filtering accuracy of PTDF. This is in agreement with the previous studies that almost all improved PTDF focus on the densification of ground seed points. For PTDF and the proposed algorithm, the RMSE of the total error is 4.78% and 1.51% with the change of θ and d, respectively.
It illustrates that the proposed algorithm is insensitive to the parameter selection even though the algorithm is parameter-dependent.

Parameter Analysis
For ground seed points of the 15 reference samples, the simulation time and the total error at different cloth resolution (CR) values (ranging from 0.3 m to 1.5 m) are shown in Figure 22. The results show that a small CR value usually has large time consumption, and the time cost is relatively stable after the CR value exceeds 0.8 m. The accuracies are relatively stable, and the average total error is less than 10% when CR ranges from 0.3 m to 1.1 m, and the accuracies are the highest around 0.6 m. In view of the balance of efficiency and accuracy, the value of this parameter is usually set from 0.8 m to 1.1 m. In terms of the principle of optimal accuracy, the value of this parameter is usually set from 0.5 m to 0.7m.  After ground seed points are obtained, the remaining ground points are extracted using PTD. To analyze the sensitivity of the parameters (θ and d) on the improvement of filtering accuracy, S11 was selected as a representative in this section. We adjusted the two parameters as follows: θ ∈ {10°, 15°, 20°, 25°, 30°, 35°, 40°, 45°, 50°, 55°, 60°} , interval θ = 5° and d ∈ {0.5m, 1m, 1.5m} , interval d = 0.5m. For each parameter configuration, the total error was calculated and illustrated in Figure 23. It shows that the total error of the proposed algorithm is lower than PTDF under all parameter configurations, which indicates obtaining more ground seed points is an important way to improve the filtering accuracy of PTDF. This is in agreement with the previous studies that almost all improved PTDF focus on the densification of ground seed points. For PTDF and the proposed algorithm, the RMSE of the total error is 4.78% and 1.51% with the change of θ and d, respectively. It illustrates that the proposed algorithm is insensitive to the parameter selection even though the algorithm is parameter-dependent.

Four Filters with Higher Accuracy
Among the four filters with higher accuracy, the two filters developed by Mongus et al. [56] and Pingel et al. [24] used the optimal parameter values for each sample, but the two filters developed by Yang et al. [47] and Hu et al. [34] used the same parameter values or adaptive parameter values for each sample. Although the accuracy of the latter two filters is slightly lower than the accuracy of the first two filters, the latter two filters are more practical. Comprehensively considering practicality and accuracy, the performance of the proposed algorithm is only poorer than the two filters developed by Yang et al. [47] and Hu et al. [34]. Specifically, Yang et al. [47] fused on two filters based on different entities (i.e., points and segments) to extract ground points from LiDAR data. The algorithm is composed of two parts: (1) Smoothness-constrained segmentation is used to separate point clouds into point entities (e.g., trees) and segment entities (e.g., buildings) based on the area of each segment. (2) Filters based on point and segment entities are used in corresponding entities, respectively. Filter based on point entities calculates the geometric properties (i.e., height difference) of each point and its neighboring points to judge each point property (ground or non-ground). Filter based on point entities usually fails to detect ground points around break lines. Filter based on segment entities calculates the height difference of each segment and its neighboring segments to judge each segment property (ground or non-ground).

Four Filters with Higher Accuracy
Among the four filters with higher accuracy, the two filters developed by Mongus et al. [56] and Pingel et al. [24] used the optimal parameter values for each sample, but the two filters developed by Yang et al. [47] and Hu et al. [34] used the same parameter values or adaptive parameter values for each sample. Although the accuracy of the latter two filters is slightly lower than the accuracy of the first two filters, the latter two filters are more practical. Comprehensively considering practicality and accuracy, the performance of the proposed algorithm is only poorer than the two filters developed by Yang et al. [47] and Hu et al. [34]. Specifically, Yang et al. [47] fused on two filters based on different entities (i.e., points and segments) to extract ground points from LiDAR data. The algorithm is composed of two parts: (1) Smoothness-constrained segmentation is used to separate point clouds into point entities (e.g., trees) and segment entities (e.g., buildings) based on the area of each segment. on point entities calculates the geometric properties (i.e., height difference) of each point and its neighboring points to judge each point property (ground or non-ground). Filter based on point entities usually fails to detect ground points around break lines. Filter based on segment entities calculates the height difference of each segment and its neighboring segments to judge each segment property (ground or non-ground). Filter based on segment entities is capable of retaining the break lines, but it has a disadvantage in removing non-ground points from the segments mixed ground and non-ground points. The area of the mixed segment is usually small, and they can be filtered better by filter based on point entities. The fusion of filters based on different entities makes use of their advantages. However, the performance of the fused filter significantly depends on segmentation results, and parameters in the segmentation process are difficult to set. Compared to the filter proposed by Yang et al. [47], the possible reason that the proposed algorithm yields lower accuracy is that we only used the filter based on point entities. Hu et al. [34] used a step factor to control the pyramidal hierarchical structure to obtain more ground points in a single iteration. In addition, the terrain surface was progressively refined using thin plate spline (TPS) instead of TIN, and the bending energy function that explicitly depicts the surface smoothness was used to calculate the adaptive threshold automatically. In contrast to TIN, TPS is constructed to generate smooth and oscillation-free trend surface, and it can produce local maxima in the neighborhood of the control points [3,31]. Above features reduce the influence of terrain gradient on the accuracy of the filtering. However, at the top level of the pyramid, cell resolution is difficult to determine. A low resolution may degrade rugged terrain features, whereas a high resolution may fail to remove large objects. Compared to the filter proposed by Hu et al. [34], the possible reason that the proposed algorithm yields lower accuracy is that TPS is more suitable for terrain expression than TIN. To summarize, compared to the proposed algorithm, the aforementioned filters have slightly higher accuracy but lower practicality, due to some parameters that are difficult to set.

Conclusions
Automatic filtering is usually necessary prior to LiDAR application for both DTM generation and feature extraction. In this study, we present a high accuracy and easy-to-use filtering algorithm. CS is used to obtain high-quality initial provisional DTM. Based on the initial provisional DTM, statistical analysis is exploited to adaptively adjust the PTD parameter thresholds under different terrain scenes. Finally, the algorithm utilizes PTD to refine the initial provisional DTM. The proposed algorithm makes the best use of the advantages of CSF and PTDF, that is, (1) CS can provide high-quality initial provisional DTM for PTD with only an easy-to-set integer parameter RT, and the high-quality initial provisional DTM can provide a good foundation for parameter threshold estimation of PTD, (2) PTD can yield promising results during refinement of the initial terrain, especially for hilltops and steep slopes. Benchmark datasets with low and high point density are used to compare the performance of the proposed algorithm with those of other publicized filtering algorithms. Overall, the proposed algorithm produces a promising filtering result and is highly automatic. The fact may help the users to use LiDAR data in their own applications. In the future, we will aim to develop a more sophisticated filter by combining optical images and LiDAR data.