Dual-directional profile filter for digital terrain model generation from airborne laser scanning data

Abstract The most important aspect of digital terrain model generation from airborne laser scanning (ALS) data is that of filtering a point cloud to obtain ground points. Numerous automatic filters have been proposed since ALS data became available. However, to filter out nonground points, a slope threshold is usually introduced to classify points into ground and nonground points; this leads to frequent over-filtering problems in cliff-like terrains. A solution to this problem is proposed, using a dual-directional slope-based filter originating from a conventional slope-based filter is proposed. This filter is designed as a directional filter in one dimension and is applied to every profile of light detection and ranging (LiDAR) points. In this process, a directional filter is first applied to the profile, and another directional filter is then applied at an angle of 180 deg from the first one. Each directional slope-based filter is complementary to the others, thus avoiding over-filtering. We utilize ISPRS LiDAR data for the test. A comparison of this filter approach with existing methods is presented. The comparison result shows that the proposed method obtains a classification accuracy that is as good as most of the compared methods, but is superior to them with regard to handling data from abrupt surfaces.


Introduction
The mapping of topographical surfaces has become extremely efficient and accurate since the availability of airborne laser scanning (ALS) systems. Based on light detection and ranging (LiDAR) technique integrated with positioning and orientation systems, ALS system scans terrains to acquire accurate three-dimensional (3-D) coordinates of scanned points. 1 The resulting densely distributed points are usually called point clouds. ALS data have several applications, including the mapping of corridors, rapid mapping, damage assessment following natural disasters, 2 forest canopy height mapping, 3 and city modeling. 4 The most common process for such applications is the generation of digital terrain models (DTMs). Over the last decade, ALS has gradually replaced photogrammetric techniques to become the new standard for DTM generation for its time and cost effectiveness. 5,6 Especially, the ability of a laser beam to penetrate a tree canopy makes ALS a superior tool for DTM generation in forest areas. 7,8 Given that DTMs accurately map the Earth's surface, a key step in DTM generation is to remove nonground points from ALS point cloud data. This process is frequently called pointcloud filtering. The reserved ground points are used for the interpolation of a grid DTM or to form a triangulated irregular network for a DTM. Several filters are directly applied to discrete raw point clouds, although some filters are designed for gridded data derived and/or interpolated from original point clouds. 9 Filters (e.g., Keqi et al.; 10 Silván-Cárdenas and Wang 11 ) applied to such gridded data typically run fast and can adopt widely available and robust image processing algorithms. 12 However, the derivation of gridded data may cause an information loss of subtle elevation changes. 13 Filters proposed for DTM generation can generally be categorized into three groups based on their approach, namely, linear regression, morphology-based, and slope-based methods. 11,13,14 The linear regression method is also known as an interpolation-based filter. This method iteratively approximates the ground surface using weighted linear regression. 13 The target surface area is first estimated, for example, the average surface between the ground and designated nonground points. The residuals, that is, the oriented distances from the surface to the points, are then calculated. 15 Each point is given a weight according to the residuals. A ground point usually has a negative residual; thus, a high weight will be assigned. Contrarily, a very low or zero weight will be assigned to a point when its residual is positive and large. The surface is finally re-estimated with weight consideration. This process is repeated until a stable situation or a maximum number of iterations is reached. The linear regression algorithm tends to inhibit high-frequency data; therefore, the generated DTMs may be over-smoothed. 14,16 Morphological filters classify points based on morphological operations, such as opening and closing. A critical issue for morphological methods is the selection of a suitable window size. 13,17 On one hand, using a large window will produce a surface with flattened protruding terrain features. 18 On the other hand, a significantly small window may fail to filter out a large group of points, such as the point group of a large building. A solution to this problem is by gradually increasing the window size of the morphological filters and filtering points according to the introduced elevation difference. 10,12 Slope-based filters remove nonground points based on the assumption that the gradients between ground and nonground points are distinctively different. Although slope-based filters preserve sharp corners on steep terrains, 19,20 a steep terrain may still be over-smoothed 21 caused by improperly setting the slope thresholds. To deal with this problem, Sithole 19 proposed the use of an adaptive slope-based (AS) filter that applies an adaptive threshold derived from a trend surface of the terrain. This modified slope filter improves filtering results, but some important topographical features still tend to be over-smoothed, such as steep drop-offs and terraced hillsides. Sithole and Vosselman 21 compared eight different filters and reported that a terrain including steep drop-offs is one of the most difficult terrains for filters and that slope-based filters have the most difficulties with such a precipitous terrain among their compared filters. In addition, a number of one-dimensional (1-D) filters have been designed to run along a scan line or a profile. These filters are called directional filters. Point clouds usually need to be segmented or organized in advance to meet the prerequisite for applying directional filters. However, in spite of this inconvenience, directional filters have the advantages of easy implementation and fast computational processes. The problem of a directional filter is that the final result depends on the direction of the measurement initially chosen. Some significant terrain features may be lost with an ineffective choice of direction. 5 For this reason, bi-directional 20 or multidirectional 5 labeling filters have been proposed. Meng et al. 5 utilized a four-directional filter with slope and elevation difference conditions to filter out nonground points. However, they did not estimate an approximate ground surface, which is critical for a slope-based filter. Current directional filters have yet to be employed for variant surface slopes.
To utilize existing filtering algorithms and solve their inherent problems, this article proposes an improved filtering algorithm. This new method combines the adaptive-slope and directional filtering concepts into a dual-directional slope-based (DS) filter. The DS filter is designed as a 1-D directional filter and is applied to every profile of LiDAR points. In this process, a directional filter is first applied to the profile, and another directional filter is then applied at an angle of 180 deg from the first one. Each directional slope-based filter is complementary to the others, thus avoiding over-filtering. Accordingly, this new algorithm retains sharp topographic features that are over-smoothed by a conventional slope-based filter and upgrades the quality of the generated DTMs. An initial filtering result was first introduced in our previously presented article. 22 In this article, we propose an improvement of a preprocessing step that co-operates with the DS filter for point-cloud filtering. We compare our proposed DS filter with several other filters. The theoretical formulas of the DS filter and the DTM generation are also elucidated.

