Using object-based analysis to derive surface complexity information for improved filtering of airborne laser scanning data

Abstract

Airborne laser scanning (ALS) is a technique used to obtain Digital Surface Models (DSM) and Digital Terrain Models (DTM) efficiently, and filtering is the key procedure used to derive DTM from point clouds. Generating seed points is an initial step for most filtering algorithms, whereas existing algorithms usually define a regular window size to generate seed points. This may lead to an inadequate density of seed points, and further introduce error type I, especially in steep terrain and forested areas. In this study, we propose the use of object-based analysis to derive surface complexity information from ALS datasets, which can then be used to improve seed point generation. We assume that an area is complex if it is composed of many small objects, with no buildings within the area. Using these assumptions, we propose and implement a new segmentation algorithm based on a grid index, which we call the Edge and Slope Restricted Region Growing (ESRGG) algorithm. Surface complexity information is obtained by statistical analysis of the number of objects derived by segmentation in each area. Then, for complex areas, a smaller window size is defined to generate seed points. Experimental results show that the proposed algorithm could greatly improve the filtering results in complex areas, especially in steep terrain and forested areas.

Keywords

airborne laser scanning / object-based analysis / surface complexity information / filtering algorithm

Introduction

Airborne Laser Scanning (ALS) has been widely used during the past decade ( Ackermann, 1999; Sithole and Vosselman, 2004) as a technique for quick acquisition of Digital Surface Models (DSM) and Digital Terrain Models (DTM). Ground and non-ground reflections are confusingly mixed together in ALS datasets. Features, such as buildings ( Sohn and Dowman, 2007; Vu et al., 2009) and vegetation ( Wagner et al., 2008; Nyström et al., 2014), can be derived from the mixed point clouds. The basic task of processing ALS data is to distinguish bare ground points from non-ground points using filtering algorithms ( Sithole and Vosselman, 2004; Meng et al., 2009; Wang and Glenn, 2009), after which the non-ground points can be further classified into finer types, such as buildings ( Sohn and Dowman, 2007; Vu et al., 2009) and trees ( Wagner et al., 2008; Nyströmet al., 2014).

As a key step to deriving terrain information from point clouds, filtering algorithms are a major focus of many research articles, such as Axelsson (2000), Meng et al. (2009), Mongus and Žalik (2012). Meng et al. (2009) developed a multi-directional filtering algorithm to combine the advantages of directional information and neighborhood information. Mongus and Žalik (2012) implemented a thin plate spline to iteratively interpolate surfaces towards the ground, and top-hat transformation was used to enhance discontinuities caused by surface objects. Yan et al. (2012) proposed an object based analysis (OBA) method for filtering ALS data, and terrain points were obtained by a logic flow method utilizing features of objects derived by segmentation.

In most algorithms, the generation of seed points is a key step ( Axelsson, 2000; Mongus and Žalik, 2012). These existing algorithms typically require a regular window size C defined by the user. The entire area of a given dataset is then divided into non-overlap areas by the same window size. It is assumed that the non-ground objects in each area are smaller than size C, and that the lowest point in each area is therefore a ground point, which is called a seed point. Most of the existing algorithms complete further operations based on these seed points. Additional knowledge about the processing area is typically required when defining a window size for the generation of seed points ( Zhang and Whitman, 2005; Mongus and Žalik, 2012). Error type I also occurs in all existing filtering algorithms, especially in steep terrain ( Sithole and Vosselman, 2004) and forested areas ( Liu, 2008).

In this paper, we aim to derive surface complexity information from ALS datasets using an object based analysis method, and then utilize this information to define the window size for the generation of seed points. We first partition point clouds into non-overlap objects by segmentation, and then surface complexity information is obtained by statistical analysis of the distribution characteristics of the derived objects.

Object based analysis and surface complexity

