- Email: [email protected]
Extraction of individual tree crowns from airborne LiDAR data in human settlements
Mathematical and Computer Modelling 58 (2013) 524–535
Contents lists available at SciVerse ScienceDirect
Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm
Extraction of individual tree crowns from airborne LiDAR data in human settlements Jiping Liu, Jing Shen ∗ , Rong Zhao, Shenghua Xu Research Center of Government Geographic Information System, Chinese Academy of Surveying and Mapping, Beijing 100830, PR China
article
info
Article history: Received 30 September 2011 Accepted 26 October 2011 Keywords: Airborne LiDAR Tree crowns Slope-based filtering Surface growing Multiple echoes Spoke wheel algorithm
abstract Extraction of individual tree crowns is meaningful for many applications. In this paper, a new method is proposed to extract individual trees from airborne LiDAR point clouds in human settlements. In the process of extraction, an improved slope-based filter is employed to separate the non-ground measurements from the ground measurements, the surface growing algorithm is utilized to segment the point clouds into segments, multiple echo information is used to distinguish the tree points from other types of non-ground measurements, and the spoke wheel algorithm is employed to get the accurate edges of each tree at last. Two datasets are employed to test the above method. Experiments show that our approach is capable of extracting more than 85% trees from the point clouds with accuracy higher than 95%, which suggests the promising applications. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction With the fast development of economy and promotion of industrialization and urbanization, the ecological environment in human settlements has been deteriorating. Therefore, more attention has been paid to the ecological environment in human settlements. Trees contribute important ecosystem services to the surrounding natural community, and they are critical to the wellbeing of both humans and wildlife in the human settlements. Particularly, the tree canopy, which consists of individual tree crowns, represents the interface where most of the fundamental interactions between the vegetation and the atmosphere take place [1]. Accurate and up-to-date information on this dynamic vegetation is of vital importance to the preservation of ecological environments. Currently, most of the statistics about the trees such as areas, tree species and tree species distribution, timber volume, increment of timber volume, and mean tree height, are estimated by measuring sample plots manually in field surveys, thus, tree inventories are rather expensive and tedious [2]. During the past several years, more effort has been put into decreasing costs by developing inventory systems that are based on remote sensing, and great progress has been made. For example, [3] delineate individual tree-crowns with very-high resolution optical imagery. However, the existing methods using optical imagery could not acquire the accurate 3D information of trees. Fortunately, the appearance of LiDAR (light detection and ranging, LiDAR) technology potentially provides a relatively novel tool for generating a unique and comprehensive mathematical description of a tree’s 3D structure [4]. The main advantage of this technique is that it provides dense, detailed and accurate point coverage over both the objects and surfaces. LiDAR is becoming a promising technique for modeling the tree canopy and thus for completing several inventory tasks. Several studies focusing on the measurement of tree canopy structure with laser technology have been published, ranging from terrestrial systems beneath the canopy [5] to airborne systems [6,2] at vertical downward angles and potentially even with space-borne sensors [7]. Furthermore, fusion of LiDAR data and other types of remote sensing data is also popular for tree extraction [8]. However, most of existing researches focus on trees in forested areas, and scarcely on trees in human
∗
Corresponding author. Tel.: +86 10 63880576. E-mail address: [email protected] (J. Shen).
0895-7177/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2011.10.071
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
525
Fig. 1. Structure element of our improved slope-based filter.
settlements. It is no doubt that the trees in human settlements are more closely related to the daily life of human beings. But extraction of trees in human settlements is a more challenging task, because the residential areas are more complex than the forested areas due to diversity of various types of natural objects and artificial objects. Pioneering work about tree extraction from LiDAR data in human settlements is done by Secord and Zakhor [9]. They proposed an approach to detecting trees in registered aerial image and LiDAR data, and SVMs is utilized to classify the segmented images. However, the accuracy of the above method not only greatly depends on the registration of image and LiDAR data, but also needs a lot of training samples, which is not practical for applications. In this paper, a method is proposed to extract individual tree crowns from airborne LiDAR point clouds in human settlements. The method is composed of three steps. First, the original LiDAR point cloud is filtered by our improved slopebased filtering method, in which the points are labeled as ground measurements and non-ground measurements, introduced in Section 2. Second, the non-ground measurements are segmented, and classification is performed based on the size of the segments and the proportion of multiple echoes, introduced in Section 3. Third, the tree points are clustered, and boundary of individual tree crown is derived and delineated by the spoke wheel algorithm, introduced in Section 3. In Section 4, we will conduct experiments and perform evaluation. At last, conclusion and discussion is made in Section 5. 2. Filtering and normalization of non-terrain measurements The goal of this section is to distinguish the non-terrain measurements from the terrain measurements, and to obtain the relative heights of non-terrain measurements from the normalized digital surface model (nDSM). An improved filtering algorithm is proposed to finish this goal. 2.1. Slope-based filtering with aid of topological analysis Digital terrain model (DTM) is useful to detect the tree points and extract the 3D information, and it is usually interpolated by ground/terrain points. Generally, the derivation of DTM is one of the main applications of airborne LiDAR, and the key step in DTM generation is filtering. Filtering is the process of classifying a LiDAR dataset into ground (bare earth) returns and non-ground returns. Many algorithms have been developed for filtering, and the approaches can be categorized into morphological [10,11], surface-based [12,13] and segment-based [14] filters. Among these filters, the morphological filter is very popular. Morphological filter is also known as a slope based filter. Since terrain slope is usually different from the slope seen between the ground and the tops of trees and buildings, this slope difference can be used to separate ground and non-ground measurements from a LIDAR dataset. Vosselman [10] developed a slope-based filter which identifies ground measurements by comparing local slopes between a LIDAR point and its neighbors. Some extensions and variants of morphological filters focus on the shape of structure element. For example, Sithole described an adaptive slope based filter with a rotationally symmetric structure element [15]. By admitting more height differences for the steep areas, the performance of this algorithm is improved in steep areas. In most of the human settlements, the terrain slopes are not very steep. Thus, the slope-based method is probably feasible after some improvements. There are three key parameters, as shown in Fig. 1, in the slope-based method:
• the height difference at same location, h1 ; • a filter radius, r; • the height difference at a filter radius, h2 . However, Vosselman’s method has an obvious disadvantage when employed for filtering points in residential areas where the sizes of buildings greatly vary. As a result, it could not filter the large buildings when filter radius r is smaller and it would filter more ground points when the filter radius r is larger. To increase the accuracy of filtering, we improved the original filter by topological analysis when a smaller r is employed. The method implemented here is as follows. (1) Overlay a rectangular mesh on the LiDAR data set. Create a 2D array, whose elements represent points falling in cells of a mesh overlaying the data set. Each point measurement pj (xj , yj , zj ) from the LiDAR data set is assigned into a cell in terms
526
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
of its x and y coordinates. If more than one point falls in the same cell, the one with the lowest elevation is selected as the array element. (2) A LiDAR survey point, p0(x0 ,y0 ,z0 ) , is classified as a ground measurement if the maximum value of slopes between this point and any other point pj within a given radius r is less than predefined threshold: s0,j =
z0 − zj − h1
x0 − xj
2
(1)
+ (y0 − yj )2
p0 ∈ ground measurements if s0,j < arctan ( h1 ), where s0,j is slope between p0 and pj (xj , yj , zj ), xj and yj represent horizontal 2 coordinates of pj (xj , yj , zj ) and zj is its elevation. (3) Topological Analysis. Once the above filtering process is finished, make connected components analysis (namely, connected component labeling). Then, the neighboring points with the same filter state are clustered as one same region. If the regions with filter state as ground meet two conditions, contained by the regions with filter state as non-ground and elevation differences between them and their adjacent non-ground regions are larger than h1 , then the regions’ filter state is relabeled as non-ground. Thus, the inner parts of large buildings can be filtered. The original filter result and our filter result are shown in Fig. 2(c) and (d), respectively, and both the filters have same filter parameters. Parameter values are set as h1 = 0.1, r = 20, h2 = 30. As the figures show, the large buildings are filtered for our proposed method while only parts of large buildings are filtered for Vosselman’s filter. h
2.2. Normalization of non-terrain measurements Once the ground measurements are acquired, various types of interpolation methods [16] such as inverse distance weighting, Kriging and triangulation with linear interpolation can be employed to derive the DTM from the ground points. In this paper, method of triangulation with linear interpolation is utilized. For any non-terrain points, the height of the point on nDSM is equal to the height of the point on DSM subtracted from the height of the point on DTM, as shown in Fig. 3. Thus, the above operations can derive the relative heights of trees to the grounds. 3. Segmentation and detection of tree points Generally, the non-terrain measurements consist of tree points, building points, wall points, bridge points, vehicle points, power line points etc. Among these non-terrain measurements, the tree points have two distinguished characteristics: the first one is the points belonging to the same tree or the points belonging to same population are scarcely coplanar; the second one is most of the tree points belong to the multiple echoes. Therefore, two operations, including automatic segmentation of LiDAR data into coplanar point clusters using a surface growing algorithm and multiple echo analysis, are employed to separate the trees points from other types of non-terrain points. 3.1. Segmentation by the surface growing algorithm Segmentation is the process of labeling each measurement in a point cloud, so that the points belonging to the same surface or region are given the same label. A number of segmentation algorithms for point cloud are available, and a literature survey of different segmentation methods has been presented in [17]. Among these approaches, we adopted the surface growing algorithm proposed by Vosselman et al. [18]. We make a short explanation of the surface growing algorithm here because of its strong relevance to our feature recognition method. The algorithm consists of the following steps. (1) Determine a seed surface. A seed surface consists of a group of nearby points that fit well to a flat plane. The algorithm selects an arbitrary unclassified point and tests if a minimum number of nearby points can be fitted to a plane. If this is the case, these points constitute the seed surface. Otherwise, another arbitrary point is tested. The key of this step is to test if a minimum number of nearby points can be well fitted to a plane, and 3D Hough transformation [18] is employed to finish this task. The principle of 3D Hough transformation is illustrated as follows. A plane in a 3D space can be described by Vosselman and Klein [19]: d = X cos α cos β + Y sin α cos β + Z sin β.
(2)
Every point p(x, y, z ) in a point cloud defines a surface in the 3D parameter space spanned by the plane parameters α , β and d. Points in the point cloud data that are coplanar will correspond to surfaces in the parameter space that intersect at a point. The α , β and d coordinates of this intersection point yield the parameter of the plane in the point cloud. Thus, the above 3D Hough transformation can be used to check if a minimum number of nearby points can be well fitted to a plane, as shown in Fig. 4(a) and (b). (2) Grow the seed surfaces. The growing can be based on one or more the following criteria.
• Proximity of points. Only points within certain distance to a seed surface, or the neighboring points (in kd-tree [20]) of a seed surface, can be added to this seed surface.
• Globally planar. For this criterion, a plane equation is determined by fitting a plane through all surface points in this seed surface. Points can only be added if the perpendicular distance to the plane is below some threshold.
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) Top view of data colored by elevation.
(c) Result of slope-based filter.
527
(b) Perspective view of the data in (a).
(d) Result of our improved slope-based filter.
(e) The non-terrain measurements derived by our method colored by their elevations. Fig. 2. Comparison of the results of the slope-based filter and our improved filter, respectively. Note the ground measurements are colored as green, and the non-ground measurements are colored as red in (c) and (d). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 5(a) shows the original laser scanned point cloud, and Fig. 5(b) shows the segmentation results of this point cloud by the surface growing algorithm. In the ideal situation, each segment should represent one of the semantic features such as roofs, walls, road segments etc. Fig. 5 also suggests that the sizes of segments can reflect the types of objects. Particularly, the sizes of segments belonging to trees, power lines and vehicles are usually smaller than those belonging to buildings. Thus, the sizes of segments can be selected as a feature to classify the point clouds.
528
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
Fig. 3. Derivation of the height of a non-terrain point from DSM and DTM.
(a) Three discrete coplanar points.
(b) 3D Hough space where the surfaces, corresponding to the above there points, intersect at a point.
Fig. 4. Voting in the 3D Hough space for plane detection.
3.2. Multiple returns analysis Katzenbeisser [21] has mentioned that one laser pulse can produce multiple reflections if the light traveling to ground illuminates several objects at different distances. Multiple reflections can occur when:
• • • •
objects are at different elevations; the laser beam has a chance to partly illuminate other objects; objects are large enough and have a sufficient reflectivity to generate echoes; the laser beam is large enough (in diameter) to illuminate several objects.
Fig. 6 shows the difference between tree points and building points. On the building areas, most of the returns are reflected as only one detectable echo, i.e., first echo and last echo are identical; also in Fig. 6(a) the single echo is reflected from low and flat vegetation. In Fig. 6(b), most of the echoes on the top of canopy are from first pulse, and then penetrated through the canopies. There are a number of surface shapes which will always provide multiple echoes, i.e., leafless deciduous trees, edges of buildings or of high trees. We now deal with the challenge of identifying tree regions by using multiple echo information associated with each segment. Laser beams used in our data collection will have a circular footprint with a radius of about 15–30 cm when they hit the surface. If more than one object of high reflectance lies within this footprint, multiple return height information will be received by the sensor. It is therefore expected that tree segments would have a higher number of multiple returns associated with them than the more planar building segments. This multiple return information is inherently point based. However, now that we have region information, we can calculate the percent of data points for each area that recorded a multiple return. By eliminating all regions where the ratio of points with multiple returns to single returns is smaller, we
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) Perspective view of the segmented point cloud in Fig. 3(e), and the points are colored by segment number.
(b) Top view of the segmented point cloud in (a) colored by segment number.
(c) The subset of singular returns colored by segment number.
(d) The subset of non-last returns colored by segment number.
529
(e) Detected tree points colored by elevation. Fig. 5. Segmentation, multiple echo analysis and classification of a point cloud. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
identify a large number of tree regions. In this work, we chose to retain segments with more than 20% multiple returns. However, this number is again dependent on the sampling technology.
530
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) Returns information of a building.
(b) Returns information of trees. Fig. 6. Differences of returns between buildings and trees.
4. Single tree detection and description Airborne LiDAR datasets with point densities of 1 point/m2 or denser offer the possibility of detecting, segmenting and modeling individual tree crowns. Single tree detection and 3D reconstruction allow counting of trees in a certain area and determination of parameters such as individual tree height, crown diameter, crow area and 3D crown shape. Researchers have made great progress. For example, Koch et al. [2] presents an approach to delineate single trees automatically in LiDAR data in deciduous and mixed temperate forests. First, in rasterized laser data, possible tree tops are detected with a local maximum filter. Afterward, the crowns are delineated with a combination of a pouring algorithm, knowledge-based assumptions on the shape of trees, and a final detection of the crown-edges by searching vectors starting from the trees’ tops. However, the above classical method could not be used in human settlements due to the various sizes and heights of the trees in human settlements. In this paper, we present a new approach to detect single trees in human settlements. Our method includes clustering, detecting local maximum with elevation, and identifying the points belonging to one tree. 4.1. Clustering In human settlements, the spatial distributions and configurations of trees vary very much. For example, the crowns of several trees may touch each other and form a continuous canopy, while some trees may be isolated from the others. Thus, we should first make a clustering analysis. In detail, construct the Delaunay triangulated irregular network (D-TIN) of tree points only using the horizontal coordinates, and then delete the long edges in the D-TIN, as shown in Fig. 7(a). Subsequently, make a connected component analysis using the above derived D-TIN data. As a result, the adjacent tree points are labeled as one cluster, and the points belonging to one isolated tree are also labeled as one cluster, as shown in Fig. 7(b). Once the clusters are determined, only the points belonging to the first pulse are remained for subsequent operations, because only the first pulse may come from the surface of a tree and the other types of pulse are sure to not come from the surface of a tree. However, we should note that the other types of pulse are quite meaningful to the clustering analysis. 4.2. Finding tree tops Suppose that the tree tops are the points whose heights are local maximums and the sizes of crowns of trees are only relevant to mean heights of trees. Another two facts are that the size of a tree crown is strongly correlated with the height of the tree belonging to a specific species, and the trees in one cluster generally belong to only one species. Thus, we first
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) Delaunay TIN of the tree points.
(c) Three tree tops are detected from data in (b).
531
(b) Two clusters are determined by connected component analysis, and only first pulses are retained.
(d) Extraction of edges of trees.
Fig. 7. Clustering, detecting tree tops, extracting edges of trees.
calculate the mean height of each cluster, and then derive a searching radius, Rtop , according to the mean height for each cluster. Subsequently, we will build a 2D kd-tree [20] for the points of each cluster only using the horizontal coordinates, and the kd-tree is capable of fast and exact searching the nearest neighbors to a fixed radius bound. Furthermore, if horizontal distance of any two adjacent tree tops is smaller than a preset value, then one of the tree tops must be deleted, which avoids detection of tree tops that belong to one same tree. As a result, we will find the points whose heights are larger than any of their nearest neighbors to a fixed radius. 4.3. Extracting tree edges Theoretically, each tree top has a local sphere of influence, which is anisotropic and directionally circular. That is, along most directions, the branches of the sphere of influence are approximately circular. Hu et al. [22] proposed the spoke wheel operator to delineate the boundary of the sphere of influence. Hu’s spoke and the spoke wheel SW are used for detecting a footprint of a pixel p, as shown in Fig. 8(a). A spoke is a line segment with a length of m (a preset value) pixels, and m = 90 pixels in Fig. 8(a). A spoke wheel is a sequence of spokes S (ϕi , m) (i = 0, . . . , 4n − 1) with common initial point p and evenly spaced angles ϕi = π i/2n, and n = 18 in Fig. 8(a). The set of pixels in a spoke wheel SW centered at the pixel p with 4n spokes is denoted by SW (p, n, m). The intersection between a spoke and a tree edge provides useful information to determine the local sphere of influence around a pixel. However, we do not know the exact location of a tree edge in advance. To search for the intersection of a spoke and the edge of a road, we start from p, move in the direction of the spoke, and observe the absolute elevation differences between p and the pixels along the spoke. The differences are small when the pixels are near to p; however, they may become larger when the pixels are far away from p. Let Si be the spoke at pixel
532
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) The configuration of a spoke wheel.
(b) Extracted tree edges by a spoke wheel algorithm operated on a tree image.
Fig. 8. Spoke wheel algorithm and its application on extraction of tree edges.
(a) Perspective view of pseudo-grid.
(b) Top view of pseudo-grid. Fig. 9. Pseudo-grid generation.
p. The cutting point, denoted by Ci , on Si is the first pixel such that
|I (Cj ) − I (p)| ≥ σi (j = 0, . . . , m, 0 ≤ i < 4n)
(3)
where σi is the intensity standard deviation on SW (p, n, m). Therefore, the threshold in Eq. (3) is adaptive. Subsequently, we connect the cutting points on all spokes around a pixel p in a counterclockwise direction, which results in a closed polygon. This represents the footprint of the pixel p, denoted by F (p). For each detected tree top, a temporary pseudo-grid is generated from the point cloud around the tree top as follows. (1) Suppose that a tree top and its preset searching radius Rtop are already known, clip out the points whose horizontal distance between them and the tree top is less than Rtop . (2) The size of pseudo-grid is calculated with average point density of laser point data. Once pseudo-grid is created, the raw laser point data is assigned to each pseudo-grid shown in Fig. 9. If there is more than one point located in the same grid, only the highest one is retained. (3) If there is no point in a grid, the central point of the grid is added into the grid, and the elevation of the central point is assigned to one of the nearest points in the point cloud. Once the pseudo-grid is generated, the above spoke wheel algorithm is performed on the pseudo-grid image, and tree edges are extracted, as shown in Figs. 7(d) and 8(b). 5. Experiments and performance evaluation In this section, two datasets are employed to extract individual tree crowns, and the performance of our proposed algorithm is evaluated.
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) Point cloud of Enschede colored by elevation.
533
(b) The 396 extracted trees by the proposed method.
(c) The extracted edges of trees are in the same color as the trees in (b). Fig. 10. Extraction of individual tree crowns from data in Enschede. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
5.1. Testing data Two point clouds are used to test our proposed approach. The first data is over the town of Enschede, the Netherlands with 700 m width and 1100 m length, as shown in Fig. 10(a). The dataset was recorded with the FLI-MAP 400 system with 20 points/m2 , and maximum number of returns per pulse reaches to 4. In this first site, there are 423 trees by manual interpretation. The second data is over the suburban of Nanning, Guangxi Province, China with 460 m width and 1000 m length, as shown in Fig. 11(a). The dataset was acquired with the RIEGL LMS-Q560 system with 10 points/m2 , and maximum number of returns per pulse reaches to 4. In this second site, there are 96 trees by manual interpretation. 5.2. Results and performance evaluation The extracted trees and their edges of the first dataset are shown in Fig. 10(b) and (c), respectively. The number of extracted trees is 397 in this dataset, and visual inspection suggests that there are 30 mistakes in this test. Among the mistakes, most of them mix more than one tree into a tree, as shown in Fig. 12(a), and some of them are the powerlines. As far as the second dataset is concerned, the extracted trees and their edges are shown in Fig. 11(b) and (c), respectively. There are 87 trees being detected, and visual inspection suggests that there are 5 false detections in this test. All the false detections mistake a cluster with more than one trees as one isolated tree. The two tests also tell us that if the number of points of one tree is not large enough, our algorithm will miss them. The two tests suggest that the performance of our approach is excellent. As far as the accuracy is concerned, 92.44% and 94.25% are reached for the first test and the second test, respectively. As far as the completeness is concerned, 86.76% and 85.41% are acquired for the above two tests, respectively. The statistics tell that our approach has a promising application in practice. However, our method also has some disadvantages. First, the proposed method will omit the small trees if they only have several points in the dataset. Second, our method will mistake adjacent trees as one tree if the searching radius is not well preset. As last, our method will mistake some other objects as trees if the objects also produce multiple echoes. 6. Conclusions and discussion Extraction of individual tree crowns from airborne LiDAR point clouds is very meaningful for many applications. However, most of the existing methods are oriented to forested areas, and they are not fitful to extraction of the trees from point clouds
534
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) Point cloud of Nanning colored by elevation.
(b) The 87 extracted trees by the proposed method.
(c) The extracted edges of trees are in the same color as the trees in (b). Fig. 11. Extraction of individual tree crowns from data in Naning, China. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(a) Side view of points belonging to two adjacent trees but mistaken as one.
(b) Some power line points are mistaken as trees.
Fig. 12. False detection in extraction of trees.
in human settlements. In these areas, the scenes are more complex than the forested areas. Thus, a method is proposed to extract individual tree crowns from airborne LiDAR point clouds in human settlements. In our method, an improved slopebased filter is employed to separate the non-ground measurements from the ground measurements, the surface growing algorithm is utilized to segment the point clouds into segments, the multiple echo analysis is used to distinguish the tree points from other types of non-ground measurements, and the spoke wheel algorithm is employed to get the edges of each tree. Two datasets are used to test our proposed method. Experiments show that our approach is capable of extracting more than 85% trees from the point clouds with the accuracy no less than 95%, which suggest its promising applications in engineering. However, the proposed method fails to extract the small trees, and will mistake mixed trees as one isolated tree. Future work will focus on the improvement of the proposed method. Acknowledgments This research was funded by the Basic Research Fund of Chinese Academy of Surveying and Mapping under Grant Nos 7771025, 7771109 and 7771017, the Project for Young Scientist Fund sponsored by the Natural Science Foundations of
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
535
China under Grant Nos 41001280 and 40901195, and the China Postdoctoral Science Foundation under Grant 2010047038, respectively. References [1] G.G. Parker, Structure and microclimate of forest canopies, in: M. Lowman, N. Nadkarni (Eds.), Forest Canopies, Academic Press, California, USA, 1995, pp. 73–106. [2] B. Koch, U. Heyder, H. Weinacker, Detection of individual tree crowns in airborne LiDAR data, Photogrammetric Engineering & Remote Sensing 72 (4) (2006) 357–363. [3] L. Wang, A multi-scale approach for delineating individual tree-crowns with very-high resolution imagery, Photogrammetric Engineering & Remote Sensing 76 (4) (2010) 371–378. [4] D.V. der Zande, W. Hoet, I. Jonckheere, J. van Aardt, P. Coppin, Influence of measurement set-up of ground-based LiDAR for derivation of tree structure, Agricultural and Forest Meteorology 141 (2006) 147–160. [5] F.M. Danson, D. Hetherington, F. Morsdorf, B. Koetz, B. Allgöwer, Forest canopy gap fraction from terrestrial laser scanning, IEEE Geoscience and Remote Sensing Letters 4 (1) (2007) 157–160. [6] J. Hyyppä, M. Schardt, H. Haggrén, B. Koch, U. Lohr, H.U. Scherrer, R. Paananen, H. Luukkonen, M. Ziegler, H. Hyyppä, U. Pyysalo, H. Friedländer, J. Uuttera, S. Wagner, M. Inkinen, A. Wimmer, A. Kukko, E. Ahokas, M. Karjalainen, HIGH-SCAN: the first European-wide attempt to derive singletree information from laserscanner data, The Photogrammetric Journal of Finland 17 (2001) 58–68. [7] H. Duong, R. Lindenbergh, N. Pfeifer, G. Vosselman, ICESat full waveform altimetry compared to airborne laser scanning over the Netherlands, IEEE Transactions on Geoscience and Remote Sensing 47 (10) (2009) 3365–3378. [8] T.L. Erdody, L.M. Moskal, Fusion of LiDAR and imagery for estimating forest canopy fuels, Remote Sensing of Environment 114 (4) (2010) 725–737. [9] J. Secord, A. Zakhor, Tree detection in urban regions using aerial LiDAR and image data, IEEE Geoscience and Remote Sensing Letters 4 (2) (2007) 196–200. [10] G. Vosselman, Slope based filtering of laser altimetry data, International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences 33 (Part B3/2) (2000) 935–942. [11] K. Zhang, S. Chen, D. Whitman, M. Shyu, J. Yan, C. Zhang, A progressive morphological filter for removing non-ground measurements from airborne LIDAR data, IEEE Transactions on Geoscience and Remote Sensing 41 (4) (2003) 872–882. [12] P. Axelsson, DEM generation from laser scanner data using adaptive tin models, International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences 33 (Part B4) (2000) 110–117. [13] K. Kraus, N. Pfeifer, Determination of terrain models in wooded areas with airborne laser scanner data, ISPRS Journal of Photogrammetry and Remote Sensing 53 (4) (1998) 193–203. [14] G. Sithole, Segmentation and classification of airborne Laser scanner data, Delft, ISBN 9061322928, p. 184, Publications on Geodesy of the Netherlands Commission of Geodesy, vol. 59, May 24, 2005. [15] G. Sithole, Filtering of laser altimetry data using a slope adaptive filter, In: IAPRS 2001, vol. XXXIV, 3/W4, Annapolis (MD), USA, 2001, pp. 203–210. [16] M.J. De Smith, M.F. Goodchild, P.A. Longley, Geospatial Analysis: A Comprehensive Guide to Principles, Techniques and Software Tools, second ed., Troubador Publishing Ltd., 2007. [17] M. Wang, Y.H. Tseng, Automatic segmentation of LiDAR data into coplanner point clusters using an octree-based split-and-merge algorithm, Photogrammetric Engineering and Remote Sensing 76 (4) (2010) 407–420. [18] G. Vosselman, B.G.H. Gorte, G. Sithole, T. Rabbani, Recognizing structure in laser scanning point clouds, International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences 46 (8/W2) (2004) 33–38. [19] G. Vosselman, R. Klein, Visualisation and structuring of point clouds, in: G. Vosselman, H.-G. Maas (Eds.), Airborne and Terrestrial Laser Scanning, Whittles Publising, ISBN: 978-1904445-87-6, 2010, pp. 43–79. [20] S. Arya, D.M. Mount, N.S. Netanyahu, R. Silverman, A.Y. Wu, An optimal algorithm for approximate nearest neighbor searching in fixed dimensions, Journal of the ACM 45 (1998) 891–923. [21] R. Katzenbeisser, Echo detection, Technica Note TopoSys. Toposys homepage, www.toposys.com, 2003 (accessed 07.06.11). [22] J. Hu, A. Razdan, J.C. Femiani, Road network extraction and intersection detection from aerial images by tracking road footprints, IEEE Transaction on Geoscience and Remote Sensing 45 (12) (2007) 4144–4157.
Contents lists available at SciVerse ScienceDirect
Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm
Extraction of individual tree crowns from airborne LiDAR data in human settlements Jiping Liu, Jing Shen ∗ , Rong Zhao, Shenghua Xu Research Center of Government Geographic Information System, Chinese Academy of Surveying and Mapping, Beijing 100830, PR China
article
info
Article history: Received 30 September 2011 Accepted 26 October 2011 Keywords: Airborne LiDAR Tree crowns Slope-based filtering Surface growing Multiple echoes Spoke wheel algorithm
abstract Extraction of individual tree crowns is meaningful for many applications. In this paper, a new method is proposed to extract individual trees from airborne LiDAR point clouds in human settlements. In the process of extraction, an improved slope-based filter is employed to separate the non-ground measurements from the ground measurements, the surface growing algorithm is utilized to segment the point clouds into segments, multiple echo information is used to distinguish the tree points from other types of non-ground measurements, and the spoke wheel algorithm is employed to get the accurate edges of each tree at last. Two datasets are employed to test the above method. Experiments show that our approach is capable of extracting more than 85% trees from the point clouds with accuracy higher than 95%, which suggests the promising applications. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction With the fast development of economy and promotion of industrialization and urbanization, the ecological environment in human settlements has been deteriorating. Therefore, more attention has been paid to the ecological environment in human settlements. Trees contribute important ecosystem services to the surrounding natural community, and they are critical to the wellbeing of both humans and wildlife in the human settlements. Particularly, the tree canopy, which consists of individual tree crowns, represents the interface where most of the fundamental interactions between the vegetation and the atmosphere take place [1]. Accurate and up-to-date information on this dynamic vegetation is of vital importance to the preservation of ecological environments. Currently, most of the statistics about the trees such as areas, tree species and tree species distribution, timber volume, increment of timber volume, and mean tree height, are estimated by measuring sample plots manually in field surveys, thus, tree inventories are rather expensive and tedious [2]. During the past several years, more effort has been put into decreasing costs by developing inventory systems that are based on remote sensing, and great progress has been made. For example, [3] delineate individual tree-crowns with very-high resolution optical imagery. However, the existing methods using optical imagery could not acquire the accurate 3D information of trees. Fortunately, the appearance of LiDAR (light detection and ranging, LiDAR) technology potentially provides a relatively novel tool for generating a unique and comprehensive mathematical description of a tree’s 3D structure [4]. The main advantage of this technique is that it provides dense, detailed and accurate point coverage over both the objects and surfaces. LiDAR is becoming a promising technique for modeling the tree canopy and thus for completing several inventory tasks. Several studies focusing on the measurement of tree canopy structure with laser technology have been published, ranging from terrestrial systems beneath the canopy [5] to airborne systems [6,2] at vertical downward angles and potentially even with space-borne sensors [7]. Furthermore, fusion of LiDAR data and other types of remote sensing data is also popular for tree extraction [8]. However, most of existing researches focus on trees in forested areas, and scarcely on trees in human
∗
Corresponding author. Tel.: +86 10 63880576. E-mail address: [email protected] (J. Shen).
0895-7177/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2011.10.071
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
525
Fig. 1. Structure element of our improved slope-based filter.
settlements. It is no doubt that the trees in human settlements are more closely related to the daily life of human beings. But extraction of trees in human settlements is a more challenging task, because the residential areas are more complex than the forested areas due to diversity of various types of natural objects and artificial objects. Pioneering work about tree extraction from LiDAR data in human settlements is done by Secord and Zakhor [9]. They proposed an approach to detecting trees in registered aerial image and LiDAR data, and SVMs is utilized to classify the segmented images. However, the accuracy of the above method not only greatly depends on the registration of image and LiDAR data, but also needs a lot of training samples, which is not practical for applications. In this paper, a method is proposed to extract individual tree crowns from airborne LiDAR point clouds in human settlements. The method is composed of three steps. First, the original LiDAR point cloud is filtered by our improved slopebased filtering method, in which the points are labeled as ground measurements and non-ground measurements, introduced in Section 2. Second, the non-ground measurements are segmented, and classification is performed based on the size of the segments and the proportion of multiple echoes, introduced in Section 3. Third, the tree points are clustered, and boundary of individual tree crown is derived and delineated by the spoke wheel algorithm, introduced in Section 3. In Section 4, we will conduct experiments and perform evaluation. At last, conclusion and discussion is made in Section 5. 2. Filtering and normalization of non-terrain measurements The goal of this section is to distinguish the non-terrain measurements from the terrain measurements, and to obtain the relative heights of non-terrain measurements from the normalized digital surface model (nDSM). An improved filtering algorithm is proposed to finish this goal. 2.1. Slope-based filtering with aid of topological analysis Digital terrain model (DTM) is useful to detect the tree points and extract the 3D information, and it is usually interpolated by ground/terrain points. Generally, the derivation of DTM is one of the main applications of airborne LiDAR, and the key step in DTM generation is filtering. Filtering is the process of classifying a LiDAR dataset into ground (bare earth) returns and non-ground returns. Many algorithms have been developed for filtering, and the approaches can be categorized into morphological [10,11], surface-based [12,13] and segment-based [14] filters. Among these filters, the morphological filter is very popular. Morphological filter is also known as a slope based filter. Since terrain slope is usually different from the slope seen between the ground and the tops of trees and buildings, this slope difference can be used to separate ground and non-ground measurements from a LIDAR dataset. Vosselman [10] developed a slope-based filter which identifies ground measurements by comparing local slopes between a LIDAR point and its neighbors. Some extensions and variants of morphological filters focus on the shape of structure element. For example, Sithole described an adaptive slope based filter with a rotationally symmetric structure element [15]. By admitting more height differences for the steep areas, the performance of this algorithm is improved in steep areas. In most of the human settlements, the terrain slopes are not very steep. Thus, the slope-based method is probably feasible after some improvements. There are three key parameters, as shown in Fig. 1, in the slope-based method:
• the height difference at same location, h1 ; • a filter radius, r; • the height difference at a filter radius, h2 . However, Vosselman’s method has an obvious disadvantage when employed for filtering points in residential areas where the sizes of buildings greatly vary. As a result, it could not filter the large buildings when filter radius r is smaller and it would filter more ground points when the filter radius r is larger. To increase the accuracy of filtering, we improved the original filter by topological analysis when a smaller r is employed. The method implemented here is as follows. (1) Overlay a rectangular mesh on the LiDAR data set. Create a 2D array, whose elements represent points falling in cells of a mesh overlaying the data set. Each point measurement pj (xj , yj , zj ) from the LiDAR data set is assigned into a cell in terms
526
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
of its x and y coordinates. If more than one point falls in the same cell, the one with the lowest elevation is selected as the array element. (2) A LiDAR survey point, p0(x0 ,y0 ,z0 ) , is classified as a ground measurement if the maximum value of slopes between this point and any other point pj within a given radius r is less than predefined threshold: s0,j =
z0 − zj − h1
x0 − xj
2
(1)
+ (y0 − yj )2
p0 ∈ ground measurements if s0,j < arctan ( h1 ), where s0,j is slope between p0 and pj (xj , yj , zj ), xj and yj represent horizontal 2 coordinates of pj (xj , yj , zj ) and zj is its elevation. (3) Topological Analysis. Once the above filtering process is finished, make connected components analysis (namely, connected component labeling). Then, the neighboring points with the same filter state are clustered as one same region. If the regions with filter state as ground meet two conditions, contained by the regions with filter state as non-ground and elevation differences between them and their adjacent non-ground regions are larger than h1 , then the regions’ filter state is relabeled as non-ground. Thus, the inner parts of large buildings can be filtered. The original filter result and our filter result are shown in Fig. 2(c) and (d), respectively, and both the filters have same filter parameters. Parameter values are set as h1 = 0.1, r = 20, h2 = 30. As the figures show, the large buildings are filtered for our proposed method while only parts of large buildings are filtered for Vosselman’s filter. h
2.2. Normalization of non-terrain measurements Once the ground measurements are acquired, various types of interpolation methods [16] such as inverse distance weighting, Kriging and triangulation with linear interpolation can be employed to derive the DTM from the ground points. In this paper, method of triangulation with linear interpolation is utilized. For any non-terrain points, the height of the point on nDSM is equal to the height of the point on DSM subtracted from the height of the point on DTM, as shown in Fig. 3. Thus, the above operations can derive the relative heights of trees to the grounds. 3. Segmentation and detection of tree points Generally, the non-terrain measurements consist of tree points, building points, wall points, bridge points, vehicle points, power line points etc. Among these non-terrain measurements, the tree points have two distinguished characteristics: the first one is the points belonging to the same tree or the points belonging to same population are scarcely coplanar; the second one is most of the tree points belong to the multiple echoes. Therefore, two operations, including automatic segmentation of LiDAR data into coplanar point clusters using a surface growing algorithm and multiple echo analysis, are employed to separate the trees points from other types of non-terrain points. 3.1. Segmentation by the surface growing algorithm Segmentation is the process of labeling each measurement in a point cloud, so that the points belonging to the same surface or region are given the same label. A number of segmentation algorithms for point cloud are available, and a literature survey of different segmentation methods has been presented in [17]. Among these approaches, we adopted the surface growing algorithm proposed by Vosselman et al. [18]. We make a short explanation of the surface growing algorithm here because of its strong relevance to our feature recognition method. The algorithm consists of the following steps. (1) Determine a seed surface. A seed surface consists of a group of nearby points that fit well to a flat plane. The algorithm selects an arbitrary unclassified point and tests if a minimum number of nearby points can be fitted to a plane. If this is the case, these points constitute the seed surface. Otherwise, another arbitrary point is tested. The key of this step is to test if a minimum number of nearby points can be well fitted to a plane, and 3D Hough transformation [18] is employed to finish this task. The principle of 3D Hough transformation is illustrated as follows. A plane in a 3D space can be described by Vosselman and Klein [19]: d = X cos α cos β + Y sin α cos β + Z sin β.
(2)
Every point p(x, y, z ) in a point cloud defines a surface in the 3D parameter space spanned by the plane parameters α , β and d. Points in the point cloud data that are coplanar will correspond to surfaces in the parameter space that intersect at a point. The α , β and d coordinates of this intersection point yield the parameter of the plane in the point cloud. Thus, the above 3D Hough transformation can be used to check if a minimum number of nearby points can be well fitted to a plane, as shown in Fig. 4(a) and (b). (2) Grow the seed surfaces. The growing can be based on one or more the following criteria.
• Proximity of points. Only points within certain distance to a seed surface, or the neighboring points (in kd-tree [20]) of a seed surface, can be added to this seed surface.
• Globally planar. For this criterion, a plane equation is determined by fitting a plane through all surface points in this seed surface. Points can only be added if the perpendicular distance to the plane is below some threshold.
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) Top view of data colored by elevation.
(c) Result of slope-based filter.
527
(b) Perspective view of the data in (a).
(d) Result of our improved slope-based filter.
(e) The non-terrain measurements derived by our method colored by their elevations. Fig. 2. Comparison of the results of the slope-based filter and our improved filter, respectively. Note the ground measurements are colored as green, and the non-ground measurements are colored as red in (c) and (d). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 5(a) shows the original laser scanned point cloud, and Fig. 5(b) shows the segmentation results of this point cloud by the surface growing algorithm. In the ideal situation, each segment should represent one of the semantic features such as roofs, walls, road segments etc. Fig. 5 also suggests that the sizes of segments can reflect the types of objects. Particularly, the sizes of segments belonging to trees, power lines and vehicles are usually smaller than those belonging to buildings. Thus, the sizes of segments can be selected as a feature to classify the point clouds.
528
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
Fig. 3. Derivation of the height of a non-terrain point from DSM and DTM.
(a) Three discrete coplanar points.
(b) 3D Hough space where the surfaces, corresponding to the above there points, intersect at a point.
Fig. 4. Voting in the 3D Hough space for plane detection.
3.2. Multiple returns analysis Katzenbeisser [21] has mentioned that one laser pulse can produce multiple reflections if the light traveling to ground illuminates several objects at different distances. Multiple reflections can occur when:
• • • •
objects are at different elevations; the laser beam has a chance to partly illuminate other objects; objects are large enough and have a sufficient reflectivity to generate echoes; the laser beam is large enough (in diameter) to illuminate several objects.
Fig. 6 shows the difference between tree points and building points. On the building areas, most of the returns are reflected as only one detectable echo, i.e., first echo and last echo are identical; also in Fig. 6(a) the single echo is reflected from low and flat vegetation. In Fig. 6(b), most of the echoes on the top of canopy are from first pulse, and then penetrated through the canopies. There are a number of surface shapes which will always provide multiple echoes, i.e., leafless deciduous trees, edges of buildings or of high trees. We now deal with the challenge of identifying tree regions by using multiple echo information associated with each segment. Laser beams used in our data collection will have a circular footprint with a radius of about 15–30 cm when they hit the surface. If more than one object of high reflectance lies within this footprint, multiple return height information will be received by the sensor. It is therefore expected that tree segments would have a higher number of multiple returns associated with them than the more planar building segments. This multiple return information is inherently point based. However, now that we have region information, we can calculate the percent of data points for each area that recorded a multiple return. By eliminating all regions where the ratio of points with multiple returns to single returns is smaller, we
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) Perspective view of the segmented point cloud in Fig. 3(e), and the points are colored by segment number.
(b) Top view of the segmented point cloud in (a) colored by segment number.
(c) The subset of singular returns colored by segment number.
(d) The subset of non-last returns colored by segment number.
529
(e) Detected tree points colored by elevation. Fig. 5. Segmentation, multiple echo analysis and classification of a point cloud. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
identify a large number of tree regions. In this work, we chose to retain segments with more than 20% multiple returns. However, this number is again dependent on the sampling technology.
530
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) Returns information of a building.
(b) Returns information of trees. Fig. 6. Differences of returns between buildings and trees.
4. Single tree detection and description Airborne LiDAR datasets with point densities of 1 point/m2 or denser offer the possibility of detecting, segmenting and modeling individual tree crowns. Single tree detection and 3D reconstruction allow counting of trees in a certain area and determination of parameters such as individual tree height, crown diameter, crow area and 3D crown shape. Researchers have made great progress. For example, Koch et al. [2] presents an approach to delineate single trees automatically in LiDAR data in deciduous and mixed temperate forests. First, in rasterized laser data, possible tree tops are detected with a local maximum filter. Afterward, the crowns are delineated with a combination of a pouring algorithm, knowledge-based assumptions on the shape of trees, and a final detection of the crown-edges by searching vectors starting from the trees’ tops. However, the above classical method could not be used in human settlements due to the various sizes and heights of the trees in human settlements. In this paper, we present a new approach to detect single trees in human settlements. Our method includes clustering, detecting local maximum with elevation, and identifying the points belonging to one tree. 4.1. Clustering In human settlements, the spatial distributions and configurations of trees vary very much. For example, the crowns of several trees may touch each other and form a continuous canopy, while some trees may be isolated from the others. Thus, we should first make a clustering analysis. In detail, construct the Delaunay triangulated irregular network (D-TIN) of tree points only using the horizontal coordinates, and then delete the long edges in the D-TIN, as shown in Fig. 7(a). Subsequently, make a connected component analysis using the above derived D-TIN data. As a result, the adjacent tree points are labeled as one cluster, and the points belonging to one isolated tree are also labeled as one cluster, as shown in Fig. 7(b). Once the clusters are determined, only the points belonging to the first pulse are remained for subsequent operations, because only the first pulse may come from the surface of a tree and the other types of pulse are sure to not come from the surface of a tree. However, we should note that the other types of pulse are quite meaningful to the clustering analysis. 4.2. Finding tree tops Suppose that the tree tops are the points whose heights are local maximums and the sizes of crowns of trees are only relevant to mean heights of trees. Another two facts are that the size of a tree crown is strongly correlated with the height of the tree belonging to a specific species, and the trees in one cluster generally belong to only one species. Thus, we first
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) Delaunay TIN of the tree points.
(c) Three tree tops are detected from data in (b).
531
(b) Two clusters are determined by connected component analysis, and only first pulses are retained.
(d) Extraction of edges of trees.
Fig. 7. Clustering, detecting tree tops, extracting edges of trees.
calculate the mean height of each cluster, and then derive a searching radius, Rtop , according to the mean height for each cluster. Subsequently, we will build a 2D kd-tree [20] for the points of each cluster only using the horizontal coordinates, and the kd-tree is capable of fast and exact searching the nearest neighbors to a fixed radius bound. Furthermore, if horizontal distance of any two adjacent tree tops is smaller than a preset value, then one of the tree tops must be deleted, which avoids detection of tree tops that belong to one same tree. As a result, we will find the points whose heights are larger than any of their nearest neighbors to a fixed radius. 4.3. Extracting tree edges Theoretically, each tree top has a local sphere of influence, which is anisotropic and directionally circular. That is, along most directions, the branches of the sphere of influence are approximately circular. Hu et al. [22] proposed the spoke wheel operator to delineate the boundary of the sphere of influence. Hu’s spoke and the spoke wheel SW are used for detecting a footprint of a pixel p, as shown in Fig. 8(a). A spoke is a line segment with a length of m (a preset value) pixels, and m = 90 pixels in Fig. 8(a). A spoke wheel is a sequence of spokes S (ϕi , m) (i = 0, . . . , 4n − 1) with common initial point p and evenly spaced angles ϕi = π i/2n, and n = 18 in Fig. 8(a). The set of pixels in a spoke wheel SW centered at the pixel p with 4n spokes is denoted by SW (p, n, m). The intersection between a spoke and a tree edge provides useful information to determine the local sphere of influence around a pixel. However, we do not know the exact location of a tree edge in advance. To search for the intersection of a spoke and the edge of a road, we start from p, move in the direction of the spoke, and observe the absolute elevation differences between p and the pixels along the spoke. The differences are small when the pixels are near to p; however, they may become larger when the pixels are far away from p. Let Si be the spoke at pixel
532
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) The configuration of a spoke wheel.
(b) Extracted tree edges by a spoke wheel algorithm operated on a tree image.
Fig. 8. Spoke wheel algorithm and its application on extraction of tree edges.
(a) Perspective view of pseudo-grid.
(b) Top view of pseudo-grid. Fig. 9. Pseudo-grid generation.
p. The cutting point, denoted by Ci , on Si is the first pixel such that
|I (Cj ) − I (p)| ≥ σi (j = 0, . . . , m, 0 ≤ i < 4n)
(3)
where σi is the intensity standard deviation on SW (p, n, m). Therefore, the threshold in Eq. (3) is adaptive. Subsequently, we connect the cutting points on all spokes around a pixel p in a counterclockwise direction, which results in a closed polygon. This represents the footprint of the pixel p, denoted by F (p). For each detected tree top, a temporary pseudo-grid is generated from the point cloud around the tree top as follows. (1) Suppose that a tree top and its preset searching radius Rtop are already known, clip out the points whose horizontal distance between them and the tree top is less than Rtop . (2) The size of pseudo-grid is calculated with average point density of laser point data. Once pseudo-grid is created, the raw laser point data is assigned to each pseudo-grid shown in Fig. 9. If there is more than one point located in the same grid, only the highest one is retained. (3) If there is no point in a grid, the central point of the grid is added into the grid, and the elevation of the central point is assigned to one of the nearest points in the point cloud. Once the pseudo-grid is generated, the above spoke wheel algorithm is performed on the pseudo-grid image, and tree edges are extracted, as shown in Figs. 7(d) and 8(b). 5. Experiments and performance evaluation In this section, two datasets are employed to extract individual tree crowns, and the performance of our proposed algorithm is evaluated.
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) Point cloud of Enschede colored by elevation.
533
(b) The 396 extracted trees by the proposed method.
(c) The extracted edges of trees are in the same color as the trees in (b). Fig. 10. Extraction of individual tree crowns from data in Enschede. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
5.1. Testing data Two point clouds are used to test our proposed approach. The first data is over the town of Enschede, the Netherlands with 700 m width and 1100 m length, as shown in Fig. 10(a). The dataset was recorded with the FLI-MAP 400 system with 20 points/m2 , and maximum number of returns per pulse reaches to 4. In this first site, there are 423 trees by manual interpretation. The second data is over the suburban of Nanning, Guangxi Province, China with 460 m width and 1000 m length, as shown in Fig. 11(a). The dataset was acquired with the RIEGL LMS-Q560 system with 10 points/m2 , and maximum number of returns per pulse reaches to 4. In this second site, there are 96 trees by manual interpretation. 5.2. Results and performance evaluation The extracted trees and their edges of the first dataset are shown in Fig. 10(b) and (c), respectively. The number of extracted trees is 397 in this dataset, and visual inspection suggests that there are 30 mistakes in this test. Among the mistakes, most of them mix more than one tree into a tree, as shown in Fig. 12(a), and some of them are the powerlines. As far as the second dataset is concerned, the extracted trees and their edges are shown in Fig. 11(b) and (c), respectively. There are 87 trees being detected, and visual inspection suggests that there are 5 false detections in this test. All the false detections mistake a cluster with more than one trees as one isolated tree. The two tests also tell us that if the number of points of one tree is not large enough, our algorithm will miss them. The two tests suggest that the performance of our approach is excellent. As far as the accuracy is concerned, 92.44% and 94.25% are reached for the first test and the second test, respectively. As far as the completeness is concerned, 86.76% and 85.41% are acquired for the above two tests, respectively. The statistics tell that our approach has a promising application in practice. However, our method also has some disadvantages. First, the proposed method will omit the small trees if they only have several points in the dataset. Second, our method will mistake adjacent trees as one tree if the searching radius is not well preset. As last, our method will mistake some other objects as trees if the objects also produce multiple echoes. 6. Conclusions and discussion Extraction of individual tree crowns from airborne LiDAR point clouds is very meaningful for many applications. However, most of the existing methods are oriented to forested areas, and they are not fitful to extraction of the trees from point clouds
534
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
(a) Point cloud of Nanning colored by elevation.
(b) The 87 extracted trees by the proposed method.
(c) The extracted edges of trees are in the same color as the trees in (b). Fig. 11. Extraction of individual tree crowns from data in Naning, China. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(a) Side view of points belonging to two adjacent trees but mistaken as one.
(b) Some power line points are mistaken as trees.
Fig. 12. False detection in extraction of trees.
in human settlements. In these areas, the scenes are more complex than the forested areas. Thus, a method is proposed to extract individual tree crowns from airborne LiDAR point clouds in human settlements. In our method, an improved slopebased filter is employed to separate the non-ground measurements from the ground measurements, the surface growing algorithm is utilized to segment the point clouds into segments, the multiple echo analysis is used to distinguish the tree points from other types of non-ground measurements, and the spoke wheel algorithm is employed to get the edges of each tree. Two datasets are used to test our proposed method. Experiments show that our approach is capable of extracting more than 85% trees from the point clouds with the accuracy no less than 95%, which suggest its promising applications in engineering. However, the proposed method fails to extract the small trees, and will mistake mixed trees as one isolated tree. Future work will focus on the improvement of the proposed method. Acknowledgments This research was funded by the Basic Research Fund of Chinese Academy of Surveying and Mapping under Grant Nos 7771025, 7771109 and 7771017, the Project for Young Scientist Fund sponsored by the Natural Science Foundations of
J. Liu et al. / Mathematical and Computer Modelling 58 (2013) 524–535
535
China under Grant Nos 41001280 and 40901195, and the China Postdoctoral Science Foundation under Grant 2010047038, respectively. References [1] G.G. Parker, Structure and microclimate of forest canopies, in: M. Lowman, N. Nadkarni (Eds.), Forest Canopies, Academic Press, California, USA, 1995, pp. 73–106. [2] B. Koch, U. Heyder, H. Weinacker, Detection of individual tree crowns in airborne LiDAR data, Photogrammetric Engineering & Remote Sensing 72 (4) (2006) 357–363. [3] L. Wang, A multi-scale approach for delineating individual tree-crowns with very-high resolution imagery, Photogrammetric Engineering & Remote Sensing 76 (4) (2010) 371–378. [4] D.V. der Zande, W. Hoet, I. Jonckheere, J. van Aardt, P. Coppin, Influence of measurement set-up of ground-based LiDAR for derivation of tree structure, Agricultural and Forest Meteorology 141 (2006) 147–160. [5] F.M. Danson, D. Hetherington, F. Morsdorf, B. Koetz, B. Allgöwer, Forest canopy gap fraction from terrestrial laser scanning, IEEE Geoscience and Remote Sensing Letters 4 (1) (2007) 157–160. [6] J. Hyyppä, M. Schardt, H. Haggrén, B. Koch, U. Lohr, H.U. Scherrer, R. Paananen, H. Luukkonen, M. Ziegler, H. Hyyppä, U. Pyysalo, H. Friedländer, J. Uuttera, S. Wagner, M. Inkinen, A. Wimmer, A. Kukko, E. Ahokas, M. Karjalainen, HIGH-SCAN: the first European-wide attempt to derive singletree information from laserscanner data, The Photogrammetric Journal of Finland 17 (2001) 58–68. [7] H. Duong, R. Lindenbergh, N. Pfeifer, G. Vosselman, ICESat full waveform altimetry compared to airborne laser scanning over the Netherlands, IEEE Transactions on Geoscience and Remote Sensing 47 (10) (2009) 3365–3378. [8] T.L. Erdody, L.M. Moskal, Fusion of LiDAR and imagery for estimating forest canopy fuels, Remote Sensing of Environment 114 (4) (2010) 725–737. [9] J. Secord, A. Zakhor, Tree detection in urban regions using aerial LiDAR and image data, IEEE Geoscience and Remote Sensing Letters 4 (2) (2007) 196–200. [10] G. Vosselman, Slope based filtering of laser altimetry data, International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences 33 (Part B3/2) (2000) 935–942. [11] K. Zhang, S. Chen, D. Whitman, M. Shyu, J. Yan, C. Zhang, A progressive morphological filter for removing non-ground measurements from airborne LIDAR data, IEEE Transactions on Geoscience and Remote Sensing 41 (4) (2003) 872–882. [12] P. Axelsson, DEM generation from laser scanner data using adaptive tin models, International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences 33 (Part B4) (2000) 110–117. [13] K. Kraus, N. Pfeifer, Determination of terrain models in wooded areas with airborne laser scanner data, ISPRS Journal of Photogrammetry and Remote Sensing 53 (4) (1998) 193–203. [14] G. Sithole, Segmentation and classification of airborne Laser scanner data, Delft, ISBN 9061322928, p. 184, Publications on Geodesy of the Netherlands Commission of Geodesy, vol. 59, May 24, 2005. [15] G. Sithole, Filtering of laser altimetry data using a slope adaptive filter, In: IAPRS 2001, vol. XXXIV, 3/W4, Annapolis (MD), USA, 2001, pp. 203–210. [16] M.J. De Smith, M.F. Goodchild, P.A. Longley, Geospatial Analysis: A Comprehensive Guide to Principles, Techniques and Software Tools, second ed., Troubador Publishing Ltd., 2007. [17] M. Wang, Y.H. Tseng, Automatic segmentation of LiDAR data into coplanner point clusters using an octree-based split-and-merge algorithm, Photogrammetric Engineering and Remote Sensing 76 (4) (2010) 407–420. [18] G. Vosselman, B.G.H. Gorte, G. Sithole, T. Rabbani, Recognizing structure in laser scanning point clouds, International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences 46 (8/W2) (2004) 33–38. [19] G. Vosselman, R. Klein, Visualisation and structuring of point clouds, in: G. Vosselman, H.-G. Maas (Eds.), Airborne and Terrestrial Laser Scanning, Whittles Publising, ISBN: 978-1904445-87-6, 2010, pp. 43–79. [20] S. Arya, D.M. Mount, N.S. Netanyahu, R. Silverman, A.Y. Wu, An optimal algorithm for approximate nearest neighbor searching in fixed dimensions, Journal of the ACM 45 (1998) 891–923. [21] R. Katzenbeisser, Echo detection, Technica Note TopoSys. Toposys homepage, www.toposys.com, 2003 (accessed 07.06.11). [22] J. Hu, A. Razdan, J.C. Femiani, Road network extraction and intersection detection from aerial images by tracking road footprints, IEEE Transaction on Geoscience and Remote Sensing 45 (12) (2007) 4144–4157.