Slope-Based Filter
Slope-based filtering of ALS data was proposed by Vosselman. 23 A slope filter is similar to a 3-D filtering mask cone (Fig. 1) defined by two parameters, w and s. w is the window size of the mask cone, and s is the slope that defines the cone's opening angle. The cone peak is moved to each testing point and is then used to check whether a neighboring point lies within or below the cone. If such a point is found, the testing point is recognized as a nonground point. Hence, this slopebased filter can be used to check any other points existing under the cone. If the check is positive, then the testing point is labeled as nonground; otherwise, the testing point is considered to be a ground point [see Fig. 1 Terrain slope in the real world naturally varies in different locations; thus, using a fixed slope parameter is impractical. For this reason, Sithole 19 proposed an AS filter that changes the slope threshold subject to the gradient over the local area covered by the filter. Preprocessing is usually needed to obtain an initial approximate ground surface in which the gradient can be estimated.

Dual-Directional Slope-Based Filter
Although the AS filter has enough flexibility to deal with surfaces having various slopes, overfiltering problems in steep terrain areas still remain because this filter has a built-in assumption that the surface slopes do not vary rapidly. Figure 2 shows a schematic showing a possible  over-filtering case for a terrain with very steep drop-offs. Under this circumstance, the interpolated DTM based on the resulting ground points will be over-smoothed because of the loss of some ground points near the most irregular areas [ Fig. 2(c)]. To avoid over-filtering, we propose a DS filter, which is a modified version of the AS filter. The DS filter is a 1-D AS filter but is divided into two equal parts. Each part becomes a directional filter [ Fig. 3(a)]. When applying these two directional filters to a profile, two sets of filtering results are obtained. Figure 3(b) shows a schematic illustrating the application of the DS filter to a cliff-like terrain. The union of the two resulting data sets is the final result of the DS filter [ Fig. 3(c)]. The ground points near the cliff-like terrain are preserved.

Filtering Algorithm of Dual-Directional Slope-Based Filter
The DS filter works along a 1-D space. To apply the DS filter, reorganization of point clouds must be prepared in advance. To this end, we segment point clouds into many strips in both the horizontal (east to west) and vertical (north to south) axes, as shown in Figs. 4(a) and 4(b). The point clouds are arranged into a grid. Each row and column of the grid is considered as a strip. The DS filter is then applied to each strip, which is also considered as a LiDAR profile. Filtering using the DS filter along either the horizontal or vertical axis should generally be sufficient because the cliff topography is kept either in the profile along the horizontal or vertical axis [e.g., Fig. 4(c-i)]. However, a special case will occur when the breaking line is parallel along either the vertical or horizontal axis [Figs. 4(c-ii) and 4(c-iii)]. In such cases, the DS filter can only demonstrate its advantage when it is applied to the profile in a direction perpendicular to the direction of the broken line. To avoid the possibility of ineffective filtering in such unusual cases, both the vertical and horizontal directions are selected in this study to segment point At each strip, the DS filtering algorithm will result in two ground point data sets by the two DS filters with opposite directions. As a result, four sets of filtered points are obtained, that is two sets (Ground 1 , Ground 2 ) from the horizontal strips and two sets (Ground 3 , Ground 4 ) from the vertical strips. The final result Ground Final of the DS filtering algorithm is the union of the filtered ground point data sets, which is expressed as follows:

Digital Terrain Model Generation
The DS algorithm needs a rough ground surface, which is used for a variable slope definition. A number of preprocessing methods (ground seed searching and ground region growing) have been utilized to obtain such a rough ground surface. The preprocessing derives an initial estimated ground surface to obtain the local slope for use by the DS filter. Obvious nonground points (e.g., large buildings and tall trees) are removed in ground region growing. Such a removal can ease the difficulties of filtering when applying the DS filter to point clouds. Figure 5 shows the procedure for the proposed DTM generation. After the ground seed searching and ground region growing, the remaining candidate ground points are then input to the DS filter. The DS filter's resulting output comprises the final classified ground points. The process is thoroughly explained in the following subsections.

Ground Seed Searching
The very low points, which are denoted as noise points, are first removed in this process. The noise points are typically relatively lower than the other near points. Noise points do not usually group together. These noise points can be classified using the conditions that the noise points have relatively low elevation and are isolated. We utilize the TerraScan module for the noise point classification. When the noise points are removed, a conventional slopebased filter 23 is applied to determine the best ground points to use as ground seeds. The conventional slope filter utilizes a fixed cone, with a shape controlled by a slope and its window size, which is applied to all point clouds. The slope must be set <0.3 based on our experiments to ensure that the searched points are ground points. The window size of the cone is the key parameter in this process, which must be set to a value larger than the longest side of objects (e.g., buildings). An urban area typically has large and dense buildings; thus, a wide window size should be set. However, a narrow window size is chosen to include numerous ground seeds. Based on our experiments and experiences, the window size is set at 20 to 70 m in cities and at 10 to 20 m in forest areas.

Ground Region Growing and Dual-Directional Slope-Based Filtering
Once ground seeds are established, ground region growing can be conducted. This process refers to the determination of all possible ground points, in which some nonground points may be included. To implement this process, point clouds are segmented into numerous strips in both the horizontal and vertical directions in advance, as shown in Figs. 4(a) and 4(b). Thus, the neighboring relation between points can be established in each strip either in a horizontal or vertical direction. The region growing begins at each ground seed. If the elevation difference ΔH between the seed point p o and its neighborhood p i is smaller than a threshold (T), then the neighborhood point (p i ) is labeled as a new seed. The process continues until the difference exceeds the threshold or meets another seed. The threshold T is set according to the variations of the topographic relief. A large value of T is used in steep terrains (3 m in our experiments); otherwise, T is set not greater than 2 m. This process can generally remove obvious nonground points with elevations that are very high relative to the ground, that is, points representing such nonground elements as buildings or trees. All the selected ground seeds in this process may be considered as candidate ground points and will be further processed using the DS filters in the subsequent step. The candidate ground points are segmented into numerous LiDAR profiles, as mentioned in Sec. 2.3 and Fig. 4. The DS filtering algorithm is then applied to the generated profiles. Two parameters are needed for performing the DS algorithm. One parameter is the window size of the DS filter, which is set the same as the window size of the conventional filter used in ground seed searching (Sec. 3.1). The other parameter is the adaptive slope, which is derived based on the ground seeds covered by the window of the DS filter. After the DS filtering, the ground points are obtained.