Object Based Image Analysis (OBIA), which in a geographic context is referred to as Geographic Object Based Image Analysis (GEOBIA), has been developed alongside increased spatial resolution of remote sensing image datasets obtained by satellite sensors, such as IKONOS (launched in 1999), QuickBird (2001), and OrbView (2003) ( Blaschke, 2010). We will use the term ‘OBA’ (Object Based Analysis) in this paper. According to the existing literature, objects can be derived from point clouds by segmentation ( Yan et al., 2012). Level 1 features derived from a single object include area, perimeter, shape, etc. ( Liu et al., 2010). Level 2 features address spatial relations between two objects, such as embeddedness, proximity, and adjacency. Level 3 features are composed of spatial patterns in which more than two objects are involved ( Liu et al., 2008). Existing papers demonstrate that by utilizing features derived from remote sensing datasets and building on appropriate rules, OBA can obtain promising results in filtering ( Yan et al., 2012) and classification ( Ke et al., 2010). In this study, we will utilize level 3 features to identify the surface complexity of point clouds.

If the area of a derived object is small, we define it as a tiny object, and then we define surface complexity as the ratio of the total area covered by tiny objects. For instance, as illustrated in Fig. 1, there are four areas marked by red boxes labeled a, b, c and d. It is clear that area a is covered by dense trees; area b contains buildings; area c is covered by sparse trees, and area d is bare land. We define area a and area c as complex areas because they are covered by small objects, such as trees and other tiny objects. Areas b and d are considered as simple areas since area b is covered by buildings with flat surfaces and area d represents an area of flat terrain.

From the surface complexity definition and illustration in Fig. 1, it is clear that buildings do not exist in complex areas, such as areas a and c, due to their flat or sloping surfaces. Therefore, a smaller window size can be used for seed point generation to filter these areas. Area a is located in steep forested areas. If a smaller window size is used for the generation of seed points to do filtering, Error Type I ( Sithole and Vosselman, 2004) of the filtering result will be improved. This is the main idea of our study.

Method

Two steps are used to obtain surface complexity information from an ALS dataset. First, segmentation is carried out to derive objects based on a grid index, and second, surface complexity information is obtained by the statistical analysis of tiny objects.

Segmentation

We first build a grid index, as described in detail in Yan et al. (2012). A data grid D_M filled with points of local minimum z-values is then obtained, followed by the development of an Edge and Slope Restricted Region Growing (ESRRG) method to derive objects from the data grid D_M.

A Local Average Height Differential Matrix (LAHDM) is defined for the proposed ESRRG method as

