A Multi-Constraint Combined Method for Ground Surface Point Filtering from Mobile LiDAR Point Clouds

: Point cloud ﬁltering is an essential preprocessing step in 3D (three-dimensional) LiDAR (light detection and ranging) point cloud processing. The ﬁltering of mobile LiDAR scanning point clouds is much more challenging due to their non-uniform distribution, the large-scale of missing data areas and the existence of both large size objects and small land features. This paper proposes a new ﬁltering method that combines range constraint, slope constraint and angular position constraint to ﬁlter ground surface points from mobile LiDAR point clouds. Firstly, a cylindrical coordinate system (CCS) is established for each block of mobile LiDAR point clouds and three attributes of mobile LiDAR points, i.e., the angular position attribute (AA), longitudinal distance attribute (LA) and range attribute (RA), are computed. Then, the mobile LiDAR point clouds are structured into a grid according to the AA and LA. Finally, the point clouds are ﬁltered by a single cross-section ﬁlter (SCSF) using range constraint and slope constraint, followed by a multiple cross-section ﬁlter (MCSF) using range constraint and angular position constraint. Five datasets are used to validate the proposed method. The experimental results show that the proposed new ﬁltering method achieves an average type I error, type II error, and total error of 1.426%, 1.885%, and 1.622%, respectively. ﬁrst computes attributes, i.e., the angular position attribute, longitudinal distance attribute and range attribute, of MLS point clouds. A grid is then constructed using the angular position attribute and longitudinal distance attribute to establish the connectivity between point clouds. Single cross-section ﬁltering (SCSF) using the range constraint and slope constraint and multiple cross-section ﬁltering (MCSF) using angular position constraint and range constraint are sequentially applied to the constructed grid to get the ground surface points. Five datasets were used to validate the performance of the proposed method. The experimental results show that the average type I error, type II error, and total error of the proposed ﬁltering method in the ﬁve datasets were 1.426%, 1.885% and 1.622%, respectively. The average type I error, type II error and total error of the proposed ﬁltering method are relatively low compared with CSF and PMF.


Introduction
Mobile LiDAR scanning (MLS) is a state-of-the-art form of three-dimensional (3D) spatial data acquisition technology and has been widely used in road surface mapping [1], tunnel mapping [2], and road inventory [3], etc. The usage of MLS data can be mainly divided into two categories: topographical mapping of the road corridor and road object information extraction. In topographical mapping of the road corridor, ground surface points should be extracted from the raw MLS data. In road object information extraction, ground surface points are first obtained (for road marking extraction) or removed (for traffic sign and pole-like object extraction) to reduce the data volume and facilitate object extraction. In both categories of MLS data usage, ground surface point filtering plays an important role in data processing.
Many filtering algorithms have been proposed to cope with the ground surface point filtering problem of 3D LiDAR point clouds. The existing filtering methods can mainly be categorized into four groups: (1) slope-based filtering methods; (2) morphology-based filtering methods; (3) surface-based filtering methods; and (4) segmentation-based filtering methods. Vosselman [4] first proposed the slope-based filtering method, in which a point is classified as a ground point if all of the points within a distance threshold have a height difference smaller than the allowed maximum difference. Instead of using the same filter kernel for an entire dataset, Sithole [5] modified the slope-based filtering method in such a manner that the threshold varied with respect to the slope of the terrain and better performance was achieved in steep areas. Local slope in [4] was refined by the local linear regression criterion in [6]. The minimum height was first found in the neighborhood and points were assigned to the ground if their heights were compatible with the linear regression criterion. Then the approximate DEM (digital elevation model) was calculated according to the weighted mean of the ground points. Finally, the points with distances to the approximate DEM that met the threshold criterion were classified as final ground points.
The morphology-based filtering methods conduct morphological operations on the point clouds and use the height difference between points before and after morphological operations to judge whether the point is a ground point or not. The morphological operator "opening" was applied several times with different window sizes to filter ground points from airborne LiDAR scanning (ALS) data by Kilian et al. [7]. A window was moved over the dataset and weights depending on the window size were assigned to points higher than the deepest point in the applied window. The ground surface was interpolated under the consideration of the weights. Zhang et al. [8] developed a progressive morphological filter to remove non-ground points from ALS data. By gradually increasing the window size, non-ground points were removed progressively. Chen et al. [9] developed a similar progressive morphological filter but the slope was no longer assumed as constant. Pingel et al. [10] implemented a morphological filter with a linearly-increasing window and simple slope thresholding. Mongus et al. [11] filtered ground surface points from ALS data with features mapped from differential morphological profiles.
There are two types of surface-based filtering methods: surface interpolation methods and progressive TIN (triangulated irregular network) densification (PTD) methods. Kraus and Pfeifer [12] and Pfeifer et al. [13] used the full ALS data with equal weights to interpolate an initial surface, and the weight of each point was then updated based on the distance to the interpolated surface. The surface was iteratively interpolated until the weights of all points remained unchanged. Points within the buffer of the final interpolated surface were classified as ground points. An active shape model that acts like a rubber cloth with elasticity and rigidity was used to model the ground surface in [14,15]. The model was fitted by minimizing the energy function and any points within the buffer of the fitted model were accepted as ground surface points. Similar to the active shape model method, cloth simulation was used to interpolate the ground surface in [16]. Hu et al. [17] interpolated thin plate spline surfaces toward the ground in their adaptive surface filter (ASF). Bending energy was embedded in ASF to self-adapt the filter threshold automatically. Axelsson [18] proposed a PTD method to separate ground surface points from non-ground surface points. The seed ground points were used to create an initial TIN structure, and the points meeting distance and angle criteria were added to densify the TIN progressively. A two-step progressive TIN densification method was proposed by Sohn and Dowman [19]. In the downward densification step, four corner points were first chosen and triangulated, then the lowest points in each triangle were added to the triangulation. The downward densification process was repeated until no more points were found beneath each triangle. In the upward densification step, a buffer was defined for each triangle. The point within the buffer that forms the flattest tetrahedral was added to the triangulation. The final TIN represents the ground surface. Zhang and Lin [20] improved the traditional PTD filter by adding segmentation using smoothness constraint (SUSC) between TIN construction and densification to preserve ground points in areas with steep terrain. Zhao et al. [21] improved the traditional PTD filter by combining it with a morphological method to provide more ground seed points, adding simulated ground points to improve the quality of the initial TIN and performing downward densification before upward densification to deal with slope variations.
The segmentation-based filtering methods classify a segment rather than a single point as ground or non-ground based on local geometrical relations such as height, slope or curvature in the neighborhood. Brovelli et al. [22,23] proposed a method to filter ALS data according to the edges detected in the interpolated DSM (digital surface model). The points inside the convex hull determined by the edge cells were classified as objects in case their height was greater than or equal to the computed mean edge height. Tóvári and Pfeifer [24] extended the point-based filtering method [12] to a segmentation-based method to achieve better performance. The raw ALS data were first segmented into segments by region growing, and then a surface interpolation method similar to [12] was used to interpolate the ground surface.
The filtering algorithms mentioned above were mainly developed for ALS data. Compared with ALS data, MLS data have a much higher point density, which usually greatly decreases the efficiency of filtering algorithms. Besides, the point density of MLS data is not uniform in the data extent because of the different scanning distances, making it confusing to determine the window size and grid size in morphology-based and surface-based filtering algorithms, respectively. In some works [25,26], the authors attempted to calculate the MLS point density accurately to minimize the influence of the non-uniform distribution on MLS data processing. Last but not least, large-scale missing data areas exist inside the MLS data extent because of occlusion. The interpolation of large-scale missing data areas decreases the reliability of data, especially when ground points are mixed with points from vegetation and buildings [4]. Filtering algorithms for MLS data considering the data characteristics are needed.
However, fewer filtering methods are designed specifically for MLS data, taking into consideration the data characteristics. Wu et al. [27] proposed a method to extract ground surface points from MLS data by organizing data into vertical profiles. Candidate ground points were extracted from the vertical profiles by an alpha shape algorithm and ground points were refined by an elevation variance filter. Nurunnabi et al. [28] proposed a locally-weighted regression technique to fit the ground surface iteratively. Points within the specific height buffer of the final fitted ground surface were identified as ground surface points. Zhou et al. [29] proposed a scanline-based segmentation method to filter ground surface points from an ordered MLS point cloud.
In contrast to intentionally filtering ground surface points, in some works the ground points are extracted or removed first in order to identify road object information from MLS data. In these works, the points on the pavement are extracted or eliminated, leaving the ground points outside the pavement unhandled. For example, Guan et al. [30] proposed a method to extract road surface points using curbstone information on both sides of the road. Zheng et al. [31] removed the points on or around the ground by a piecewise elevation histogram segmentation method. These methods are mainly designed to extract or remove the points inside the road pavement, and they are not suitable for applications requiring full ground surface points within the data extent.
Besides filtering ground surface points from MLS data, in some works, the authors try to label the point clouds or classify the whole point clouds into different categories directly. Usually, ground surface points represent one of the categories. Yang and Dong [32] proposed a shape-based method to segment the MLS point clouds of urban scenes into different objects. Luo et al. [33] proposed a patch-based semantic labeling framework to label each point of colorized mobile LiDAR point clouds. Both semantic labeling and direct classification are much more challenging and time-consuming due to the scene complexity and the huge volume of MLS point clouds.
This paper proposes a new filtering method specifically designed to filter ground surface points from MLS data. In our method, a cylindrical coordinate system is established for each block of the mobile LiDAR point cloud and three attributes of mobile LiDAR point clouds, i.e., angular position attribute (AA), longitudinal distance attribute (LA) and range attribute (RA), are first computed. Then, the MLS point clouds are structured into a grid according to AA and LA. The MLS point clouds are quasi-uniformly distributed and large-scale missing data areas are avoided in such a grid. Unlike the existing filtering methods that use slope differences or elevation differences between ground surface points and non-ground surface points, the single cross-section filter (SCSF) using range constraint and slope constraint, and the multiple cross-sections filter (MCSF) using range constraint and angular position constraint are designed to filter ground surface points from MLS data.
The remainder of this paper is organized as follows. The datasets used to validate the proposed filtering method are first described in Section 2, followed by the detailed description of the proposed filtering method. Section 3 presents the experimental results, and the proposed filtering method is discussed in Section 4. Finally, Section 5 concludes this paper.

MLS Datasets
The mobile LiDAR point cloud data for the Jincheng highway in Beijing, China were collected by a mobile LiDAR system equipped with a RIEGL VUX-1 laser scanner. The speed of the mobile LiDAR system during data acquisition was about 40 km/h. The angular step width of laser scanning was 0.02 • and the distance between two scan lines was about 0.06 m. The total data coverage was about 57 km ( Figure 1). For convenience, five subsets with different ground surface types and different non-ground objects were selected as the test data from the whole dataset ( Figure 1).

MLS Datasets
The mobile LiDAR point cloud data for the Jincheng highway in Beijing, China were collected by a mobile LiDAR system equipped with a RIEGL VUX-1 laser scanner. The speed of the mobile LiDAR system during data acquisition was about 40 km/h. The angular step width of laser scanning was 0.02° and the distance between two scan lines was about 0.06 m. The total data coverage was about 57 km ( Figure 1). For convenience, five subsets with different ground surface types and different non-ground objects were selected as the test data from the whole dataset ( Figure 1). Different ground surface types and non-ground objects were contained in the selected test data. Figure 2 illustrates the cross-section views of the five datasets. The selected five subsets all have flat road pavements with moving cars on them (the pavements lie inside the red boxes in Figure 2). Besides a flat road pavement, Subset 1, Subset 3 and Subset 4 contain gently sloping ground surfaces, while Subset 2 and Subset 5 contain steep ground surfaces (① in Figure 2b and ② in Figure 2e, respectively). All five subsets contain trees and cars, especially Subset 1 and Subset 2. Lamps exist in Subset 1. Buildings exist in Subset 1 and Subset 4. Power lines exist in Subset 2, Subset 3 and Subset 5. The details of the five selected subsets of data are described in Table 1. Different ground surface types and non-ground objects were contained in the selected test data. Figure 2 illustrates the cross-section views of the five datasets. The selected five subsets all have flat road pavements with moving cars on them (the pavements lie inside the red boxes in Figure 2). Besides a flat road pavement, Subset 1, Subset 3 and Subset 4 contain gently sloping ground surfaces, while Subset 2 and Subset 5 contain steep ground surfaces ( 1 in Figure 2b and 2 in Figure 2e, respectively). All five subsets contain trees and cars, especially Subset 1 and Subset 2. Lamps exist in Subset 1. Buildings exist in Subset 1 and Subset 4. Power lines exist in Subset 2, Subset 3 and Subset 5. The details of the five selected subsets of data are described in Table 1.

MLS Datasets
The mobile LiDAR point cloud data for the Jincheng highway in Beijing, China were collected by a mobile LiDAR system equipped with a RIEGL VUX-1 laser scanner. The speed of the mobile LiDAR system during data acquisition was about 40 km/h. The angular step width of laser scanning was 0.02° and the distance between two scan lines was about 0.06 m. The total data coverage was about 57 km ( Figure 1). For convenience, five subsets with different ground surface types and different non-ground objects were selected as the test data from the whole dataset ( Figure 1). Different ground surface types and non-ground objects were contained in the selected test data. Figure 2 illustrates the cross-section views of the five datasets. The selected five subsets all have flat road pavements with moving cars on them (the pavements lie inside the red boxes in Figure 2). Besides a flat road pavement, Subset 1, Subset 3 and Subset 4 contain gently sloping ground surfaces, while Subset 2 and Subset 5 contain steep ground surfaces (① in Figure 2b and ② in Figure 2e, respectively). All five subsets contain trees and cars, especially Subset 1 and Subset 2. Lamps exist in Subset 1. Buildings exist in Subset 1 and Subset 4. Power lines exist in Subset 2, Subset 3 and Subset 5. The details of the five selected subsets of data are described in Table 1.  Besides x, y, and z coordinates, a precise GPS (global positioning system) timestamp is another attribute of each MLS point. The trajectory data, including a precise GPS timestamp and the position and attitude of the mobile LiDAR system platform at that time point, are delivered together with the MLS point clouds.  Besides x, y, and z coordinates, a precise GPS (global positioning system) timestamp is another attribute of each MLS point. The trajectory data, including a precise GPS timestamp and the position and attitude of the mobile LiDAR system platform at that time point, are delivered together with the MLS point clouds.

Computation of MLS Point Cloud Attributes
In order to compute the attributes of the MLS point clouds, i.e., AA, LA and RA, the cylindrical coordinate system was used. AA, LA and RA correspond to the angular position, longitudinal distance and radial distance in cylindrical coordinate system, respectively. The cylindrical coordinate system can be seen as an extension of a two-dimensional (2D) polar coordinate system. In a polar coordinate system, a point on a plane is determined by a distance from a reference point (origin) and an angle from a reference direction [34]. The distance from the reference point is called the radial distance, and the angle from the reference direction is called the angular position (angular coordinate, or azimuth). Figure 3a is an illustration of a polar coordinate system. In Figure 3a, A is the reference direction and O is the reference point (origin). The coordinates of point K1 in the polar coordinate system is determined by the radial distance from O to K1 (ρ1) and the angle from OA to OK1 (ϕ1). The coordinates of point K2 are determined by the distance from O to K2 (ρ2) and the angle from OA to OK2 (ϕ2).
A 3D cylindrical coordinate system combines the 2D polar coordinate system with an additional z-axis vertically or L-axis horizontally (this paper uses the L-axis horizontally) [35]. To construct a 3D cylindrical coordinate system, the reference plane R should first be obtained. On the reference plane R, the reference direction A and the reference point O are specified to construct a polar coordinate system. The longitudinal axis L starts from reference point O and is perpendicular to the reference plane R, pointing in the positive direction. Figure 3b is an illustration of a 3D cylindrical coordinate system. In Figure 3b, T is an assistant plane that contains point C and is parallel to the reference plane R. The intersection of T and longitudinal axis L ( O in Figure 3b) is the reference point of the polar coordinate system on plane T, and the reference direction is the same as that in plane R. The position of point C in a 3D cylindrical coordinate system can be denoted by the ordered triplet (ρ, ϕ, d), in which ρ and ϕ are the radial distance and angular position in the polar coordinate system, respectively, while d is the longitudinal position distance between the assistant plane T and reference plane R.

Computation of MLS Point Cloud Attributes
In order to compute the attributes of the MLS point clouds, i.e., AA, LA and RA, the cylindrical coordinate system was used. AA, LA and RA correspond to the angular position, longitudinal distance and radial distance in cylindrical coordinate system, respectively. The cylindrical coordinate system can be seen as an extension of a two-dimensional (2D) polar coordinate system. In a polar coordinate system, a point on a plane is determined by a distance from a reference point (origin) and an angle from a reference direction [34]. The distance from the reference point is called the radial distance, and the angle from the reference direction is called the angular position (angular coordinate, or azimuth). Figure 3a is an illustration of a polar coordinate system. In Figure 3a, A is the reference direction and O is the reference point (origin). The coordinates of point K1 in the polar coordinate system is determined by the radial distance from O to K1 (ρ1) and the angle from OA to OK1 (φ1). The coordinates of point K2 are determined by the distance from O to K2 (ρ2) and the angle from OA to OK2 (φ2).
A 3D cylindrical coordinate system combines the 2D polar coordinate system with an additional z-axis vertically or L-axis horizontally (this paper uses the L-axis horizontally) [35]. To construct a 3D cylindrical coordinate system, the reference plane R should first be obtained. On the reference plane R, the reference direction A and the reference point O are specified to construct a polar coordinate system. The longitudinal axis L starts from reference point O and is perpendicular to the reference plane R, pointing in the positive direction. Figure 3b is an illustration of a 3D cylindrical coordinate system. In Figure 3b, T is an assistant plane that contains point C and is parallel to the reference plane R. The intersection of T and longitudinal axis L ( in Figure 3b) is the reference point of the polar coordinate system on plane T, and the reference direction is the same as that in plane R. The position of point C in a 3D cylindrical coordinate system can be denoted by the ordered triplet ( , , ), in which and are the radial distance and angular position in the polar coordinate system, respectively, while is the longitudinal position distance between the assistant plane T and reference plane R. MLS point clouds are usually divided into blocks and the points are ordered by timestamp in each block. To establish a cylindrical coordinate system for each block, the plane that contains the first point P1 (earliest on timestamp) and is perpendicular to the trajectory data is selected as reference plane R. The intersection of reference plane R and the trajectory data is defined as the reference point (origin) O, and the z-direction is defined as the reference direction A. In this way, the polar coordinate system in the reference plane is established. The longitudinal axis L of the 3D cylindrical coordinate system starts at reference point O and is perpendicular to the reference plane R, with the direction along the trajectory direction as the positive direction.
Since the trajectory of MLS point cloud data can hardly be a straight line and the longitudinal axis L will not be parallel to the trajectory data, the conventional cylindrical coordinate system MLS point clouds are usually divided into blocks and the points are ordered by timestamp in each block. To establish a cylindrical coordinate system for each block, the plane that contains the first point P1 (earliest on timestamp) and is perpendicular to the trajectory data is selected as reference plane R. The intersection of reference plane R and the trajectory data is defined as the reference point (origin) O, and the z-direction is defined as the reference direction A. In this way, the polar coordinate system in the reference plane is established. The longitudinal axis L of the 3D cylindrical coordinate system starts at reference point O and is perpendicular to the reference plane R, with the direction along the trajectory direction as the positive direction.
Since the trajectory of MLS point cloud data can hardly be a straight line and the longitudinal axis L will not be parallel to the trajectory data, the conventional cylindrical coordinate system illustrated in Figure 3b should be modified to adapt it to the attributes of the computation of MLS point clouds. First, the assistant plane T used to figure out the cylindrical coordinates of point C will no longer be parallel to the reference plane R. Instead, T is redefined as the plane that contains point C and is perpendicular to the trajectory data. The intersection of plane T and trajectory data is the reference point (origin) of the polar coordinate system on plane T. The z-direction is still the reference direction of the polar coordinate system on plane T. In this way, the radial distance and the angular position of a point in an MLS point clouds can be computed. Second, the longitudinal position d is not the distance between plane T and plane R, but the trajectory length between plane R and plane T. Figure 4a is an illustration of the cylindrical coordinate system used to compute the attributes of MLS point clouds.
Remote Sens. 2017, 9, x FOR PEER REVIEW 7 of 19 illustrated in Figure 3b should be modified to adapt it to the attributes of the computation of MLS point clouds. First, the assistant plane T used to figure out the cylindrical coordinates of point C will no longer be parallel to the reference plane R. Instead, T is redefined as the plane that contains point C and is perpendicular to the trajectory data. The intersection of plane T and trajectory data is the reference point (origin) of the polar coordinate system on plane T. The z-direction is still the reference direction of the polar coordinate system on plane T. In this way, the radial distance and the angular position of a point in an MLS point clouds can be computed. Second, the longitudinal position is not the distance between plane T and plane R, but the trajectory length between plane R and plane T. Figure 4a is an illustration of the cylindrical coordinate system used to compute the attributes of MLS point clouds.  As mentioned above, AA, LA and RA correspond to the angular position, longitudinal distance and radial distance in a cylindrical coordinate system, respectively. To compute the attributes of a point P of MLS point clouds, denoted by ( , , , ), the Cartesian coordinates of the reference point in plane T, denoted as ( , , , , , ), are obtained by making perpendicular to the trajectory. The timestamp t of point P is used to find out the trajectory data between t − dt and t + dt so as to obtain the reference point in the polar coordinate system more efficiently. Then the RA of P (denoted as ) can be obtained by calculating the Euclidean distance between point P and reference point ,according to Equation (1). As illustrated in Figure 4b, the AA of P (denoted as ) can be estimated by the RA and the z difference between P and through the trigonometric function as described in Equation (2). The LA of P is calculated by accumulating the distance from to along the trajectory.

Filtering Method
The attributes of MLS point clouds, i.e., AA, LA and RA, are first computed according to Section 2.2. Then, a point in the MLS point clouds is denoted as ( , , , , , , ) where , , , , , , denotes timestamp, x coordinate, y coordinate, z coordinate, range attribute, angular position attribute, and longitudinal distance attribute, respectively. In order to improve the efficiency of the filtering algorithm, a grid is first constructed to establish the connectivity between the MLS points based on AA and LA. In the constructed grid, a column is a cross-section of the road. To extract the ground surface points, a single cross-section filter (SCSF) using range constraint and slope constraint is first applied to each As mentioned above, AA, LA and RA correspond to the angular position, longitudinal distance and radial distance in a cylindrical coordinate system, respectively. To compute the attributes of a point P of MLS point clouds, denoted by (t, x, y, z), the Cartesian coordinates of the reference point O in plane T, denoted as (x 0 , y 0 , z 0 ), are obtained by making PO perpendicular to the trajectory. The timestamp t of point P is used to find out the trajectory data between t − dt and t + dt so as to obtain the reference point in the polar coordinate system more efficiently. Then the RA of P (denoted as ρ) can be obtained by calculating the Euclidean distance between point P and reference point O , according to Equation (1). As illustrated in Figure 4b, the AA of P (denoted as ϕ) can be estimated by the RA and the z difference between P and O through the trigonometric function as described in Equation (2). The LA of P is calculated by accumulating the distance from O to O along the trajectory.
on the left side of trajectory . (2)

Filtering Method
The attributes of MLS point clouds, i.e., AA, LA and RA, are first computed according to Section 2.2. Then, a point in the MLS point clouds is denoted as (t, x, y, z, ρ, ϕ, d) where t, x, y, z, ρ, ϕ, d denotes timestamp, x coordinate, y coordinate, z coordinate, range attribute, angular position attribute, and longitudinal distance attribute, respectively. In order to improve the efficiency of the filtering algorithm, a grid is first constructed to establish the connectivity between the MLS points based on AA and LA. In the constructed grid, a column is a cross-section of the road. To extract the ground surface points, a single cross-section filter (SCSF) using range constraint and slope constraint is first applied to each cross-section. By using SCSF, most of the non-ground points will be filtered out. To refine the ground surface points, multiple cross-sections are used to fit the local surface. The points with an angular position smaller than the minimum angular position or larger than the maximum angular position of that cross-section and the points with a range far away from the fitted local surface are filtered out. The workflow of the filtering method proposed in this paper is illustrated in Figure 5. cross-section. By using SCSF, most of the non-ground points will be filtered out. To refine the ground surface points, multiple cross-sections are used to fit the local surface. The points with an angular position smaller than the minimum angular position or larger than the maximum angular position of that cross-section and the points with a range far away from the fitted local surface are filtered out. The workflow of the filtering method proposed in this paper is illustrated in Figure 5.

Grid Construction
A grid g is constructed to establish the connectivity between MLS points based on angular position attributes and longitudinal distance attributes. The grid g is constructed by the following three steps: • Finding the minimum angular position and the minimum longitudinal distance .

•
Determining the angular resolution and distance resolution .

•
Assigning the MLS point ( , , , , , , ) to the grid cell with row index and column index computed using Equations (3) and (4). If more than one point is assigned to the same grid index, the point with largest range attribute is preserved.

= Floor
, Since the angular step-width of a LiDAR scanner during a project will seldom be changed, the point clouds are quasi-uniformly distributed along the angular position direction even though the angular position is not strictly the same as the scanning angle in laser scanner data acquisition geometry. Besides, as MLS point clouds are usually divided into small blocks and the speed of the platform during a block can be regarded as steady, the MLS points are also quasi-uniformly distributed along the longitudinal distance direction. Compared with the grid created in the XY horizontal plane, the advantage of the grid created in this paper based on the angular position attribute and the longitudinal distance attribute is that the point clouds are quasi-uniformly distributed in the grid and large-scale missing data areas are avoided in such a grid. Figure 6 is an illustration of mobile LiDAR point clouds and the grid created based on the angular position attribute and the longitudinal distance attribute. Figure 6a shows MLS point clouds colored by elevation, while Figure 6b shows the established grid based on the angular position attribute and the longitudinal distance attribute colored by range attribute. In Figure 6b, white cells represent empty cells. It can be seen from Figure 6b that only a few empty cells exist in the grid, which indicates that the point clouds are quasi-uniformly distributed in such a grid and no large-scale empty cells exist in such a grid.

Grid Construction
A grid g is constructed to establish the connectivity between MLS points based on angular position attributes and longitudinal distance attributes. The grid g is constructed by the following three steps:

•
Finding the minimum angular position ϕ min and the minimum longitudinal distance d min .

•
Determining the angular resolution R a and distance resolution R d .

•
Assigning the MLS point P(t, x, y, z, ρ, ϕ, d) to the grid cell with row index cell row and column index cell col computed using Equations (3) and (4). If more than one point is assigned to the same grid index, the point with largest range attribute is preserved.
Since the angular step-width of a LiDAR scanner during a project will seldom be changed, the point clouds are quasi-uniformly distributed along the angular position direction even though the angular position is not strictly the same as the scanning angle in laser scanner data acquisition geometry. Besides, as MLS point clouds are usually divided into small blocks and the speed of the platform during a block can be regarded as steady, the MLS points are also quasi-uniformly distributed along the longitudinal distance direction. Compared with the grid created in the XY horizontal plane, the advantage of the grid created in this paper based on the angular position attribute and the longitudinal distance attribute is that the point clouds are quasi-uniformly distributed in the grid and large-scale missing data areas are avoided in such a grid. Figure 6 is an illustration of mobile LiDAR point clouds and the grid created based on the angular position attribute and the longitudinal distance attribute. Figure 6a shows MLS point clouds colored by elevation, while Figure 6b shows the established grid based on the angular position attribute and the longitudinal distance attribute colored by range attribute. In Figure 6b, white cells represent empty cells. It can be seen from Figure 6b that only a few empty cells exist in the grid, which indicates that the point clouds are quasi-uniformly distributed in such a grid and no large-scale empty cells exist in such a grid.

Single Cross-Section Filtering (SCSF)
As illustrated in Figure 7, there are two typical types of road cross-section: (1) road cross-sections with a roadside slope lower than the pavement; and (2) road cross-sections with a roadside slope higher than the pavement. O in Figure 7a,b represents the origin of the polar coordinate system used to compute the three attributes of MLS point clouds and the blue lines of dashes represent the range attributes of the MLS point clouds. In Figure 7a, the ground surface and objects between ② and ③ are occluded. The range attributes of points from ① to ② and from ③ to ④ will keep increasing obviously. In Figure 7b, the range attributes of points from ① to ② will keep increasing obviously. For points between ② and ③ in Figure 7b, if the next point is outside the circle which centers at O with a radius equal to the range attribute of the current point, the range attributes of the points will keep increasing. In the constructed grid g, a column is a cross-section of the road. For the points of a cross-section, the range attributes of the ground surface points on both sides of the minimum range point (the point with minimum range attribute in the cross-section) will keep increasing. Figure 8 is an illustration of a cross-section taken from real data. As illustrated in Figure 8b, the range attributes of the ground surface points on both sides of the minimum range point keep increasing. This characteristic can be used to filter ground surface points from MLS point clouds cross-section by cross-section.

Single Cross-Section Filtering (SCSF)
As illustrated in Figure 7, there are two typical types of road cross-section: (1) road cross-sections with a roadside slope lower than the pavement; and (2) road cross-sections with a roadside slope higher than the pavement. O in Figure 7a,b represents the origin of the polar coordinate system used to compute the three attributes of MLS point clouds and the blue lines of dashes represent the range attributes of the MLS point clouds. In Figure 7a, the ground surface and objects between 2 and 3 are occluded. The range attributes of points from 1 to 2 and from 3 to 4 will keep increasing obviously. In Figure 7b, the range attributes of points from 1 to 2 will keep increasing obviously. For points between 2 and 3 in Figure 7b, if the next point is outside the circle which centers at O with a radius equal to the range attribute of the current point, the range attributes of the points will keep increasing.

Single Cross-Section Filtering (SCSF)
As illustrated in Figure 7, there are two typical types of road cross-section: (1) road cross-sections with a roadside slope lower than the pavement; and (2) road cross-sections with a roadside slope higher than the pavement. O in Figure 7a,b represents the origin of the polar coordinate system used to compute the three attributes of MLS point clouds and the blue lines of dashes represent the range attributes of the MLS point clouds. In Figure 7a, the ground surface and objects between ② and ③ are occluded. The range attributes of points from ① to ② and from ③ to ④ will keep increasing obviously. In Figure 7b, the range attributes of points from ① to ② will keep increasing obviously. For points between ② and ③ in Figure 7b, if the next point is outside the circle which centers at O with a radius equal to the range attribute of the current point, the range attributes of the points will keep increasing. In the constructed grid g, a column is a cross-section of the road. For the points of a cross-section, the range attributes of the ground surface points on both sides of the minimum range point (the point with minimum range attribute in the cross-section) will keep increasing. Figure 8 is an illustration of a cross-section taken from real data. As illustrated in Figure 8b, the range attributes of the ground surface points on both sides of the minimum range point keep increasing. This characteristic can be used to filter ground surface points from MLS point clouds cross-section by cross-section. In the constructed grid g, a column is a cross-section of the road. For the points of a cross-section, the range attributes of the ground surface points on both sides of the minimum range point (the point with minimum range attribute in the cross-section) will keep increasing. Figure 8 is an illustration of a cross-section taken from real data. As illustrated in Figure 8b, the range attributes of the ground surface points on both sides of the minimum range point keep increasing. This characteristic can be used to filter ground surface points from MLS point clouds cross-section by cross-section. To filter out non-ground surface points in the nth cross-section, we first found the kth point ( , , , , , , ) with a minimum range attribute in the cross-section. The filtering is then accomplished by two steps: (1) Searching from to ( ) , the points with values that continue to increase from are classified as ground surface points; and (2) Searching from to ( ), the points with values that continue to increase from are classified as ground surface points.
( ) and ( ) are the maximum and minimum value of the crosssection, respectively. The filtering criteria can be summarized by Equation (5). In Equation (5), max{ , … , } and max{ , … , } are the maximum value of the ground surface points within ( , and , ), respectively.
After filtering using the range constraint, the slope constraint is further applied to each crosssection since with the slope constraint it is easy to filter out points on vertical objects and building roofs. The kth point with minimum is undoubtedly regarded as a ground surface point. Then the slope constraint is accomplished by two steps: (1) from to ( ), the slope value of each point ( , , , , , , ) is estimated by Equation (6)  To filter out non-ground surface points in the nth cross-section, we first found the kth point (t k n , ρ k n , ϕ k n , d k n , x k n , y k n , z k n ) with a minimum range attribute in the cross-section. The filtering is then accomplished by two steps: (1) Searching from ϕ k n to ϕ Max (n), the points with ρ values that continue to increase from ρ k n are classified as ground surface points; and (2) Searching from ϕ k n to ϕ Min (n), the points with ρ values that continue to increase from ρ k n are classified as ground surface points. ϕ Max (n) and ϕ Min (n) are the maximum and minimum ϕ value of the cross-section, respectively. The filtering criteria can be summarized by Equation (5). In Equation (5), max{ρ i+1 n , . . . , ρ k n } and max{ρ k n , . . . , ρ i−1 n } are the maximum ρ value of the ground surface points within (ϕ i n ,ϕ k n ] and [ϕ k n , ρ i n ), respectively.
After filtering using the range constraint, the slope constraint is further applied to each cross-section since with the slope constraint it is easy to filter out points on vertical objects and building roofs. The kth point with minimum ρ is undoubtedly regarded as a ground surface point. Then the slope constraint is accomplished by two steps: (1) from ϕ k n to ϕ Max (n), the slope value of each point (t  ) .

Multiple Cross-Section Filtering (MCSF)
In SCSF, not all points that satisfy Equation (5) are ground surface points. Figure 9b is an illustration of the result of SCSF. All points in Figure 9b satisfy Equation (5): the green points are ground surface points but the purple and black points are not ground surface points. There are two main types of non-ground points that exist in the result of SCSF: (1) the tree points or vegetation points in the margin of the cross-section (black points in Figure 9b); and (2) the tree points or vegetation points between ground surface points (purple points in Figure 9b).

Multiple Cross-Section Filtering (MCSF)
In SCSF, not all points that satisfy Equation (5) are ground surface points. Figure 9b is an illustration of the result of SCSF. All points in Figure 9b satisfy Equation (5): the green points are ground surface points but the purple and black points are not ground surface points. There are two main types of non-ground points that exist in the result of SCSF: (1) the tree points or vegetation points in the margin of the cross-section (black points in Figure 9b); and (2) the tree points or vegetation points between ground surface points (purple points in Figure 9b).  Figure 9c). The non-ground surface points between ground surface points could be eliminated by comparing the range attribute with the true surface range in the same grid cell. The range of non-ground surface points is smaller than the true surface range. The range of the true surface of a cross-section can be fitted by the nearby N1 cross-sections. For the ith grid cell in the nth cross-section, the fitted range is the maximum range of the ith grid cell of all the N1 nearby cross-sections. That means can be obtained by Equation (7). Then the fitted range is filtered by SCSF and the range of the filtered-out grid cell is interpolated by the nearby grid cell using linear interpolation. The minimum ground point angular position The marginal non-ground surface points in the nth cross-section could be eliminated by the angular position constraint. If the minimum angular position of ground point ϕ min (n) and the maximum angular position of ground point ϕ max (n) of the nth cross-section are known, the marginal points could be eliminated by eliminating the points with an angular position less than ϕ min (n) and greater than ϕ max (n) (the red dash line in Figure 9c). The non-ground surface points between ground surface points could be eliminated by comparing the range attribute with the true surface range in the same grid cell. The range of non-ground surface points is smaller than the true surface range.
The range of the true surface of a cross-section can be fitted by the nearby N1 cross-sections. For the ith grid cell in the nth cross-section, the fitted range is the maximum range of the ith grid cell of all the N1 nearby cross-sections. That means ρ i f n can be obtained by Equation (7). Then the fitted range is filtered by SCSF and the range of the filtered-out grid cell is interpolated by the nearby grid cell using linear interpolation. The minimum ground point angular position ϕ min (n) and the maximum ground point angular position ϕ max (n) of the nth cross-section is estimated by the nearby N2 fitted ground surface cross-sections according to Equations (8) and (9). The non-ground surface points between ground surface points could be eliminated in two steps: (1) for the ith point between ϕ k n and ϕ max (n), if ρ i n is smaller than the range of the (i − 1)th grid cell in the fitted ground surface, the point is classified as a non-ground point and eliminated; and (2) for the ith point between ϕ k n and ϕ min (n), if ρ i n is smaller than the range of the (i + 1)th grid cell in the fitted ground surface, the point is classified as a non-ground point and eliminated. The criteria for eliminating non-ground surface points between ground surface points is described in Equation (10).