Test Data and Filters
The DS filter was compared with some known filters, including the filter indicated in the commercial software TerraSolid. Test data were downloaded from the ISPRS website (http://www.itc .nl/isprswgIII-3/filtertest/index.html). Fifteen data sets extracted from a variety of terrains (S11 to S42 in urban areas with an average point density range of 1.0 to 1.5 m and S51 to S71 in rural areas with an average point density range of 2.0 to 3.5 m) were tested. Ground truth data (results of manual filtering) and the test results obtained using eight various methods are also provided on the website. We chose all eight filtering results for comparison to our algorithm. Among these methods, Axelsson's algorithm yields better results than the others in most of the 15 data sets in terms of total error. A detailed description of the chosen filters and the test data set can be found in Sithole and Vosselman's study. 21 The classification accuracy of our algorithm was also compared with that of some recently proposed filtering algorithms. We applied the AS filter and TerraSolid module (TerraScan) to the test data. The filtering principle of TerraScan utilizes Axelsson's algorithm. 24 The advantages and shortcomings of the DS filter will be elucidated based on comparison to the results of the AS filter and TerraScan.

Accuracy Assessment
Each point is classified as ground or nonground after filtering. A true positive (tp) occurs if a ground point is correctly classified as ground; a false negative (fn), which is also known as an omission error, occurs if a ground point is classified as nonground; a false positive (fp), which is also known as a commission error, occurs if a nonground point is classified as ground; a true negative (tn) occurs if a nonground point is correctly classified as a nonground. The commission errors usually cause a rough surface in the resulted DTMs and sometimes lead to a significant difference between the generated DTM and the real ground surface when a commission point's elevation is higher than the ground. The omission error usually decreases the details of the ground surface, especially when the omission points represent a specific characteristic of the ground surface. In this study, the classification results are judged in terms of type I error, type II error, total error, and Cohen's kappa. All the assessment indices can be obtained using Eq. (2): where N ¼ tp þ tn þ fp þ fn is the total number of inputs, Po ¼ tpþtn N is the observed agreement, and Pc ¼ ðtpþfpÞÃðtpþfnÞþðfnþtnÞÃðfpþtnÞ is the chance agreement.

Results and Analysis
The employed parameter values of the DS algorithm and the errors in each test site are shown in Table 1. The results indicate that the classification is generally better in urban areas (S11 to S42) than in rural areas (S51 to S71). This result could be due to the complexity of the surface features in the test area. Type I error, type II error, total error, and kappa values for the 15 samples processed using the eight methods in the report and our algorithm are compared in Tables 2 through 5, respectively. Based on the average of the four evaluation indices, the DS algorithm has the best performance in terms of type I error (3.42%), total error (4.32%), and kappa values (84.66%). Type II error is generally inverse to type I error for all the compared algorithms, of which the DS filtering algorithm has the worst result (9.28%). Such a result can be explained by the fact that filtering errors easily occur in the landscape boundary. The DS filter also needs geometric information in multidirections to judge a point, that is, either ground or nonground. When points locate near the landscape boundary, regions outside the test areas where no data were sensed will not provide the geometric information for the DS filtering algorithm to judge the points. Thus, filtering errors, either type I or II, could occur. Figures 6(d) and 6(h) show several example cases in which errors occurred near the boundaries when using the DS filter. For all the compared filters, type II error is generally low (under 10% for each of the 15 data sets). However, type I error results can reveal the filtering ability of the eight methods. Types I and II error results indicate that all the filters are capable of removing more commission errors than omission errors (considering that a low type II error necessarily means few commission errors, and all the filters have low type II error for the 15 data sets). The results also show that the total error correlates well with type I error rather than type II error. When all the filters have similar performances for type II error, type I error becomes the key factor in determining their final classification results, that is, their total error. We also compared the total error and kappa values of the DS filtering algorithm with some other leading edge presented algorithms after 2009. The comparison result shown in Table 6 indicates that our algorithm is as good as the compared algorithms, except Pingel et al.'s algorithm, that achieved the lowest total error and the highest kappa value.
To further test the merit of the proposed DS filter with regard to handling ragged terrain, we applied and compared AS, TerraScan (using Axelsson's algorithm), and the DS filter to study S52 and S53 which contain some ragged terrains. The spatial distributions of types I and II error for the three filters are shown in Fig. 6. The blue points in Fig. 6 indicate type I errors, which are mostly distributed on the ragged terrain. The results reveal that the AS algorithm has considerable type I errors [Figs. 6(b) and 6(f)]. Several of these type I errors are successfully avoided by the DS filter [Figs. 6(d) and 6(h)]. This result demonstrates the effectiveness of the modification from the AS filter to the DS filter by adding directional information. The TerraScan algorithm also performs better than the AS filter. However, the cliff in TerraScan algorithm's results [Figs. 6(c) and 6(g)] erodes more than that in the DS algorithm's results [Figs. 6(d) and 6(h)]. We also investigated the 3-D ground surfaces of S53 using the filtered ground points shown in Fig. 7. The ground surfaces were generated using the reference (a), AS-filtered (b), TerraScan-filtered (c), and DS-filtered ground points (d). The resulting ground surface generated using the DS algorithm can preserve sharp surfaces compared with that generated using AS and TerraScan algorithms. Therefore, the DS filter is more capable of keeping the ground points near the breaking lines rather than TerraScan and AS filters that tend to over-filter such ground points. Identifying numerous possible ground points to preserve the most key features of terrain relief is the ultimate aim of an automatic ALS data filter for DTM generation. Several well-designed filters can filter out most nonground points, but frequently classify some key-feature points  located on ragged terrains into nonground points, which results in over-filtering. The proposed DS filter aims to overcome the over-filtering problem without losing filtering accuracy. This new filtering algorithm uses DS filters in vertical and horizontal axes. Each directional slope-based filter complements the other; thus, the over-filtered ground points can be re-included. In this way, the DS filter keeps the good properties of a slope-based filter and can adapt to rapidly varying topographical reliefs, such as mountains or hills. The experimental results based on the evaluation of total errors and kappa values show a significant improvement in terms of the accuracy and quality of the generated DTMs using the proposed filter. These results are especially promising for extremely uneven terrains. A possible limitation of our proposed DS filter may occur when a dense forest with very few ground points available in the data set is scanned. The DS filter relies on the ground points to  Fig. 7 Three-dimensional ground surface rendering of S53 by using the data from (a) reference ground points (considered as true ground surface), (b) AS ground points, (c) TerraScan ground points, and (d) DS ground points. The resulting ground surface generated using the DS algorithm can preserve sharp surface compared with that generated using AS and TerraScan algorithms (see the breaking lines inside the red polygon).
generate initial ground slope estimation for the use of the DS filter. However, insufficient ground points provide an unreliable slope parameter. The DS filter also needs ground points to discriminate the nonground points. As the number of ground points in a given area decreases, the commission errors increase. However, this limitation is common and is found in most other filters. A possible solution is to increase the overlapping ratio of flying strips when a DTM is wanted for a problematic terrain, such as in a dense forest.