INDIVIDUAL TREE OF URBAN FOREST EXTRACTION FROM VERY HIGH DENSITY LIDAR DATA

Airborne LiDAR (Light Detection and Ranging) data have a high potential to provide 3D information from trees. Most proposed methods to extract individual trees detect points of tree top or bottom firstly and then using them as starting points in a segmentation algorithm. Hence, in these methods, the number and the locations of detected peak points heavily effect on the process of detecting individual trees. In this study, a new method is presented to extract individual tree segments using LiDAR points with 10cm point density. In this method, a two-step strategy is performed for the extraction of individual tree LiDAR points: finding deterministic segments of individual trees points and allocation of other LiDAR points based on these segments. This research is performed on two study areas in Zeebrugge, Bruges, Belgium (51.33° N, 3.20° E). The accuracy assessment of this method showed that it could correctly classified 74.51% of trees with 21.57% and 3.92% underand over-segmentation errors respectively. * Corresponding author


INTRODUCTION
Tree extraction from remotely sensed data, in addition to applications in mapping and three-dimensional modelling, is worthwhile from the aspect of environmental term and urban green space management.So far, many studies on the tree extraction have been done using aerial and satellite images (Kalapala, 2014;Immitzer et al., 2012;Larsen, 2007).With the increase in the spatial and spectral sensor resolutions, computation power and reducing in the cost of data gathering, the level of detail for automatic individual trees extraction has been increased.However, the vertical structure of trees, tangled crowns and automatic detection in the level of details of remote sensing data, are tree extraction challenges that researchers are faced with them.Nevertheless, aerial image interpretation is impeded by different spectroradiometric distortions caused by vignetting effects (caused by more light reaching the centre of an image than the edges), atmospheric variations, sun-target-sensor geometry or view/illumination geometry and topographyinduced illumination variations (Lillesand et al. 2008).Airborne LiDAR can penetrate in tree crown, provide the geometric information and show some internal structure of the trees, and also these capabilities may increase by increasing in point density.With the high ability of airborne LiDAR to provide three-dimensional information directly and to receive multiple return signals from vegetation, it has become a viable alternative to image.The single tree detection using airborne LiDAR data was primarily proposed by Brandtberg (1999), Hyyppä and Inkinen (1999) and Samberg and Hyyppä, (1999), which subsequently became an accepted topic for investigators (Hyyppä et al. 2008).Most published algorithms detect individual trees from an interpolated raster surface from LiDAR points that hit on the tree canopy surface (Zhang et al., 2015).Table 1 shows some of these methods.For example, in the developed algorithm by Koch et al., a Digital Crown Height Model (DCHM) is first calculated by subtracting the height value of the Digital Terrain Model (DTM) at each pixel from the height value of the Digital Surface Model (DSM).In the smoothed DCHM, tree tops are then detected with a local maximum filter.The Pouring algorithm is then used to obtain crowns (Koch et al., 2006 The CHM can have inherent errors and uncertainties from a number of sources.For example, spatial error can be introduced during the interpolation process from the point cloud to the gridded height model (Guo et al., 2010), which can decrease the accuracy of tree segmentations and relevant measurements (Li et al., 2012).To increase the accuracy of detection, Persson et al. applied a Gaussian function to elevation model (Persson et al. 2002), but the use of smoothing filter cause to estimate the tree height incorrectly (Tiede et al., 2005).As the amount of softening increases, the possibility of losing trees top increases.This can cause the indiscrimination of small trees (undersegmentation error) and vice versa (e.g. the very low softening level can cause to create additional tops and over-segmentation error).In the study of Kumar, the impact of the degree of softening on the number of detected peaks is clearly visible on table 2. Hence, in these algorithms, the number and the locations of detected peaks heavily effect on the process of detecting individual trees (Kumar, 2012).As small trees are dominated by larger trees in dense forests, it is impossible to find them by pixel analysis (Reitberger et al., 2009).In the study of Duncanson et al. intermediate trees are over predicted (commission error) and understory trees are often undetected (omission error) (Duncanson et al., 2014).To overcome this problem, Alonzo et al., tried to control the distance between the detected peaks by applying one or more thresholds (Alonzo et al., 2014).But one or more thresholds may not enough to identify a variety of trees (variety of species, age, environmental conditions and growth rate) and so trees don't be correctly identified.Unlike the above methods, some methods use LiDAR-point clouds directly to extract individual trees points.Table 3 shows some of these methods.It's hard to define top points of deciduous or young trees under the other tree canopy in the point cloud (Wang et al., 2008).According to the three-dimensional coordinates of LiDAR points, Li et al. used a distance threshold to tree segmentation.But it's hard to define this threshold, especially in dense forests.Inappropriate threshold may cause to under-and/or over-segmentation errors (Li et al., 2012).Another problem is that points must participate in the segmentation and classification processes one by one, that in itself rises the computation time.Lu et al. offered a method to detect tree trunks points based on higher levels of intensity.In their method, the tree segmentation begins from bottom of the trees.Consequently, their method will be influenced by bushes in some areas (Lu et al., 2014).
All above mentioned methods search seed points as tree top or bottom to start segmentation algorithm.Hence, the number and the locations of the selected seed points effect on the detection process directly.In our study, a new method to extract individual trees points from LiDAR data with no need to extract seed points is developed.The main idea of this method consists in extracting deterministic segments of individual tree points and investigating membership of other points to them.

GENERAL SCHEMA
As mentioned above, a two-step strategy is performed for the extraction of individual tree LiDAR points.These steps are (1) finding the Deterministic Segments of Individual Trees points (DSITs), and (2) the allocation of other LiDAR points based on these DSITs.Figure 1 indicates a presentation of the different phases of this proposed approach.In continue of this paper, details of this method are explained.
Figure 1.The flowchart of proposed algorithm

Data Pre-processing
Pre-processing of LiDAR data in this study consisted two phases: (1) Because of using urban forest point cloud for tree segmentation, we first separated urban forest point clouds class from other urban objects class (such as roads and buildings) using a city plan.
(2) The surface of the ground is not flat and our algorithm is based on the relative spacing of each point, therefore we first need to eliminate the effect of land.For this purpose, a digital elevation model (DEM) was first generated by interpolating the ground points using ordinary kriging (Guo et al., 2010).The LiDAR point cloud was then normalized by subtracting the DEM height from the LiDAR point cloud (Lee et al., 2010).According to the general shape of trees, a cylindrical shape is considered as a DSIT (Figure 2).These cylinders are defined based on a direction and a radius which should be extracted for each tree.

Obtaining Cylinders' Direction
Tree growth of each tree is in line with its tree trunk.Hence, the central axis of the cylinder should be defined along the tree trunk.The steps of extracting tree trunk in this paper are implied as following steps: firstly, the bushes and the ground points are removed from the point cloud by applying a maximum and minimum height thresholds which are defined experimentally (Figure 3a) based on the study area.In this research, these thresholds were equal to 1.4 m and 4.5 m.Then, in order to classify the resulted point cloud to trunk and non-trunk classes, third thresholding is performed on intensity of each point (Figure 3b).This experimental threshold is obtained from this fact that the number of echoes of a pulse related to the tree trunks and its branches are usually less than the number of echoes related to the leaves.Hence, the intensity of the trunk and branches are higher than the intensity of the leaves (Lu et al. 2014).With this thresholding, the leaves points are removed from the tree trunk class as far as possible.In this research, the intensity threshold was equal to 25. OPTICs algorithm is a clustering method based on density data; this means that wherever the data density is high, there may be a pattern.The input parameters of this algorithm are minimum points and the neighbourhood radius parameters, which are indicated with K and ε respectively.A point is considered as a core point, if there are at least K points in the ε radius for it.
Connected core points represent a pattern and the representative of this pattern is a core point that point density around it, is more than the rest.
The output of the OPTICs algorithm is a curve that represents this structure in two-dimensions.In this curve vertical axis, RD, is each point's reachability distance and horizontal axis, order, is the order of points in the reachability graph.Every valley in this curve represents a cluster.We used a Matlab code of OPTICs algorithm that set the value of ε, automatically (Daszykowski et al. 2002).In this study, an appropriate value for K is considered by examining the half of trunk class at plot 1 (Figure 4a).The number 8 was obtained for K by trial and error method.The resulted plot of OPTICs algorithm has been shown in Figure 4b.By putting any point into cluster of its nearest representative, tree trunk point's clusters are obtained (Figure 5a).The obtained Trunk Clusters are accompanying with some bushes points and should be eliminated to extract final individual tree trunks.This elimination is done by a thresholding operation on their Euclidean distance from mean of each cluster.The final trunks point's clusters are obtained after this operation (Figure 5b).In the last step of obtaining cylinders' direction, by extracting the trunk of single trees, their directions are extracted using principal component analysis (PCA). Figure 5c shows extracted trunk direction of two neighbouring trees.

Obtaining Cylinders' Radii
According to the trees growth, as the diameter of a tree trunk is larger, its crown is more expanded.Crowns are tangled with each other in dense areas.Therefore, the radii of the cylinders are defined as a function of two parameters: distance and diameter.If the tree diameter is large or tree is farther from its neighbour trees, the radius of its cylinder should be defined large.As mentioned above, we express the steps of our method for two neighbouring tree.Cylinder radius of the tree 1 and 2 (r 1 and r 2 respectively), are defined by the following equations: If trunk directions of two neighbour trees are approximately parallel to each other (Figure 6a), their radiuses are determined based on equations 1: and (1) Where = distances between centres of tree trunks , = trunk diameter for tree 1 and 2 respectively.
If two neighbour trees are getting away from each other (Figure 6b), their radiuses are determined based on equations 2: and (2) Where , = trunks direction angles for target and neighbour trees, respectively = is determined based on equations 3 ( ) and (3) Where = perpendicular distance of i-th tree trunk top point from intersection line of perpendicular planes on tree trunk centre line.If two neighbour trees are approaching to each other and trees trunk centre lines are intersecting (Figure 6c), same equation with previous case is applied, but we limit the cylindrical zones with a break plane.This plane creates and angles with plane perpendicular to the central line of the target and neighbour tree respectively.Thus, this plane allows changing cylinder radius along with height changes.This radius change helps cylinders to keep the concept of DSITs.Finally, as shown in the figure 6d, if two neighbour trees are approaching to each other and trees trunk centre lines are skew, then we define the truncated cones instead of cylinders.In this case we define radii of these truncated cones as equations 4: , , and Where = nearest and the farthest neighbour tree trunk respectively, , = radius of the truncated cones for target tree, , = radius of the truncated cones for neighbour tree.
By having DSITs for each pair of neighbouring trees:

DSITs Inner Points Extraction
In continue of our proposed algorithm, for each neighbouring trees pair, other tree points out of defined DSITs are extracted.The steps of this section are developed as follows: 1. Firstly, tree trunk centre of gravity of each cluster is obtained and a Triangulated Irregular Network (TIN) model is constructed on them.This TIN model is used to determine trees topology and extract neighbour trees.In addition to trees topology, TIN model is applied to restrict the number of points which are under classification process.
2. A trunk of a single tree is selected as the target tree, randomly, and its neighbours are sorted from nearest to the farthest distance.Then, the target tree and its nearest neighbour are selected as a trees pair and their DSITs are defined.The points between the trees pair that are inside of DSIT target/neighbour tree , may be related to the target/neighbour tree and indicated as initial target/ neighbour points (neighbour points) (Figure 7).  4. Neighbour points must be removed from point cloud for the target tree and be added to the point cloud for the next target tree.Initial target points must be involved in the process of segmentation for the next neighbour tree.
5. The steps 2, 3 and 4 are repeated for the target tree until the last of its neighbour tree.
6.The latest initial target points are the final target points and must be extracted as that target tree points.
7. The steps 2 to 6 are run until the last target tree.

TESTS AND ACCURACY ASSESSMENT
Our method is performed on LiDAR data (IEEE GRSS Data Fusion Contest, 2015) which was acquired in Zeebruges, Belgium (51.33°N, 3.20° E; 20 m above sea level) (Figure 9), on March 13, 2011.The point density for the LiDAR sensor was approximately 65points/m², which is related to point spacing of approximately 10cm.The study sites are plot 1 and 2 with an area of 4900 m2 and 3300 m2 respectively.The trees of both plots are deciduous.
Figure 9.The location of study areas In order to validate of our proposed method, we compared the classification result with manual extracted trees from aerial ortho-images (Figure 10).For this purpose, the ortho-images are  4. Segmentation results compared to aerial image superimposed on the LiDAR point cloud.Then, the points of each tree are extracted by an expert operator.This comparison shows that our algorithm, correctly classify 74.51% of trees with 21.57% and 3.92% under and over segmentation errors respectively (Table 4).

CONCLUSIONS AND FUTURE WORK
In this study, we proposed a new method to extract individual trees points from LiDAR point data.The basis of this method is to define cylindrical areas for each tree and to allocate points to each tree based on these areas.This process allowed us to determine the class of large number of points and to use a process to segment less number of points in non-deterministic segments at the same time.The results of the evaluation show that the proposed method has good potential to extract individual trees points Future work could be focused on improving this algorithm in mixed forests.Tree crowns have a very complex structure, and change from one species to another and from one place to another largely.In individual tree extraction operations, there are some points in non-deterministic areas between trees that membership verification process of them is complex especially in very high density LiDAR data of semi-intensive urban forests.
Prior knowledge of trees in the study area, helps to choose a better shape for DSITs and then to calculate DSITs parameters.
In our scheme, the deterministic and non-deterministic area around each tree is separated by cylindrical zones.Other geometric volumes can be used.The tuning processes in the proposed methodology is partly dependent on the type and height of the trees in the study area.For example, 1.4 m and 4 m threshold values for maximum and minimum height thresholding to extract initial tree trunks.Obviously these thresholds are area specific and needs re-tuning for trees in different environments.Nevertheless, much research work still needs to be conducted with a view to achieving complete elimination of the tuning stage from the entire individual tree extraction process.
The used LiDAR data only had intensity and threedimensional coordinates of the points.Full-waveform analysis of LiDAR data can reach better results.

Figure 3 .
Figure 3. Initial classifications for trunk class extraction

Figure 4 .
Figure 4. Trunks point cloud of half of trees of Plot 1 and (b) output of OPTICs algorithm that was carried out on this point cloud.

Figure 5 .
Figure 5. (a) The results of clustering with OPTICs algorithm, (b) the final trunks and (c) trunk direction of two neighbouring trees

Figure 6 .
Figure 6.Visual expression of equations parameters: (a) trees are parallel, (b) trees are getting away from each other and (c and d) trees are approaching to each other.

Figure 7 .
Figure 7.The DSITs definition for couple tree

Figure 8 .
Figure 8. Allocation of points outside of the Deterministic Segment of Individual Trees (DSITs) to their trees.

Figure 10 .
Figure 10.Aerial image of the study areas (Plot 1 and Plot 2) A general view of part of the results is presented in Fig. 11.

Figure 11 .
Figure 11.Final segmentation results related to trees from half of plot 1.
. Different individual tree detection methods based on LiDAR-derived raster surface

Table 3 .
Different individual tree detection methods based on LiDAR-point cloud This contribution has been peer-reviewed.doi:10.5194/isprsarchives-XLI-B3-337-2016 Figure 2. The schematic concept of a DSIT for two trees