(1) \mathrm{LAHDM}(pt)=z_{pt}-{\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \sum z_{p{t}^{\text{'}}}},p{t}^{\text{'}}\in N_{pt},

where pt indicates a point in the data grid D_M, pt^’ denotes a neighboring point of pt, and N_pt is the set of all points neighboring pt. Then, if LAHDM(pt) is greater than the user-defined height threshold, pt is defined as an edge point. In our paper, the height threshold is set to be 1.0 m empirically.

Given a point pt(x_pt, y_pt, z_pt) and one of its neighboring points pt₁(x₁, y₁, z₁), the slope criterion check function of pt and pt₁ ( {\theta }_{pt-p{t}_{\mathrm{1}}} ) is defined as

(2) {\theta }_{pt-p{t}_{\mathrm{1}}}=atan\left({\displaystyle \frac{\left|\Delta z\right|}{\sqrt{{(\Delta x)}^{\mathrm{2}}+{(\Delta y)}^{\mathrm{2}}}}}\right)=atan\left({\displaystyle \frac{|z_{\mathrm{1}}-z_{pt}|}{\sqrt{{(x_{\mathrm{1}}-x_{pt})}^{\mathrm{2}}+{(y_{\mathrm{1}}-y_{pt})}^{\mathrm{2}}}}}\right).

The segmentation procedure starts from the first point in D_M, and each point with no object ID is considered as a candidate point of a new object. If a point is accepted as part of the object, it will be labeled with the ID of the corresponding object, and its neighboring points (four-connected neighborhood) will be put into a queue, regarded as candidate points belonging to the same object. Each object is derived when its queue is empty. The segmentation is then continued on the remaining points that have no object ID.

A matrix M_id is utilized to identify different objects by a series of adjacent points with the same ID. If a cell of the data grid D_M belongs to object o, the corresponding position of M_id will be filled with the ID of object id_o. A hash table is then created to store the object using its ID as the key.

The main procedure is shown in Fig. 2.

In this procedure, function ‘get_lahdm’ is used to derive the LAHDM criteria information as defined in Eq. (1). If an accepted point is an edge point, its neighboring points no longer need to be checked to see whether they should be added to the segmentation queue. In addition, the function ‘check_enter’ is used to check whether a neighboring point pt₁(x₁, y₁, z₁) of an accepted point pt(x_pt, y_pt, z_pt) should be accepted as part of the same object. The check function (f) as

(3) f=\{\begin{array}{l}true\hspace{1em}{\theta }_{pt-p{t}_{\mathrm{1}}}\le {\theta }_{th}\\ falseotherwise\end{array},

where {\theta }_{pt-p{t}_{\mathrm{1}}} indicates the slope between pt₁ and pt, as defined in Eq. (2) and is an appropriate threshold for checking whether pt₁ should be accepted. If {\theta }_{pt-p{t}_{\mathrm{1}}} is smaller than the chosen threshold, pt₁ will be accepted or otherwise refused.

Deriving surface complexity information

Based on the objects obtained by segmentation, the distribution information of these objects can be calculated by statistical analysis of the area in which they are contained. In this paper, we utilize the distribution information to identify the surface complexities of different areas.

The analyzed area is first partitioned into N_R areas with a regular grid C (in Eq. (4)). In most filtering algorithms, the regular grid C is used for seed point generation, and is usually defined as the biggest building size in the analyzed area. In our study, the default grid dimensions are set to 60 m by 60 m. In addition, each area R_ij could be defined by Eq. (5).

(4) N_R=floor\left({\displaystyle \frac{X_{\max }-X_{\min }}{C}}+\mathrm{1.0}\right)*floor\left({\displaystyle \frac{Y_{\max }-Y_{\min }}{C}}+\mathrm{1.0}\right),

(5) R_{ij}=\left\{\begin{array}{l}(x,y)|x\in [X_{\min }+j*C, X_{\min }+(j+\mathrm{1})*C]\wedge \\ \hspace{1em}\hspace{1em}y\in [Y_{\min }+i*C, Y_{\min }+(i+\mathrm{1})*C]\end{array}\right\}.

The surface complexity could be further defined by Eq. (6), where th_area stands for the threshold area of tiny objects. In our paper, th_area is set to be 10.0 m² empirically. The objects’ area rate r_area in each area is defined in Eq. (6).

(6) r_{area}^{ij}={\displaystyle \frac{{\displaystyle \sum are{a}_{ob{j}_k}}}{C^{\mathrm{2}}}},INTERSECT(ob{j}_k,R_{ij})=true\wedge are{a}_{ob{j}_k}\le t{h}_{area}.

A threshold Area_th is defined on the basis of statistical information on tiny objects in order to decide whether or not an area is complex (Eq. (7)). This threshold Area_th was set empirically to 0.2.

(7) f=\{\begin{array}{l}true\hspace{1em}t{h}_{area}^{ij}\le Are{a}_{th}\\ falseotherwise\end{array}.

We have therefore derived surface complexity information using the method described above. As illustrated in Fig. 3(a), the red boxes indicate those areas identified as complex areas. This identification is based on the assumption that if an area is composed of many tiny objects and contains no buildings, it is a complex area. After defining the complex areas, a smaller window size is used for seed point generation, as illustrated in Fig. 3(b). As a result, we utilize the surface complexity information to improve the generation of seed points.

Computing complexity

Assuming that the total number of laser points is N, the complexity for building the grid index is O(N). Points can be effectively indexed with the complexity of O(1) ( Yan et al., 2012). By using the grid index ( Yan et al., 2012), there are W*H points with the lowest Z value of each cell for segmentation (W and H is width and height of the grid respectively). The computational complexity of LAHDM is O(W*H). The ESRRG algorithm is similar to that of typical region-growing algorithms. According to Shih (2010), the computational complexity is O(W*H*lg(W*H)).

Assuming that there are N_obj objects derived after segmentation, the computational complexity for deriving surface complexity information is O(N_obj). N_obj is much smaller than W*H, and therefore, we can conclude that the total time complexity for our algorithm is O(N)+O(W*H)+O(W*H*lg(W*H))+O(N_obj)≈O(W*H*lg(W*H)). The segmentation process is therefore the most time-consuming step.

Experiments

The experiments presented here illustrate the utilization of surface complexity information for improving the filtering result. The first dataset we used to demonstrate our algorithm was from the western section of Salzburg (13°03′E, 47°43′N), Austria (Fig. 1). This area consists of a variety of complex surfaces, including steep hill slopes, buildings, vegetation, etc.

Most of the existing filtering algorithms generate seed points by regular window size. These seed points are usually considered as ground points. The remaining points are processed by iterative comparison with these seed points with different methods and criteria. One example is the use of the adaptive triangulated irregular network filtering algorithm ( Axelsson, 2000), a promising and stable filtering algorithm according to the ISPRS filtering test ( Sithole and Vosselman, 2004). In addition, part of the new algorithm is based on Yan et al. (2012), and thus, we compare our filtering results to those obtained using the methods of both Axelsson (2000) and Yan et al. (2012).

We implemented the adaptive triangulated irregular network filtering algorithm ( Axelsson, 2000), and used window size C = 60 m to generate seed points. We obtained the filtering result as shown in Fig. 4(a). Several complex areas are marked with blue boxes labeled a, b, etc. It is clear that the terrain of the complex parts is suppressed, which means that numerous ground points in these areas are wrongly identified as non-ground points. We utilized the surface complexity method to regenerate seed points in these areas using a smaller window size, usually 20 m in our experiments. As shown in Fig. 4(b), the filtering results over similar areas improved significantly using our method.

Figure 5 shows the profiles of the complex areas marked by the blue boxes, and it becomes clear that the suppressed areas in a1 to a5 are improved in b1 to b5. In d1 and e1, some ground points of the highest parts of the complex areas are wrongly classified as non-ground points, while in d2 and e2 the filtering results are improved significantly using our method. Figure 5 illustrates that the filtering result is greatly improved using our method, especially at the edges of complex terrain areas.

Another two datasets (named ‘site1’ and ‘site2’) used to verify our algorithm are located in the rural region of the western part of Salzburg (13°03′E, 47°43′N), Austria. There are 1,824,683 and 1,109,030 points in the two sites respectively. We obtained the ground points of each site by manual editing, with the result considered as the reference filtering result.

According to Sithole and Vosselman (2004), Error Type I, Error Type II and Error Total information are introduced to verify the performance of a filtering algorithm. Error Type I (eI), Error Type II (eII) and Error Total (eT) are defined in Eq. (8) as follows:

(8) \{\begin{array}{l}eI={\displaystyle \frac{b}{a+b}}\\ eII={\displaystyle \frac{c}{c+d}}\\ eT={\displaystyle \frac{b+c}{a+b+c+d}}\end{array},

where a is the number of ground points that have been correctly identified as ground points; b is the number of ground points that have been incorrectly identified as non-ground points; c is the number of non-ground points that have been incorrectly identified as ground points; and d is the number of non-ground points that have been correctly identified as non-ground points. It is suggested that filtering should aim to minimize Error Type I because Error Type II is easier to edit manually ( Sithole and Vosselman, 2004).

The value of Error Type I, Error Type II, and Error Type Total of the two sites are given in Table 1.We highlight the ‘eI (%)’ and the ‘eT (%)’ column, effectively showing that our proposed method improves Error Type I. Error Type I of Axelsson (2000) is 11.145% and 12.585% for site 1and site2, respectively. The results of Yan et al. (2012) are 12.573% and 15.647% for site1 and site2, respectively, while the results of our method are 7.169% and 7.336%. Error Type Total of Axelsson (2000) is 10.703% and 12.677% for site1 and site2, respectively. The results of Yan et al. (2012) are 11.977% and 15.642% for site1 and site2, respectively, while the results of our method are 7.107% and 7.855%. However, Error Type II of our method seems a little higher than with the algorithm of Axelsson (2000).

Conclusions and future work

We assume an area to be complex if it is composed of many tiny objects, and that there are no buildings in the complex areas. Based on this assumption, a method to derive surface complexity information from LiDAR datasets using OBA is proposed. A grid index is initially built to organize the original point clouds. Objects are then derived by our proposed Edge and Slope Restricted Region Growing (ESRRG) segmentation algorithm based on the grid index. Due to the height difference information, it is easy to do segmentation of LiDAR point clouds with an appropriate data index and neighborhood method. Based on objects derived by segmentation, surface complexity information is obtained on the number of objects in each area by statistical analysis, which is the equivalent of level 3 features in OBA. We then utilize surface complexity information to guide the generation of seed points. The experiments demonstrate that our method could improve the filtering results, especially for Error Type I and Error Type Total in steep terrains and steep forested areas, by as much as 4%‒5% compared to the algorithm of Axelsson (2000).

The surface complexity information can be used to guide further classification of non-ground points, which are usually composed of building points and vegetation points. We assume that buildings should not appear in the complex areas. Based on objects derived by segmentation of non-ground points, buildings can then be extracted by surface complexity information and other restrictions, such as area and average height, as illustrated in Fig. 6. We will focus on this issue in our further studies.

References

[1] Ackermann F (1999). Airborne laser scanning-present status and future expectation. ISPRS J Photogramm Remote Sens, 54(2–3): 64–67

[2] Axelsson P (2000). DEM generation from laser scanner data using adaptive TIN models. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XXXIII(Part B4/1): 110–117

[3] Blaschke T (2010). Object based image analysis for remote sensing. ISPRS J Photogramm Remote Sens, 65(1): 2–16

[4] Ke Y, Quackenbush L J, Im J (2010). Synergistic use of QuickBird multispectral imagery and LIDAR data for object-based forest species classification. Remote Sens Environ, 114(6): 1141–1154

[5] Liu H X, Wang L, Sherman D, Gao Y G, Wu Q S (2010). An object-based conceptual framework and computational method for representing and analyzing coastal morphological changes. Int J Geogr Inf Sci, 24(7): 1015–1041

[6] Liu X (2008). Airborne lidar for dem generation: some critical issues. Prog Phys Geogr, 32(1): 31–49

[7] Liu Y, Guo Q H, Kelly M (2008). A framework of region-based spatial relations for non-overlapping features and its application in object based image analysis. ISPRS J Photogramm Remote Sens, 63(4): 461–475

[8] Meng X L, Wang L, Silván-Cárdenas J L, Currit N (2009). A multidirectional ground filtering algorithm for airborne LIDAR. ISPRS J Photogramm Remote Sens, 64(1): 117–124

[9] Mongus D, Žalik B (2012). Parameter-free ground filtering of LiDAR data for automatic DTM generation. ISPRS J Photogramm Remote Sens, 67(1): 1–12

[10] Nyström M, Holmgren J, Fransson J E S, Olsson H (2014). Detection of windthrown trees using airborne laser scanning. Int J Appl Earth Obs Geoinf, 30(8): 21–29

[11] Shih F Y (2010). Image Processing and Pattern Recognition: Fundamentals and Techniques. New Jersey: Wiley-IEEE Press, 157–158

[12] Sithole G, Vosselman G (2004). Experimental comparison of filter algorithms for bare-Earth extraction from airborne laser scanning point clouds. ISPRS J Photogramm Remote Sens, 59(1–2): 85–101

[13] Sohn G, Dowman I J (2007). Data fusion of high-resolution satellite imagery and LiDAR data for automatic building extraction. ISPRS J Photogramm Remote Sens, 62(1): 43–63

[14] Vu T T, Yamazaki F, Matsuoka M (2009). Multi-scale solution for building extraction from LiDAR and image data. Int J Appl Earth Obs Geoinf, 11(4): 281–289

[15] Wagner W, Hollaus M, Briese C, Ducic V (2008). 3D vegetation mapping using small-footprint full-waveform airborne laser scanners. Int J Remote Sens, 29(5): 1433–1452

[16] Wang C, Glenn N F (2009). Integrating LiDAR intensity and elevation data for terrain characterization in a forested area. IEEE Geosci Remote Sens Lett, 6(3): 463–466 doi:10.1109/LGRS.2009.2016986

[17] Yan M L, Blaschke T, Liu Y, Wu L (2012). An object-based analysis filtering algorithm for airborne laser scanning. Int J Remote Sens, 33(22): 7099–7116

[18] Zhang K, Whitman D (2005). Comparison of three algorithms for filtering airborne lidar data. Photogramm Eng Remote Sensing, 71(3): 313–324

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg