Airborne LiDAR data filtering based on geodesic transformations of mathematical morphology

: The capability of acquiring accurate and dense three-dimensional geospatial information that covers large survey areas rapidly enables airborne light detection and ranging (LiDAR) has become a powerful technology in numerous ﬁelds of geospatial applications and analysis. LiDAR data ﬁltering is the ﬁrst and essential step for digital elevation model generation, land cover classiﬁcation, and object reconstruction. The morphological ﬁltering approaches have the advantages of simple concepts and easy implementation, which are able to ﬁlter non-ground points effectively. However, the ﬁltering quality of morphological approaches is sensitive to the structuring elements that are the key factors for the ﬁltering success of mathematical operations. Aiming to deal with the dependence on the selection of structuring elements, this paper proposes a novel ﬁlter of LiDAR point clouds based on geodesic transformations of mathematical morphology. In comparison to traditional morphological transformations, the geodesic transformations only use the elementary structuring element and converge after a ﬁnite number of iterations. Therefore, this algorithm makes it unnecessary to select different window sizes or determine the maximum window size, which can enhance the robustness and automation for unknown environments. Experimental results indicate that the new ﬁltering method has promising and competitive performance for diverse landscapes, which can effectively preserve terrain details and ﬁlter non-ground points in various complicated environments.