Results
The five datasets described in Section 2.1 were filtered by the method proposed in this paper. The outlier points with high or low height values were removed manually before filtering. In grid construction, the angular resolution R a was set the same as the scanning angular step width (0.02 • ) and the distance resolution R d was set the same as the distance between two consecutive scan lines (0.06 m). The slope threshold T s was set to 45 • . In the MCSF, the number of nearby cross-sections used to interpolate the true ground surface points was set to 50 (about 3 m) and the number of nearby fitted cross-sections used to estimate minimum and maximum angular position of a cross-section was 250 (about 15 m). In order to compare the performance of the proposed filtering method with other filtering algorithms, the five datasets were also filtered by progressive morphological filtering (PMF) [8] and cloth simulation filtering (CSF) [16].
PMF uses a morphological filter to detect and remove non-ground points. The window size of the filter is gradually increased and the elevation difference threshold is determined based on the elevation variations of terrain, buildings and trees. The details of the PMF algorithm can be found in [8]. When filtering with PMF, the maximum window size, the slope value used to compute the height threshold, the initial distance, and the maximum distance were set to 20 m, 0.1, 0.05 m, and 0.5 m, respectively.
In CSF, the point clouds are first inverted and then a rigid cloth is used to cover the inverted surface. The locations of the cloth nodes are determined by analyzing the interactions between the cloth nodes and the corresponding LiDAR point. Ground surface points are obtained by comparing the original points with the generated cloth surface. The details of the CSF algorithm can be found in [16]. When filtering with CSF, the scene types of all the five datasets were set to a steep slope, and post-processing to handle steep slopes was used. The cloth grid resolution and the distance threshold were set to 1.3 and 0.02 m, respectively.
Reference ground surface points for the five datasets were created by the ground surface points filtering method in the commercial software TerraScan and through manual editing. First, the outlier-removed datasets were filtered by TerraScan. When filtering with TerraScan, the maximum building size, the terrain angle, the iteration angle and the iteration distance were set to 15 m, 45 • , 6 • , and 0.5 m, respectively. Then the filtering results were manually edited carefully to add ground surface points and remove non-ground surface points. Finally, the points with a distance smaller than 2 cm (a threshold by which almost all pavement points are classified as ground surface points) to the surface constructed by the edited ground surface points were all classified as ground surface points.
The top views and cross-section views of the filtering results created using the proposed method are shown in Figure 10. Since the MLS points in our filtering method are gridded and only one point is kept in a grid, many ground surface points are removed in the grid construction step. In order to evaluate the filtering results more reasonably, the points of filtered ground surface points were used to construct a TIN surface and the original points within 2 cm of the TIN surface were all reclassified as ground surface points. It can be seen from the top views in Figure 10 that the proposed method can correctly classify almost all pavement points as ground surface points and label objects such as cars above the pavement as non-ground points correctly. From the cross-section views in Figure 10 it can be seen that the proposed method can extract most of the ground surface points correctly even in challenging areas such as steep slope areas (Figure 10b2,d2,e2) and under dense vegetation (Figure 10a2,b2,c2). Non-ground points are seldom accepted as ground surface points in such natural landscapes.
To evaluate and compare the performance of the proposed method quantitatively, the type I error, type II error and total error [36] of the filtering results using CSF, PMF, and the proposed method were calculated based on the created reference data and are shown in Table 2. It can be seen from Table 2 that the method proposed in this paper achieves the best type I error in three of the five datasets (Subset 1, Subset 2 and Subset 5), the best type II error in three of the five datasets (Subset 1, Subset 3 and Subset 4) and the best total error in all five datasets. Compared with CSF, the proposed method performs better than CSF in type I error and total error in all five datasets. As for type II error, the proposed method performs better than CSF in gentle slope areas (Subset 1, Subset 3 and Subset 4), while performing slightly poorer in steep slope areas (Subsets 2 and 5). Compared with PMF, the proposed method achieves a comparable type I error in Subsets 1 and 3. The proposed method performs better than PMF for type I error in Subsets 2 and 5, while performing worse in Subset 4. Besides, the proposed method performs better than PMF in type II error, except in Subset 2. The results indicate that the proposed method can achieve good results for all the five datasets, and better results are achieved for ground surfaces with a gentle slope (Subset 1, Subset 3 and Subset 4) than steep ground surfaces (Subsets 2 and 5).
Remote Sens. 2017, 9, x FOR PEER REVIEW 13 of 19 The top views and cross-section views of the filtering results created using the proposed method are shown in Figure 10. Since the MLS points in our filtering method are gridded and only one point is kept in a grid, many ground surface points are removed in the grid construction step. In order to evaluate the filtering results more reasonably, the points of filtered ground surface points were used to construct a TIN surface and the original points within 2 cm of the TIN surface were all reclassified as ground surface points. It can be seen from the top views in Figure 10 that the proposed method can correctly classify almost all pavement points as ground surface points and label objects such as cars above the pavement as non-ground points correctly. From the cross-section views in Figure 10 it can be seen that the proposed method can extract most of the ground surface points correctly even in challenging areas such as steep slope areas (Figure 10b2,d2,e2) and under dense vegetation ( Figure  10a2,b2,c2). Non-ground points are seldom accepted as ground surface points in such natural landscapes.
To evaluate and compare the performance of the proposed method quantitatively, the type I error, type II error and total error [36] of the filtering results using CSF, PMF, and the proposed method were calculated based on the created reference data and are shown in Table 2. It can be seen from Table 2 that the method proposed in this paper achieves the best type I error in three of the five datasets (Subset 1, Subset 2 and Subset 5), the best type II error in three of the five datasets (Subset 1, Subset 3 and Subset 4) and the best total error in all five datasets. Compared with CSF, the proposed method performs better than CSF in type I error and total error in all five datasets. As for type II error, the proposed method performs better than CSF in gentle slope areas (Subset 1, Subset 3 and Subset 4), while performing slightly poorer in steep slope areas (Subsets 2 and 5). Compared with PMF, the proposed method achieves a comparable type I error in Subsets 1 and 3. The proposed method performs better than PMF for type I error in Subsets 2 and 5, while performing worse in Subset 4. Besides, the proposed method performs better than PMF in type II error, except in Subset 2. The results indicate that the proposed method can achieve good results for all the five datasets, and better results are achieved for ground surfaces with a gentle slope (Subset 1, Subset 3 and Subset 4) than steep ground surfaces (Subsets 2 and 5).  The computing time for attribute computation and filtering were recorded to analyze the time complexity of the proposed method. The proposed methods were implemented in C++ and run on an Intel(R) Core(TM) i7-3770 computer. The time costs are shown in Table 3. The results show that more than 80% of the time is spent on attribute computation. Finding the origin of the polar coordinates in attribute computation is time-consuming. However, the time cost is still tolerable, taking the huge amount of point numbers into consideration.   The computing time for attribute computation and filtering were recorded to analyze the time complexity of the proposed method. The proposed methods were implemented in C++ and run on an Intel(R) Core(TM) i7-3770 computer. The time costs are shown in Table 3. The results show that more than 80% of the time is spent on attribute computation. Finding the origin of the polar coordinates in attribute computation is time-consuming. However, the time cost is still tolerable, taking the huge amount of point numbers into consideration.

Accuracy
The proposed filtering method achieved the best total error in all the five datasets and achieves the best type I error and type II error in three of the five datasets. The proposed method achieved inspirational performance even in challenging areas such as steep slope areas and under dense vegetation. The average type I error (1.462%), average type II error (1.885%) and average total error (1.622%) of the proposed filtering method were relatively low compared with PMF and CSF. These results demonstrate that the proposed algorithm can be adapted to the filtering of MLS point clouds in different terrain types and different scene types.

Parameter Setting
The angular resolution R a and the distance resolution R d in grid construction can easily be set to the same as the scanning angular step width and distance between two consecutive scan lines. In this way, at least one point will be contained in a grid cell in its ideal condition and as few points as possible will be removed due to their falling in the same grid cell.
The number of nearby cross-sections N1 and N2 in MCSF are two parameters that affect the filtering results. N2 is determined by the maximum length of consecutive trees or buildings that prevent the LiDAR signal from hitting the ground surface. It is similar to the parameter maximum building size in the PTD algorithm. In order to filter out the trees and buildings on the margin of both sides of the road, we intended to set a large N2. It was much harder to choose an optimal N1, for the true ground surface could not be accurately interpolated with a small N1, while many ground points will be rejected in steep areas with a large N1 (as illustrated in Figure 11, more ground points are rejected in Figure 11b,c with the increasing N1). In 3D point cloud filtering, in order to generate the correct DEM with the filtered ground surface points, a small type II error is preferred, and having some ground surface points rejected is regarded as tolerable. In order to ensure that as few non-ground points as possible exist in the filtering result and the land details are kept as much as possible, we recommend N1 to be set to a value that covers around 3-5 m of the neighborhood.

Accuracy
The proposed filtering method achieved the best total error in all the five datasets and achieves the best type I error and type II error in three of the five datasets. The proposed method achieved inspirational performance even in challenging areas such as steep slope areas and under dense vegetation. The average type I error (1.462%), average type II error (1.885%) and average total error (1.622%) of the proposed filtering method were relatively low compared with PMF and CSF. These results demonstrate that the proposed algorithm can be adapted to the filtering of MLS point clouds in different terrain types and different scene types.

Parameter Setting
The angular resolution and the distance resolution in grid construction can easily be set to the same as the scanning angular step width and distance between two consecutive scan lines. In this way, at least one point will be contained in a grid cell in its ideal condition and as few points as possible will be removed due to their falling in the same grid cell.
The number of nearby cross-sections N1 and N2 in MCSF are two parameters that affect the filtering results. N2 is determined by the maximum length of consecutive trees or buildings that prevent the LiDAR signal from hitting the ground surface. It is similar to the parameter maximum building size in the PTD algorithm. In order to filter out the trees and buildings on the margin of both sides of the road, we intended to set a large N2. It was much harder to choose an optimal N1, for the true ground surface could not be accurately interpolated with a small N1, while many ground points will be rejected in steep areas with a large N1 (as illustrated in Figure 11, more ground points are rejected in Figure 11b,c with the increasing N1). In 3D point cloud filtering, in order to generate the correct DEM with the filtered ground surface points, a small type II error is preferred, and having some ground surface points rejected is regarded as tolerable. In order to ensure that as few nonground points as possible exist in the filtering result and the land details are kept as much as possible, we recommend N1 to be set to a value that covers around 3-5 m of the neighborhood.

THE Combination of Multiple Filtering Methods
Even though many algorithms have been proposed, ground surface point filtering from 3D LiDAR point clouds in large-scale complex scenes is still challenging. One of the trends of developing filtering algorithms is to combine multiple filtering methods to achieve better performance in challenging areas. Zhang and Lin [20] combined a segmentation-based method with a progressive TIN-densification method to achieve better performance. Yang et al. [37] combined a segmentationbased filter with a multi-scale morphological filter to remove large objects and preserve landscape details. This paper combines range constraint, slope constraint and angular position constraint to filter ground surface points from MLS point clouds and better performance was achieved compared with PMF and CSF.

Conclusions
This research proposed a new method to filter ground surface points from MLS point clouds. In contrast to other filtering methods that only use the xyz Cartesian coordinates, our method first computes attributes, i.e., the angular position attribute, longitudinal distance attribute and range attribute, of MLS point clouds. A grid is then constructed using the angular position attribute and longitudinal distance attribute to establish the connectivity between point clouds. Single cross-section filtering (SCSF) using the range constraint and slope constraint and multiple cross-section filtering (MCSF) using angular position constraint and range constraint are sequentially applied to the constructed grid to get the ground surface points. Five datasets were used to validate the performance of the proposed method. The experimental results show that the average type I error, type II error, and total error of the proposed filtering method in the five datasets were 1.426%, 1.885% and 1.622%, respectively. The average type I error, type II error and total error of the proposed filtering method are relatively low compared with CSF and PMF.
The novelty of the method proposed in this paper can be summarized in two thematic blocks. Firstly, the computed attributes of MLS point clouds. The angular position attribute, longitudinal distance attribute and range attribute provide an alternative for MLS data processing. Besides the SCSF and MCSF proposed in this research, some other algorithms may be more convenient and effective with the help of the angular position attribute, longitudinal distance attribute and range attribute. More algorithms such as pole-like object extraction and tree extraction may be proposed with the help of these attributes. Secondly, the SCSF using the range constraint and slope constraint and MCSF using the angular position constraint and range constraint. The experiments show that the combination of range constraint, slope constraint and angular position constraint can filter ground surface points from MLS point clouds effectively. This also provides an alternative to the combination of multiple filters to achieve better performance in challenging scenes.
Future research should find a more efficient algorithm to compute the attributes, since the computation is time-consuming at this time. Furthermore, the research should move forward to

THE Combination of Multiple Filtering Methods
Even though many algorithms have been proposed, ground surface point filtering from 3D LiDAR point clouds in large-scale complex scenes is still challenging. One of the trends of developing filtering algorithms is to combine multiple filtering methods to achieve better performance in challenging areas. Zhang and Lin [20] combined a segmentation-based method with a progressive TIN-densification method to achieve better performance. Yang et al. [37] combined a segmentation-based filter with a multi-scale morphological filter to remove large objects and preserve landscape details. This paper combines range constraint, slope constraint and angular position constraint to filter ground surface points from MLS point clouds and better performance was achieved compared with PMF and CSF.

Conclusions
This research proposed a new method to filter ground surface points from MLS point clouds. In contrast to other filtering methods that only use the xyz Cartesian coordinates, our method first computes attributes, i.e., the angular position attribute, longitudinal distance attribute and range attribute, of MLS point clouds. A grid is then constructed using the angular position attribute and longitudinal distance attribute to establish the connectivity between point clouds. Single cross-section filtering (SCSF) using the range constraint and slope constraint and multiple cross-section filtering (MCSF) using angular position constraint and range constraint are sequentially applied to the constructed grid to get the ground surface points. Five datasets were used to validate the performance of the proposed method. The experimental results show that the average type I error, type II error, and total error of the proposed filtering method in the five datasets were 1.426%, 1.885% and 1.622%, respectively. The average type I error, type II error and total error of the proposed filtering method are relatively low compared with CSF and PMF.
The novelty of the method proposed in this paper can be summarized in two thematic blocks. Firstly, the computed attributes of MLS point clouds. The angular position attribute, longitudinal distance attribute and range attribute provide an alternative for MLS data processing. Besides the SCSF and MCSF proposed in this research, some other algorithms may be more convenient and effective with the help of the angular position attribute, longitudinal distance attribute and range attribute. More algorithms such as pole-like object extraction and tree extraction may be proposed with the help of these attributes. Secondly, the SCSF using the range constraint and slope constraint and MCSF using the angular position constraint and range constraint. The experiments show that the combination of range constraint, slope constraint and angular position constraint can filter ground surface points from MLS point clouds effectively. This also provides an alternative to the combination of multiple filters to achieve better performance in challenging scenes.
Future research should find a more efficient algorithm to compute the attributes, since the computation is time-consuming at this time. Furthermore, the research should move forward to develop more algorithms such as pole-like object recognition with the help of the computed attributes proposed in this paper.