Introduction
Light detection and ranging (LiDAR) is capable of rapidly collecting point clouds with accurate three-dimensional coordinates of the large survey area by integrating a high-frequency laser scanner, a global positioning system and an inertial measuring unit mounted on an aircraft.In contrast to remote sensing imagery that requires massive photogrammetric post-processing for calculating elevation, LiDAR point clouds explicitly describes the detailed 3D topographic profile of the Earth surface.In particular, LiDAR pulses can partly penetrate tree canopy and are not affected by shadows and surrounding light conditions.LiDAR performs well in some special areas such as forests, marshes and deserts, where other survey technologies experience difficulty [1][2][3].Therefore, LiDAR has been widely used in numerous fields of geospatial application such as topographic mapping, land management, urban planning, natural hazard assessment, forestry, geology, civil engineering, hydrology, communications and archaeology [4][5][6][7].
Over the past decades, LiDAR data filtering has attracted a large number of researches around the world.Various significant filtering approaches have been investigated and proposed, which can be categorized into different groups, such as slope-based approaches, triangulated irregular network (TIN) densification approaches, interpolation-based approaches and morphological approaches [9,[16][17][18]22,23].Slope-based approaches compare the heights of neighboring points, and separate ground from non-ground points by assuming that objects have larger slopes or elevation differences than terrains with smooth surfaces [2,14,[24][25][26][27].Interpolation-based approaches and TIN densification approaches create an initial ground surface using seed points, and then progressively interpolate or densify to approximate the ground surface under strict distance or angle constraints [10,[28][29][30][31][32].Morphological approaches usually involve a series of morphological erosion, dilation, opening or closing operations to separate ground from non-ground points [8,18,[33][34][35].The operations can be combined to obtain different spatial features.The structuring elements adopted in those morphological operations have great influence on filtering results.In addition, some researchers combine different approaches for filtering point clouds.For instance, Hui et al. [16] combine the progressive morphological approach and multi-level interpolation approach.Yang et al. [36] combine the segment-based filtering and multi-scale morphological filtering.Chen et al. [37] embed a region growing-based segmentation process in multi-resolution hierarchical classification to optimize the search for initial ground seeds.Zhao et al. [12] combine the progressive TIN densification filtering algorithm with a morphological method for complex forest sites.
Morphological approaches have the advantages of simple concepts and easy implementation, which is able to filter non-ground points effectively [8,33,36,38,39].The definition of the filtering window, also called structuring element in mathematical morphology, is important for successful filtering.On one hand, objects can be filtered out only if the adopted structuring element is sufficiently large.On the other hand, large structuring elements are likely to cut off abrupt changes in terrain heights such as cliffs, mountain ridges and peaks.Filtering performance is sensitive to the selection of window size [12,36,37].Hence, to deal with objects with different sizes, most morphological approaches need to adopt a series of windows with different sizes for filtering progressively, and the maximum window size depends on prior knowledge of the largest object in the survey area.
Previous studies have made significant contributions, and some existing approaches have been successfully applied.However, LiDAR data filtering is still challenging because previous filtering algorithms often have some limitations, especially when facing varied landscapes with complex objects and many abrupt changes in terrain heights [1,4,10,11,15,16,21,[40][41][42].Moreover, many algorithms need to determine several parameters such as filtering windows that have significant effects on filtering success.Normally, window sizes have to be tested empirically, which is difficult for varied environments [4,11,14].A novel filter of LiDAR point clouds based on geodesic transformations of mathematical morphology is proposed in this paper to address the problem discussed above.Geodesic transformations apply an iterative process until stability has been reached, eliminating the need of defining window sizes.Thus, this algorithm adopts the elementary filtering window without considering how to select different window sizes or determine the maximum window size, which enhance the robustness and automation for unknown environments.

Methodology
The proposed morphological filtering method based on geodesic transformations is explained by the following subsections.Section 2.1 introduces the reorganizing data and removing low outliers.Morphological gradients are employed as the indicator for different changes of terrain and objects, which are presented in Section 2.2.The morphological primitive created by combining geodesic transformations and the detailed implementation procedure are explained in Section 2.3.The experimental approach for assessing the filtering algorithm is introduced in Section 2.4.

Reorganizing Data and Removing Low Outliers
Raw LiDAR point clouds are usually recorded in the time order of acquiring laser pulse echoes, and a survey strip generally contains millions of LiDAR points.Irregularly distributed point clouds need to establish efficient data organization, which is the basis of subsequent spatial analyses and calculations.This study employs a simple grid index to restructure point cloud to utilize the simplicity of regular grid and avoid the accuracy loss of interpolation.A regular index grid can be constructed based on the average point spacing and the x, y bounding box of point clouds as shown in Figure 1.Each LiDAR point is stored in the corresponding grid cell according to its coordinates.During subsequent processes such as neighboring analysis, spatial search is performed around the index grid.Then, the corresponding grid cells are determined to select the contained LiDAR points for relevant calculations.
Remote Sens. 2017, 9, 1104 3 of 20 The proposed morphological filtering method based on geodesic transformations is explained by the following subsections.Section 2.1 introduces the reorganizing data and removing low outliers.Morphological gradients are employed as the indicator for different changes of terrain and objects, which are presented in Section 2.2.The morphological primitive created by combining geodesic transformations and the detailed implementation procedure are explained in Section 2.3.The experimental approach for assessing the filtering algorithm is introduced in Section 2.4.

Reorganizing Data and Removing Low Outliers
Raw LiDAR point clouds are usually recorded in the time order of acquiring laser pulse echoes, and a survey strip generally contains millions of LiDAR points.Irregularly distributed point clouds need to establish efficient data organization, which is the basis of subsequent spatial analyses and calculations.This study employs a simple grid index to restructure point cloud to utilize the simplicity of regular grid and avoid the accuracy loss of interpolation.A regular index grid can be constructed based on the average point spacing and the x, y bounding box of point clouds as shown in Figure 1.Each LiDAR point is stored in the corresponding grid cell according to its coordinates.During subsequent processes such as neighboring analysis, spatial search is performed around the index grid.Then, the corresponding grid cells are determined to select the contained LiDAR points for relevant calculations.Low outliers often appear in point clouds, which are generally caused by the laser pulses reflected by some objects several times or the laser rangefinder malfunction [9].Low outliers influence the filtering results and the quality of DEM generation, and are removed based on the fact of that these outliers are commonly far below the surrounding points and have few neighbors with similar heights.
To remove low outliers, we utilize the morphological black top-hat (BTH) transformation, which combines two basic morphological operations: dilation and erosion.The dilation operation δ for one considered point (x, y) in the set f of point clouds is defined as [43]: where B represents the window, also referred to as structuring element (SE) in mathematical morphology, and the domains of f and B are Df and DB, respectively.On the contrary, the erosion ε is to obtain the minimum from the heights of neighboring points, which is defined as: Combining dilation and erosion, the morphological BTH transformation using SE B on the set f is defined as: Low outliers often appear in point clouds, which are generally caused by the laser pulses reflected by some objects several times or the laser rangefinder malfunction [9].Low outliers influence the filtering results and the quality of DEM generation, and are removed based on the fact of that these outliers are commonly far below the surrounding points and have few neighbors with similar heights.
To remove low outliers, we utilize the morphological black top-hat (BTH) transformation, which combines two basic morphological operations: dilation and erosion.The dilation operation δ for one considered point (x, y) in the set f of point clouds is defined as [43]: where B represents the window, also referred to as structuring element (SE) in mathematical morphology, and the domains of f and B are D f and D B , respectively.On the contrary, the erosion ε is to obtain the minimum from the heights of neighboring points, which is defined as: Combining dilation and erosion, the morphological BTH transformation using SE B on the set f is defined as: Because of the low outlier characteristics of extremely low elevation and scattering distribution, the set L of low outliers can be detected by the criteria shown in Formula (4).If the BTH transformation result of point p i is larger than a threshold T1, the point p i has local extremely low elevation in comparison to the surrounding points.Because the low outliers have few close neighboring points, the sum of neighboring points within a window is used to judge whether p i is scattered.The window is predefined by the distribution of low outliers in practical data.Since the low outliers are far lower than surrounding point clouds, the thresholds in Formula (4) can be set by trial and error.
where [BTH B (f )](p i ) is the BTH transformation result of a considered point p i , N i is the set of neighboring points within a 7 × 7 window, Z i and Z j are the elevations of points p i and p j , and thresholds T1, T2 and T3 are assigned to be 6 m, 5 m and 6 points, respectively.Empty grid cells are expected to appear due to the irregular distribution of point clouds.Thus, a filling process for the index grid is implemented, which can suppress the influence of empty grid cells.After removing low outliers, each empty grid cell is assigned with a virtual point whose elevation is equal to the minimum elevation within a 3 × 3 neighboring window because the lowest point has most possibility of representing ground.Since the cell size of the index grid is set by the average point spacing of point clouds, most empty grid cells can be filled with a 3 × 3 window.The grid index is only used for neighbor search, and all relevant calculations are performed on the points contained in the corresponding grid cells.

Calculating Morphological Gradients
Changes in elevation are the main cues in search of objects because ground surfaces are commonly continuous and smooth in most areas.Morphological gradients can be easily calculated without considering the direction of elevation changes, and thus suitable for irregularly distributed LiDAR point clouds.Morphological gradients include several numerical expressions such as the internal and external gradients [43].After the erosion operation using the elementary 3 × 3 SE B, the difference from the original height is the internal gradient ρ − , which is defined as follows.
By contrast, external gradient ρ + is the height difference from the dilation operation, and is defined as: The two morphological gradients above indicate different kinds of transition points of changing heights.Lower jump points, such as points a and c in Figure 2, commonly have larger external gradients than their neighbors, whereas higher jump points, such as points b and d in Figure 2, have larger internal gradients than their neighbors.
The two morphological gradients above indicate different kinds of transition points of changing heights.Lower jump points, such as points a and c in Figure 2, commonly have larger external gradients than their neighbors, whereas higher jump points, such as points b and d in Figure 2, have larger internal gradients than their neighbors.However, drastic terrain changes, such as point e, exhibit large gradients.To differentiate low transition points such as points a and c from steep slope points such as point e, the gradient of external gradients, which is denoted by ∆(ρ + B ), can be defined as follows: Low transition points such as points a and c have higher ∆(ρ + B ) than the slope points such as point e, due to the gradually changing of terrain elevation.

Filtering Point Clouds by Geodesic Transformations
Most of the existing morphological filtering algorithms are based on the erosion and dilation operations with specific structuring elements that move around the input data, and the sizes of structuring elements, which are the key to successful filtering, need to adjust according to the distributed objects in the survey area.In this study, we use geodesic transformations with no necessity of selecting specific structuring elements, which are restricted to the elementary erosion and dilation operations with an elementary 3 × 3 window.
A geodesic transformation involves two datasets, i.e., a marker set f and a geodesic mask g.The geodesic mask is employed as the limit to the propagation of erosion or dilation of the marker dataset.Elementary geodesic erosion, which is denoted by ε (1) g ( f ), is that the marker dataset f is eroded by the elementary window, and the eroded results are forced to remain in the scope indicated by mask g [43].
Similarly, the elementary geodesic dilation of an input set f within a geodesic mask g is defined as: After the finite iteration of geodesic transformations, the input dataset converges to a stable status due to the mask limit.The geodesic erosion of f with a mask g that iterates i times until reaching stability, which is called reconstruction by erosion, is denoted by R ε g ( f ): Similarly, the geodesic dilation of f with a mask g that iterates i times until reaching stability, which is called reconstruction by dilation, is denoted by R δ g ( f ): Morphological reconstruction operations, i.e., reconstruction by erosion and reconstruction by dilation, are combined to create a new morphological primitive that is the kernel of the proposed filtering algorithm.The cross-section profile of the utilized morphological primitive is shown in Figure 3. Except the height h and slope θ need to be predefined, the width, i.e., the size of primitives is unnecessarily considered, which can self-adaptively fit various objects (e.g., objects a and c in Figure 3) closely and naturally by combining morphological reconstruction operations.The boundaries around objects commonly exhibit abrupt changes in gradients.Thus, we determine portions a and c in Figure 3 as non-ground objects because they fit the two sides of the primitive with changing value h and are higher than the entire primitive.The gradual change of terrain, such as portion b, is retained because no drastic gradient occurs.If some portions only have one side of large change (e.g., cliffs) or some parts beneath the primitive, the portions, such as d in Figure 3, are classified as ground.The primitive proposed above is realized by a series of geodesic transformations with the following steps.
(1) The low transition point set D with abrupt gradient changes that are larger than threshold h is detected as shown by Formula ( 12) and the red points in Figure 4a.
(2) The surrounding points around D are eroded by geodesic erosion if the points have a height difference h higher than the points of D, as shown by Formula (13).Elevation after erosion is assigned to the elevation of the points of D that is added to h, as shown by the blue points in Figure 4b.
where Ni is the set of neighboring points of pi within the 3 × 3 window, and Zi and Zj are the elevations of points pi and pj, respectively.
(3) A morphological reconstruction operator by erosion proceeds using the above point set f as the marker set.During each time of geodesic erosion, for surrounding point pi around pj that is eroded on the last time of geodesic erosion, the eroded elevation of pi is assigned to the elevation of pj that is added to ∆h if pi has a height difference ∆h higher than pj, as described by the following formula.
where ∆h = pointSpacing × tanθ, pointSpacing is the average point spacing of point clouds and θ is the slope of the morphological primitive.
After the morphological reconstruction operator by erosion is accomplished, the elevations of the eroded points are assigned to the elevations of low transition points, as shown in Figure 4c.The primitive proposed above is realized by a series of geodesic transformations with the following steps.
(1) The low transition point set D with abrupt gradient changes that are larger than threshold h is detected as shown by Formula ( 12) and the red points in Figure 4a.
(2) The surrounding points around D are eroded by geodesic erosion if the points have a height difference h higher than the points of D, as shown by Formula (13).Elevation after erosion is assigned to the elevation of the points of D that is added to h, as shown by the blue points in Figure 4b.
where N i is the set of neighboring points of p i within the 3 × 3 window, and Z i and Z j are the elevations of points p i and p j , respectively.
(3) A morphological reconstruction operator by erosion proceeds using the above point set f as the marker set.During each time of geodesic erosion, for surrounding point p i around p j that is eroded on the last time of geodesic erosion, the eroded elevation of p i is assigned to the elevation of p j that is added to ∆h if p i has a height difference ∆h higher than p j , as described by the following formula.
where ∆h = pointSpacing × tanθ, pointSpacing is the average point spacing of point clouds and θ is the slope of the morphological primitive.
After the morphological reconstruction operator by erosion is accomplished, the elevations of the eroded points are assigned to the elevations of low transition points, as shown in Figure 4c.
[  () where ∆h = pointSpacing × tanθ, pointSpacing is the average point spacing of point clouds and θ is the slope of the morphological primitive.
After the morphological reconstruction operator by erosion is accomplished, the elevations of the eroded points are assigned to the elevations of low transition points, as shown in Figure 4c.(4) Terrain discontinuity features, such as portion d in Figure 4, are restored by a morphological reconstruction by dilation.The marker set is shown as the blue point in Figure 4c and determined as follows.
where Z ε i is the elevation after erosion, Z i is the original elevation, and N(p i ) is the set of neighboring points within the 3 × 3 window.
Subsequently, geodesic dilation is performed.If the terrain exhibits significant changes as shown in Figure 5, the terrain points, such as the blue points, are not retrieved to the original surface.Thus, if the eroded surface has gradual changes regarding the surrounding intact points, the points on the gradual surface are assigned to the original surface after the geodesic dilation.For a point p i ∈ f, its elevation Z δ i after dilation is defined as follows.
where Z ε i is the elevation after erosion, Z i is the original elevation, and Z j is the maximum elevation of neighboring points within the 3 × 3 window.Thresholds T 1 and T 2 are set to 0.2 × pointSpacing and 10 m, respectively.Geodesic dilation iterates until reaching stability, which is restricted to the eroded points above.(4) Terrain discontinuity features, such as portion d in Figure 4, are restored by a morphological reconstruction by dilation.The marker set is shown as the blue point in Figure 4c and determined as follows.
where Z ε i is the elevation after erosion, Zi is the original elevation, and N(pi) is the set of neighboring points within the 3 × 3 window.Subsequently, geodesic dilation is performed.If the terrain exhibits significant changes as shown in Figure 5, the terrain points, such as the blue points, are not retrieved to the original surface.Thus, if the eroded surface has gradual changes regarding the surrounding intact points, the points on the gradual surface are assigned to the original surface after the geodesic dilation.For a point pi ∈ f, its elevation Z δ i after dilation is defined as follows.
where Z ε i is the elevation after erosion, Zi is the original elevation, and Zj is the maximum elevation of neighboring points within the 3 × 3 window.Thresholds T1 and T2 are set to 0.2 × pointSpacing and 10 m, respectively.Geodesic dilation iterates until reaching stability, which is restricted to the eroded points above.The dilated points are shown in Figure 4d.Non-ground points are determined based on the difference of elevations before and after the geodesic transformations.The points are classified as objects if the elevation difference is larger than the threshold that is set to the height h of the proposed morphological primitives.
(5) The creation of the filtering primitive is based on the hypothesis that the boundary around The dilated points are shown in Figure 4d.Non-ground points are determined based on the difference of elevations before and after the geodesic transformations.The points are classified as objects if the elevation difference is larger than the threshold that is set to the height h of the proposed morphological primitives.
(5) The creation of the filtering primitive is based on the hypothesis that the boundary around objects has drastic gradient changes.In practice, several objects are partially attached to the terrain, for example, the buildings on steep slopes as shown in Figure 6a, and the bridges connected to ground surface as shown in Figure 6b.The primitive with an isotropic structuring element cannot work for these kinds of objects.To deal with this problem, the primitive with elementary directional structuring elements, that is, the 3 × 1 and 1 × 3 windows, is employed for filtering the attached objects.Directional filtering has the same geodesic erosion process as isotropic filtering, including the above steps (1)-( 3).For the following geodesic dilations, directional filtering has looser constraints than isotropic filtering to avoid eliminating ground points because the terrain continuity decreases along the 1D directions.Thus, the geodesic dilation in directional filtering replaces the Formulas ( 15) and ( 16) with Formulas ( 17) and (18).
where Z ε i is the elevation after erosion, Z i is the original elevation, Z j is the maximum elevation of neighboring points within the elementary directional window, and threshold T is set to pointSpacing.where Z ε i is the elevation after erosion, Zi is the original elevation, Zj is the maximum elevation of neighboring points within the elementary directional window, and threshold T is set to pointSpacing.

Test Data
The test data provided by the International Society for Photogrammetry and Remote Sensing (ISPRS) Working Group III/3 are employed to evaluate the performance of our algorithm.The test regions are situated in the Vaihingen/Enz test field and Stuttgart city, which consist of diverse urban and rural landscapes.The LiDAR data were collected by an Optech ALTM scanner with a point spacing of 1-1.5 m in urban areas and 2-3.5 m in rural areas.Fifteen sample data have been selected from the survey data for the contained challenging features described in Table 1, and the reference data are manually produced by ISPRS for quantitative evaluation of filtering.

Test Data
The test data provided by the International Society for Photogrammetry and Remote Sensing (ISPRS) Working Group III/3 are employed to evaluate the performance of our algorithm.The test regions are situated in the Vaihingen/Enz test field and Stuttgart city, which consist of diverse urban and rural landscapes.The LiDAR data were collected by an Optech ALTM scanner with a point spacing of 1-1.5 m in urban areas and 2-3.5 m in rural areas.Fifteen sample data have been selected from the survey data for the contained challenging features described in Table 1, and the reference data are manually produced by ISPRS for quantitative evaluation of filtering.To assess the performance of this proposed algorithm on data with different characteristics, the test data, which are derived from the Actueel Hoogtebestand Nederland (AHN) (http://www.ahn.nl), are also used for filtering evaluation.The AHN LiDAR data are obtained by the collaboration of the provinces, central government, and water authorities of Netherlands, and have been classified through an industrial procedure.The two adopted test datasets are situated at the outskirt and downtown of Meerssen City.Dataset 1 is located in the outskirt, which contains dense vegetation on slopes and mixed buildings and vegetation on abruptly changing ground surface.Dataset 2 is located in the downtown, which includes a complex distribution of diverse buildings and vegetation on uneven terrains.The density of point clouds is about 1.4 points/m 2 .Datasets 1 and 2 contain 333,729 and 357,363 points, respectively.

Parameter Setting
For filtering the data with different density, resolution or accuracy, this proposed algorithm has three parameters to set by the data characteristics, which are the cell size of the index grid, height h and slope θ of the morphological primitive.
To represent the neighboring relation between point clouds, the cell size of the index grid is set to be slightly smaller than the average point spacing of point clouds, which can avoid that too many points are contained in each grid cell.Meanwhile, the filling process of the index grid as mentioned in Section 2.1 can suppress the influence of empty grid cells.Therefore, the cell size of the index grid is set to 1 m × 1 m for urban data, and 2 m × 2 m for rural data provided by ISPRS.The size of grid cell is set to be 0.5 m × 0.5 m for the grid index of AHN data.
The height h of the morphological primitive is set to be the smallest height difference between ground and objects which can be distinguished in the point clouds.Since the AHN data have higher resolution and accuracy than the ISPRS data, very low objects such as shrub need to be distinguished from AHN data.Therefore, the height h is set to be 0.3 × pointSpacing for the ISPRS data and 0.1 × pointSpacing for the AHN data.
The slope θ of the morphological primitive is defined by the relation of the size and height of objects.There is always a trade-off between what is desirable in a classification scheme as discussed in the reference [9].If the low objects are desired to remove to great extent, the slope θ is set to be relatively small.On the contrary, the slope θ is set to be relatively large if the filtering performance biased towards preserving terrain changes.Therefore, the slope θ is set to be arctan(0.5)for the ISPRS data and arctan(0.3)for the AHN data, which have many low objects with relatively large sizes to remove.

Evaluation Method
The quantitative evaluation of filtering algorithms is commonly based on three types of errors: omission, commission and total errors [9].Omission error, which is also called type I error, is the percentage of ground portions removed as objects.Commission error, which is also called type II error, is the percentage of object portions misclassified as ground.Total error is the overall percentage of all misclassified points in point clouds.A cross-matrix of filtering result can be obtained (Table 2), and then the three types of errors are calculated as Formulas ( 19)- (21).
Type II error = c/(c + d) Total error = (b + c)/(a Table 2. Cross-matrix of filtering result.

Ground Points Object Points
Reference Ground points a b Object points c d

Results and Discussion
The main contribution of this research is the avoidance of selection of specific structuring elements that is necessary and crucial in most filtering methods such as classic morphology-based algorithms.To verify the feasibility and effectiveness of this algorithm, the results are compared with the classic morphology-based algorithms published in [44].Li et al. [44] employs the morphological top-hat transformation with progressively increased windows and the constraints along directions to reduce the removal of protruding terrain features when using large windows.Comparing with other popular approaches, Li et al. [44] obtains the reliable and competitive performance for varied survey areas.Additionally, the comparison is also made with the filtering performance of the widely used commercial software Terrasolid TerraScan, which requires manual tuning of relevant parameters for different environments [40,45].The comparison of filtering results is shown as Figures 7-9.
Remote Sens. 2017, 9, 1104 10 of 20 reduce the removal of protruding terrain features when using large windows.Comparing with other popular approaches, Li et al. [44] obtains the reliable and competitive performance for varied survey areas.Additionally, the comparison is also made with the filtering performance of the widely used commercial software Terrasolid TerraScan, which requires manual tuning of relevant parameters for different environments [40,45].The comparison of filtering results is shown as Figures 7-9.The comparison shows that the proposed method tends to suppress omission error, and obtain a relatively lower average total error.However, the commission error is high compared to the other methods.Thus, our method classifies more non-ground as ground points than the other two methods, while fewer ground points are removed from the dataset as they are classified as non-ground points.In the present case, our filtering method is biased in favor of minimizing Type I errors as shown in Figure 10.The comparison shows that the proposed method tends to suppress omission error, and obtain a relatively lower average total error.However, the commission error is high compared to the other methods.Thus, our method classifies more non-ground as ground points than the other two methods, while fewer ground points are removed from the dataset as they are classified as non-ground points.In the present case, our filtering method is biased in favor of minimizing Type I errors as shown in Figure 10.The comparison shows that the proposed method tends to suppress omission error, and obtain a relatively lower average total error.However, the commission error is high compared to the other methods.Thus, our method classifies more non-ground as ground points than the other two methods, while fewer ground points are removed from the dataset as they are classified as non-ground points.In the present case, our filtering method is biased in favor of minimizing Type I errors as shown in Figure 10.Filtering sample 41 with a cluster of low outliers caused by multi-path error has the largest omission error.A hypothesis of this filter is that low outliers scatter, and the amount of close low outliers is supposed to be smaller than six points, which is true for most data.However, the cluster of low outliers in sample 41 contains more than six points, which result in the remaining low outliers.During filtering with directional structuring element, the ground points between the middle low outliers and the data edge are incorrectly eliminated.Figure 13 shows the two results with and without directional structuring elements.The omission, commission and total errors when filtering without the directional structuring element are 15.42%, 3.00% and 9.20%, respectively.The omission error in the middle is eliminated as shown in Figure 13b.Filtering sample 41 with a cluster of low outliers caused by multi-path error has the largest omission error.A hypothesis of this filter is that low outliers scatter, and the amount of close low outliers is supposed to be smaller than six points, which is true for most data.However, the cluster of low outliers in sample 41 contains more than six points, which result in the remaining low outliers.During filtering with directional structuring element, the ground points between the middle low outliers and the data edge are incorrectly eliminated.Figure 13 shows the two results with and without directional structuring elements.The omission, commission and total errors when filtering without the directional structuring element are 15.42%, 3.00% and 9.20%, respectively.The omission error in the middle is eliminated as shown in Figure 13b.Filtering sample 41 with a cluster of low outliers caused by multi-path error has the largest omission error.A hypothesis of this filter is that low outliers scatter, and the amount of close low outliers is supposed to be smaller than six points, which is true for most data.However, the cluster of low outliers in sample 41 contains more than six points, which result in the remaining low outliers.During filtering with directional structuring element, the ground points between the middle low outliers and the data edge are incorrectly eliminated.Figure 13 shows the two results with and without directional structuring elements.The omission, commission and total errors when filtering without the directional structuring element are 15.42%, 3.00% and 9.20%, respectively.The omission error in the middle is eliminated as shown in Figure 13b.Besides the above comparison between the proposed method and the classic morphology-based algorithm in Reference [44], Table 3 shows the three types of average errors of other algorithms performed against ISPRS datasets.Jahromi et al. [46] employed the filter based on artificial neural networks.Susaki [47] filtered the urban datasets by an adaptive slope filtering method.Zhang and Lin [21] improved the progressive TIN densification filtering method by a point cloud segmentation using smoothness constraint.Mongus and Zalik [40] employed a thin plate spline surface interpolation in hierarchical filtering to increase the robustness for varied landscapes.Chen et al. [48] improved the algorithm of Mongus and Zalik [40] to produce low total error.Pingel et al. [42] developed an improved classic morphology-based filter.As shown in Table 3, our proposed method has comparable performance in comparison to other algorithms.The most important is that our proposed method makes it unnecessary to select different window sizes or determine the maximum window size which is necessary and crucial in other filtering methods such as Pingel et al. [42], Mongus and Zalik [40] and Chen et al. [48].The method proposed by Arefi and Hahn [38] employs morphological grayscale reconstruction based on geodesic dilation for filtering LiDAR point clouds, which also avoid the selection of windows.However, there are some differences between our proposed method and the method by Arefi and Hahn [38] as follows: (1) Arefi and Hahn [38] convert the LiDAR data into an image by the nearest neighbor interpolation which cause the loss of accuracy of original data.We use a simple grid index to reorganize point cloud to utilize the simplicity of regular grid and avoid the accuracy loss of interpolation.
(2) For morphological reconstruction in Arefi and Hahn [38], the input image is used as the mask set, and the marker set is generated by subtracting a offset value from the input image.A sequence of marker sets is obtained by varied offset values to get objects with different heights.The series of offset values need to predefine by the prior knowledge, which are difficult to obtain on changing terrain such as steep slopes.For our proposed method, the mask set and marker set are determined by the analysis of the gradient features of objects and ground.Therefore, our method is adaptive for various objects and changing terrain.
(3) The attached objects such as bridges cannot be eliminated effectively by Arefi and Hahn [38].Our proposed method can remove the attached objects successfully due to the use of geodesic transformations along certain directions.
With the comparison of true DEM and the DEM derived by this filtering algorithm, as shown in Figures 14 and 15, complex objects are effectively eliminated, and the terrain features with drastic changes are preserved.This new filtering algorithm achieves reasonable performance at the rugged outskirt area and at the complicated city area.The error distribution and filtering results of point clouds for AHN Test Site 1 are shown in Figures 16a and 17.The obtained omission, commission, and total errors are 1.44%, 17.72% and 6.00%, respectively.The dense forest at the top left is removed successfully from the steep slopes.The buildings and vegetation are mostly eliminated with exception of a few very low objects.The main omission error occurs at the bottom left where a sharp ridge appears.The gradient changes of the ridge exceed the slope defined by the proposed morphological primitive.Thus, the significant sharp portions of the ridge are eliminated.The error distribution and filtering results of point clouds for AHN Test Site 1 are shown in Figures 16a and 17.The obtained omission, commission, and total errors are 1.44%, 17.72% and 6.00%, respectively.The dense forest at the top left is removed successfully from the steep slopes.The buildings and vegetation are mostly eliminated with exception of a few very low objects.The main omission error occurs at the bottom left where a sharp ridge appears.The gradient changes of the ridge exceed the slope defined by the proposed morphological primitive.Thus, the significant sharp portions of the ridge are eliminated.The error distribution and filtering results of point clouds for AHN Test Site 1 are shown in Figures 16a and 17.The obtained omission, commission, and total errors are 1.44%, 17.72% and 6.00%, respectively.The dense forest at the top left is removed successfully from the steep slopes.The buildings and vegetation are mostly eliminated with exception of a few very low objects.The main omission error occurs at the bottom left where a sharp ridge appears.The gradient changes of the ridge exceed the slope defined by the proposed morphological primitive.Thus, the significant sharp portions of the ridge are eliminated.The error distribution and filtering results of point clouds for AHN Test Site 2 are shown in Figures 16b and 18.The obtained omission, commission, and total errors are 1.28%, 13.63% and 7.25%, respectively.This dataset features the mixed objects, such as buildings with varied shapes and sizes, different vegetation among buildings, and small objects (e.g., cars).Various objects are effectively  The error distribution and filtering results of point clouds for AHN Test Site 2 are shown in Figures 16b and 18.The obtained omission, commission, and total errors are 1.28%, 13.63% and 7.25%, respectively.This dataset features the mixed objects, such as buildings with varied shapes and sizes, different vegetation among buildings, and small objects (e.g., cars).Various objects are effectively The error distribution and filtering results of point clouds for AHN Test Site 2 are shown in Figures 16b and 18.The obtained omission, commission, and total errors are 1.28%, 13.63% and 7.25%, respectively.This dataset features the mixed objects, such as buildings with varied shapes and sizes, different vegetation among buildings, and small objects (e.g., cars).Various objects are effectively filtered as shown by comparing true non-ground points and non-ground points obtained by this filtering approach.In addition, the uneven ground surface is retained well as shown in Figure 18d.filtered as shown by comparing true non-ground points and non-ground points obtained by this filtering approach.In addition, the uneven ground surface is retained well as shown in Figure 18d.

Conclusions
Diverse objects and changing terrains in practical environments pose significant challenges to automatic LiDAR data filtering.The selection of different window sizes and the determination of the maximum window size are crucial problems in many filtering algorithms such as classic morphology-based approaches.Filtering quality is sensitive to the window selection, which imposes considerable limitation in unknown survey areas.To deal with the dependence on the selection of windows, a novel filtering method of LiDAR point clouds is proposed based on the geodesic transformations of mathematical morphology in this paper study.The advantage of geodesic transformations is the avoidance of selecting windows for filtering.The morphological primitive created by combining geodesic transformations has no fixed shape and is self-adaptive for diverse objects.Thus, this algorithm adopts only the elementary window, avoiding the selection of various window sizes, which enhances robustness and reduces the need for a priori scene knowledge.Our experiments indicate that this new filtering method achieves a promising and comparable performance when faced with diverse landscapes, which can effectively preserve the changing terrain features and filter non-ground points in various complicated environments.

Conclusions
Diverse objects and changing terrains in practical environments pose significant challenges to automatic LiDAR data filtering.The selection of different window sizes and the determination of the maximum window size are crucial problems in many filtering algorithms such as classic morphology-based approaches.Filtering quality is sensitive to the window selection, which imposes considerable limitation in unknown survey areas.To deal with the dependence on the selection of windows, a novel filtering method of LiDAR point clouds is proposed based on the geodesic transformations of mathematical morphology in this paper study.The advantage of geodesic transformations is the avoidance of selecting windows for filtering.The morphological primitive created by combining geodesic transformations has no fixed shape and is self-adaptive for diverse objects.Thus, this algorithm adopts only the elementary window, avoiding the selection of various window sizes, which enhances robustness and reduces the need for a priori scene knowledge.Our experiments indicate that this new filtering method achieves a promising and comparable performance when faced with diverse landscapes, which can effectively preserve the changing terrain features and filter non-ground points in various complicated environments.

Figure 1 .
Figure 1.Organizing point clouds with a grid index: (a) input point clouds; (b) constructing a grid index on point clouds; and (c) filling empty grid cells with virtual points.

Figure 1 .
Figure 1.Organizing point clouds with a grid index: (a) input point clouds; (b) constructing a grid index on point clouds; and (c) filling empty grid cells with virtual points.

Figure 3 .
Figure 3. Example of the morphological primitive for filtering: (a) input set f; and (b) the morphological primitive self-adaptively fitting different objects.

Figure 3 .
Figure 3. Example of the morphological primitive for filtering: (a) input set f ; and (b) the morphological primitive self-adaptively fitting different objects.

Figure 4 .Figure 4 .
Figure 4. Example of the morphological filter by combining geodesic transformations: (a) original point set (the red points are the low transition points); (b) the process of reconstruction by erosion (the blue points represent the marker set for morphological reconstruction operator by erosion); (c) (c) Remote Sens. 2017, 9, 1104 7 of 20 the results of reconstruction by erosion (the blue points represent the marker set for morphological reconstruction operator by dilation); and (d) the results of reconstruction by dilation.

Figure 5 .
Figure 5. Example of filtering ground surface of abrupt changes (the blue points are not retrieved to the original ground surface).

Figure 5 .
Figure 5. Example of filtering ground surface of abrupt changes (the blue points are not retrieved to the original ground surface).

Figure 6 .
Figure 6.Example of attached objects: (a) buildings on steep slopes; and (b) bridges connected to ground surface.

Figure 6 .
Figure 6.Example of attached objects: (a) buildings on steep slopes; and (b) bridges connected to ground surface.

Figure 7 .
Figure 7.Comparison of type I errors for filtering the International Society for Photogrammetry and Remote Sensing (ISPRS) datasets.Figure 7. Comparison of type I errors for filtering the International Society for Photogrammetry and Remote Sensing (ISPRS) datasets.

Figure 7 .
Figure 7.Comparison of type I errors for filtering the International Society for Photogrammetry and Remote Sensing (ISPRS) datasets.Figure 7. Comparison of type I errors for filtering the International Society for Photogrammetry and Remote Sensing (ISPRS) datasets.

Figure 7 .
Figure 7.Comparison of type I errors for filtering the International Society for Photogrammetry and Remote Sensing (ISPRS) datasets.

Figure 8 .
Figure 8.Comparison of type II errors for filtering ISPRS datasets.

Figure 9 .
Figure 9.Comparison of total errors for filtering ISPRS datasets.

Figure 8 .
Figure 8.Comparison of type II errors for filtering ISPRS datasets.

Figure 7 .
Figure 7.Comparison of type I errors for filtering the International Society for Photogrammetry and Remote Sensing (ISPRS) datasets.

Figure 8 .
Figure 8.Comparison of type II errors for filtering ISPRS datasets.

Figure 9 .
Figure 9.Comparison of total errors for filtering ISPRS datasets.

Figure 9 .
Figure 9.Comparison of total errors for filtering ISPRS datasets.

Figure 10 .
Figure 10.Comparison of results of filtering Sample 52: (a) results of Li et al. [44]; and (b) results of this proposed method.

Figure 10 .
Figure 10.Comparison of results of filtering Sample 52: (a) results of Li et al. [44]; and (b) results of this proposed method.

Figure 11 .
Figure 11.Error distribution of filtering urban samples 11-42 of ISPRS data.

Figure 11 .
Figure 11.Error distribution of filtering urban samples 11-42 of ISPRS data.

Figure 12 .
Figure 12.Error distribution of filtering rural samples 51-71 of ISPRS data.

Figure 13 .
Figure 13.Results of the proposed method filtering sample 41: (a) filtering with elementary isotropic structuring element and directional structuring elements; and (b) filtering with elementary isotropic structuring element.

Figure 12 .
Figure 12.Error distribution of filtering rural samples 51-71 of ISPRS data.

Figure 12 .
Figure 12.Error distribution of filtering rural samples 51-71 of ISPRS data.

Figure 13 .
Figure 13.Results of the proposed method filtering sample 41: (a) filtering with elementary isotropic structuring element and directional structuring elements; and (b) filtering with elementary isotropic structuring element.

Figure 13 .
Figure 13.Results of the proposed method filtering sample 41: (a) filtering with elementary isotropic structuring element and directional structuring elements; and (b) filtering with elementary isotropic structuring element.

Figure 14 .
Figure 14.Results of the proposed method filtering the Actueel Hoogtebestand Nederland (AHN) Test Site 1: (a) digital surface model (DSM) generated from raw point clouds (the dark line denotes profile position); (b) digital elevation model (DEM) generated from true ground points; (c) DEM generated from ground points derived by the proposed filter; (d) LiDAR points along the profile in (a); (e) LiDAR points along the profile in (b); and (f) LiDAR points along the profile in (c).

Figure 15 .
Figure 15.Results of the proposed method filtering AHN Test Site 2: (a) DSM generated from raw point clouds (the dark line denotes profile position); (b) DEM generated from true ground points; (c) DEM generated from ground points derived by the proposed filter; (d) LiDAR points along the profile in (a); (e) LiDAR points along the profile in (b); and (f) LiDAR points along the profile in (c).

Figure 14 . 20 Figure 14 .
Figure 14.Results of the proposed method filtering the Actueel Hoogtebestand Nederland (AHN) Test Site 1: (a) digital surface model (DSM) generated from raw point clouds (the dark line denotes profile position); (b) digital elevation model (DEM) generated from true ground points; (c) DEM generated from ground points derived by the proposed filter; (d) LiDAR points along the profile in (a); (e) LiDAR points along the profile in (b); and (f) LiDAR points along the profile in (c).

Figure 15 .
Figure 15.Results of the proposed method filtering AHN Test Site 2: (a) DSM generated from raw point clouds (the dark line denotes profile position); (b) DEM generated from true ground points; (c) DEM generated from ground points derived by the proposed filter; (d) LiDAR points along the profile in (a); (e) LiDAR points along the profile in (b); and (f) LiDAR points along the profile in (c).

Figure 15 .
Figure 15.Results of the proposed method filtering AHN Test Site 2: (a) DSM generated from raw point clouds (the dark line denotes profile position); (b) DEM generated from true ground points; (c) DEM generated from ground points derived by the proposed filter; (d) LiDAR points along the profile in (a); (e) LiDAR points along the profile in (b); and (f) LiDAR points along the profile in (c).

Figure 16 .
Figure 16.Error distribution of the proposed method filtering AHN data: (a) result of AHN Test Site 1; and (b) result of AHN Test Site 2.

Figure 17 .
Figure 17.Results of the proposed method filtering AHN Test Site 1: (a) true non-ground points; (b) non-ground points derived by this filter; (c) true ground points; and (d) ground points derived by this filter.

Figure 16 . 20 Figure 16 .
Figure 16.Error distribution of the proposed method filtering AHN data: (a) result of AHN Test Site 1; and (b) result of AHN Test Site 2.

Figure 17 .
Figure 17.Results of the proposed method filtering AHN Test Site 1: (a) true non-ground points; (b) non-ground points derived by this filter; (c) true ground points; and (d) ground points derived by this filter.

Figure 17 .
Figure 17.Results of the proposed method filtering AHN Test Site 1: (a) true non-ground points; (b) non-ground points derived by this filter; (c) true ground points; and (d) ground points derived by this filter.

Figure 18 .
Figure 18.Results of the proposed method filtering AHN Test Site 2: (a) true non-ground points; (b) non-ground points derived by this filter; (c) true ground points; and (d) ground points derived by this filter.

Figure 18 .
Figure 18.Results of the proposed method filtering AHN Test Site 2: (a) true non-ground points; (b) non-ground points derived by this filter; (c) true ground points; and (d) ground points derived by this filter.

Table 1 .
Special features of the dataset selected by ISPRS for filter test.
Region Point Spacing Sample Data No. Special Features Urban 1.0-1.5 m 11, 12 Steep slopes with mixture of buildings and vegetation on hillside, small objects 21, 22, 23, 24 Complex buildings, large buildings, roads with bridges, tunnel, data gaps, ramp, disconnected terrain 31 Densely distributed buildings and vegetation, building with irregular roof, mixture of low and high features, low point

Table 1 .
Special features of the dataset selected by ISPRS for filter test.

Table 3 .
Three types of average errors of different filters against ISPRS samples.