DETECTION OF SINGLE TREE STEMS IN FORESTED AREAS FROM HIGH DENSITY ALS POINT CLOUDS USING 3 D SHAPE DESCRIPTORS

Airborne Laser Scanning (ALS) is a widespread method for forest mapping and management purposes. While common ALS techniques provide valuable information about the forest canopy and intermediate layers, the point density near the ground may be poor due to dense overstory conditions. The current study highlights a new method for detecting stems of single trees in 3D point clouds obtained from high density ALS with a density of 300 points/m. Compared to standard ALS data, due to lower flight height (150-200 m) this elevated point density leads to more laser reflections from tree stems. In this work, we propose a three-tiered method which works on the point, segment and object levels. First, for each point we calculate the likelihood that it belongs to a tree stem, derived from the radiometric and geometric features of its neighboring points. In the next step, we construct short stem segments based on high-probability stem points, and classify the segments by considering the distribution of points around them as well as their spatial orientation, which encodes the prior knowledge that trees are mainly vertically aligned due to gravity. Finally, we apply hierarchical clustering on the positively classified segments to obtain point sets corresponding to single stems, and perform `1-based orthogonal distance regression to robustly fit lines through each stem point set. The `1-based method is less sensitive to outliers compared to the least square approaches. From the fitted lines, the planimetric tree positions can then be derived. Experiments were performed on two plots from the Hochficht forest in Oberösterreich region located in Austria. We marked a total of 196 reference stems in the point clouds of both plots by visual interpretation. The evaluation of the automatically detected stems showed a classification precision of 0.86 and 0.85, respectively for Plot 1 and 2, with recall values of 0.7 and 0.67.


INTRODUCTION
Accurate measurements of forest structure are increasingly required across large areas to support a wide range of activities related to sustainable forest management.Cost-effective and more automated methods are needed to provide tree attribute data for forest ecosystem services.Therefore, remote sensing tools are increasingly being used to survey forest structures.However, the spatial extent and spatial resolution of a given sensor are inversely related.Airborne Laser Scanning (ALS) has become a key tool for gathering information on 3D structure of forests (Wulder et al., 2012).The derived information from ALS data can grant detailed estimations of forest characteristics and single tree analysis (Yao et al., 2012;Maltamo et al., 2012).Previous studies are showed several properties of single trees such as species, height and crown properties which can be measured with high resolution ALS data (Maltamo et al., 2012).However, limited persistence has been done on the stem detection of single trees.For instance, single tree stems have been determined from the interpolated CHM (Canopy Height Model) at the highest positions (Solberg et al., 2006) or by using hierarchical clustering for stem reflections and reconstructions with a RANSAC-based adjustment (Reitberger et al., 2007).Due to the low point density and lack of information about the reflection characteristics, minor focus has been given to tree reconstruction using laser hits on the stems (Reitberger et al., 2007;Polewski et al., 2016).In order to detect stem of trees, data with greater precision would be required to allow a more accurate representation of the actual discontinuities in the single trees (Vauhkonen, 2010).On the other hand, Terrestrial Laser Scanning (TLS) has also been proven to be a suitable method for obtaining very detailed information about geometry of trees in forests (Liang and Hyyppä, 2013).Pfeifer and Winterhalder (2004) showed a method for reconstructing the cross section of tree stems and branches from TLS data with free-form curves.However, the study indicated that expected cross section reconstruction works only satisfyingly, if the branch is covered with points from all sides which it is not possible due to occlusions.Liang et al. (2012) presented a fully automatic algorithm based on single-scan TLS data for stem detection and mapping with the overall accuracy of 73%.The stem points are established using classification based on the local covariance matrix features.The provided method is capable of giving good parameters only when the points are evenly distributed on the tree trunk; however the applicability of feature estimation for a group of points has not been taken into account.Lindberg et al. (2012) projected candidate stem points onto a 2D plane.They applied Hough transforms to locate circles, representing the potential stems.Heinzel and Huber (2016) used a 3D voxel grid transformation of the input TLS point clouds to detect tree stems using morphological operations and empirical rules.They reported detection rates of 84% to 97% for the number and location of stems depending on the tree DBH (diameter at breast height).Olofsson et al. (2014) by voxelizing the point cloud and analyzing the influence of different height layers could estimate tree stem positions with an average proportion of 87%.Wang et al. (2016) performed in the first step, RANSAC based circle fitting of projected stem points, which is later followed by RANSAC cylinder fitting in 3D space.Polewski et al. (2017) proposed a statistical framework for detecting cylinders based on accumulation and voting in parameter space.The size of the accumulation cell is determined automatically using bandwidth selection methods for kernel density estimators, which relaxes the requirement of setting this critical parameter manually or through trial-and-error.The method is applied on a dense 3D point cloud for mapping fallen tree stems.Based on the mentioned studies, the main advantage of TLS data lies in its capacity to scan a sample plot in forest horizontally and vertically in detail.However, important factors such as the occlusion effect and relatively high cost of the instrument transportation from site to site are negatively effecting the use of TLS data.Moreover, the coregistration of several scans covering a study area is an essential step in the interpretation of multi-scan data acquisition.The fully automated registration between several scans at the point level is still challenging.In standard operational applications of ALS, the flight heights are usually between 400 to 800 m, resulting in point densities up to 30 points/m 2 .However, an alternative scenario is also possible, where the flying altitude is significantly decreased to below 150 m, for example using a helicopter mounted system or even a UAV with around 50 m flight height.This can be seen as a middle ground between standard ALS and TLS techniques, trading off large area coverage for increased point density.This tradeoff is due to the lower flight altitudes associated with this technique compared to standard ALS campaigns, which results in smaller footprints.Razak et al. (2011) used high resolution DEMs (Digital Elevation Models) extracted from high density ALS data to semi-automatical recognition of morphometric landslide features even under forest canopy.Höfle et al. (2012) provided an example of high density ALS data potential to use for urban vegetation detection purposes.They used a high point density of 50 points/m 2 .Khosravipour et al. (2014) also presented an algorithm which is able to create a pit-free CHM raster using full waveform ALS data with 160 points/m 2 density.The algorithm significantly improves the accuracy of tree detection compared to local maxima based methods.The data collected by high density ALS systems is less precise in comparison with TLS.However, within an equal time frame, the area that can be investigated by utilizing high resolution ALS is significantly larger than the area investigated with TLS.Also the aforementioned TLS based methods for stem detection are not practical for the applications using ALS data, since the curvature shape of the tree trunk in ALS point clouds is not captured as detailed as in case of TLS.On the other hand, the increased point density, resulting from lower flight height can provide more details in the point clouds compared to standard ALS, enabling the use of 3D spatial descriptors to locate individual tree stems.Therefore, an automated method for stem mapping within high-density ALS data is interesting from forestry application point of view.The objectives of this study are to develop a new method for single tree stem detection based on high density ALS data using (i) point and object part level 3D shape descriptors, and (ii) 1 regularized line fitting with a prior on orientation.A further objective is (iii) to assess the accuracy of detected tree stems.This paper is motivated by the successful application of high density ALS systems for precise monitoring of vegetation and forest structure, reported in the aforementioned studies.Also, detected tree stems could be used to improve the 3D segmentation algorithm as prior knowledge in terms of the detection rate and the position of trees.
The remainder of this work is structured as follows: Section 2 focuses on the details of our approach.Section 3 shows the experiment results.Finally, the results are discussed with conclusions in Sections 4 and 5.

METHOD
We adapt the method of Polewski et al. (2015) which is originally designed for fallen tree segmentation, to detect the standing stems of single trees from unstructured high density ALS point clouds.The main goal is to detect linear structures in the ALS 3D point clouds which are likely to represent single tree stems.The output of our method is a set of 3D lines which correspond to detected stems.The pipeline describing our approach is presented in Fig. 1.Our approach is a three-tiered detection procedure at (i) point, (ii) segment and (iii) object levels.The segment term refers to the grouping of points within a fixed length cylindrical neighborhood which are likely to represent part of a tree stem.Objects refer to entire tree stems which are composed from groups of similarly aligned segments.The method proceeds as follows.First the likelihood of points belonging to a tree stem is estimated.Second, the segments containing the highest probability stem points are detected in the 3D point clouds.Finally the segments are merged through hierarchical clustering to produce single tree stems.In the following, we explain the steps of our method in detail.
Figure 1.Overview of stem detection pipeline.

Point level
In the input data depending on the forest characteristics, different objects are present such as ground vegetation, regenerations, standing tree stems etc.The focus of this step is to obtain for every point, an estimate of the probability that it belongs to a tree stem.High density ALS data can provide a range of features related mainly to the 3D structure of single trees.The features derived from 3D point clouds can be grouped into three categories: 1. Point feature histograms (PFH): a local 3D shape descriptor of the neighborhood around the target point, based on the angular relationships between adjacent surface normals.It is useful for distinguishing between different types of surface classes based on their shape (plane, cylindrical, spherical, etc) (Rusu et al., 2008;Polewski et al., 2015).
2. Covariance features: set of features derived from eigenvalues of local neighborhood covariance matrix around target point (see Weinmann et al. (2015)).3. Normalized height: the relative point heights over the Digital Terrain Model (DTM).
For the point classification level we chose Random Forest (Breiman, 2001) as a binary classifier due to its robust ability to estimate the class probability.

Segment level
This level focuses initially on generating segment candidates.For each point with high probability, a cylindrical neighborhood with constant radius rseg and height lseg is defined.Afterwards, all the points inside the cylinder space are taken into account to perform the least-squares Orthogonal Distance Regression (ODR) (Al-Subaihi and Watson, 2004).This is done by eigenanalysis of the point coordinates' covariance matrix: the ODR line's direction is the eigenvector corresponding to the maximum eigenvalue, and this line passes through the points' centroid.The ODR line becomes the segment's axis.We classify the candidate segments generated in the previous step into the 'positive' and 'negative' groups.The 'positive' group represents the segments which are really parts of tree stems, and the 'negative' contains branches, ground vegetation, and etc.The first set of segment features are derived from a modified version of Cylindrical 3d Shape Context (CSC) built around the segments' axes (Polewski et al., 2015) (see Fig. 2).A spatial distribution histogram of points within the cylindrical volume around the segment axis can then be constructed.The point counts inside the histogram bins form the features for the classification.The second set of features are calculated based on the angular deviation of the segment axes from the world Z.The segments with the deviation larger than θ thr are rejected without regard to remaining feature values.Additionally, the point probability statistics are extracted as another group of features for classification.In the current study, we looked at the points inside the segment neighborhood and we created the stem point probability histograms quantized at 0.2 intervals, resulting in 5 features.

Object level
The goal of current level is to take the stem segments with high probability from the previous step and merge them to reconstruct individual tree stems.This is based on the collinearity and mutual distances between segments.

Hierarchical clustering
In the next step of the pipeline, the representative 'positive' segments are merged together.For this purpose, first we have a combination of two distances to calculate.The aggregate distance d between segments Si and Sj is the weighted sum of angular deviation dA and the spatial distance between point centroids dC (see Eq. 1).
In the Eq. 1, S refers to the axis and S indicates the point centroid of each segment.w1 is the weight component for the spatial distance dC .Fig. 3 shows the angular deviation dA and the spatial distance dC between two segments surrounded by cylinders.In the current clustering process the distance D between two clusters Ci and Cj is defined as the largest distance d from any combination of member segments.The stopping criterion consists of two conditions.First, when considering two clusters Ci, Cj if the spatial distance dC between any pair of segments S k ∈ Ci, S l ∈ Cj is bigger than a predefined threshold dC,max then the merging of Ci and Cj is aborted.On the other hand, the ODR line is fitted to the set of points belonging to segments in Ci, Cj and the orthogonal projected distance of all points to the ODR line is calculated.If the projected distance of any point is greater than a set threshold dp,max or the angular deviation between the fitted line and the Z axis exceeds dA,max, this also terminates the merging for that cluster pair.Fig. 4 represents the stopping criterion in the hierarchical clustering phase between segments.
Figure 4.The stopping criterion for clustering the segments.The green lines represent the segments in cluster Ci and the magenta lines show the segments in cluster Cj, respectively.The red line is the line fitted to both clusters' points.The red points outside the cylinder have projected distance greater than dp,max.

Stem line fitting
In the final step, for each cluster extracted from the merging step, containing all the stem points, we execute the line fitting procedure.For the fitting, we apply orthogonal distance regression with the 1 norm of residuals as the error criterion and a prior on the line's verticality.The problem is formulated in terms of an energy minimization.The total energy for a line L, Etot(L) is presented in the Eq. 2. The total energy Etot has 2 components, the data goodness-of-fit term E d (Eq. 3) and the angular prior term Ea (Eq.4).The α element refers to the balance coefficient between energy terms.The angular prior was considered due to the prior knowledge that the tree stems grow almost always vertically.Furthermore, the 1 norm is a more robust estimator and less sensitive to outliers available in the data compared to the 2 norm (Al-Subaihi and Watson, 2004).In our experiment, the 2 norm in the segment level and the 1 norm is used in the object level, respectively.This decision is related to the computational expenses of the 1 and 2 fittings.Usually, the number of the segments to process is several orders of magnitude higher compared to the stems obtained from hierarchical clustering.On the other hand, the 1 version is much more computationally expensive than the 2 based method, which makes it intractable to apply for all segments.
In the Eq. 3, element Pi, i ∈ 1..n refers to the ith point inside the cluster.Furthermore, Z element in the Eq. 4 represents Z axis of world coordinate system.L(P ) term is the projection of point P onto line L and w indicates L's direction vector.The line L is parametrized using 4 values a = [a1, ..., a4] (Al-Subaihi and Watson, 2004) showed in the Eq.5: where t is the location parameter.Note that the z component of w is always positive with a value of 1, and therefore the angle between the world Z axis and w always lies in the interval [0; π 2 ].This allows us to drop the absolute value on the angular prior term Ea since the arccosine is guaranteed to be positive in the mentioned interval.For the optimization of the orthogonal distance fitting, a two-step method similar to Liu and Wang (2008) as well as Watson (2002) is used.The first step for the re-parametrization computes the projection of all fitted points Pj onto the current line L(a, t) to minimize the distance from Pi to L: The line positions tj corresponding to the orthogonal projection onto L can be recovered using Eq. 7, where p0 = [a1, a3, 0]: In the second step, we minimize the energy similar to Etot, but with the projected line positions fixed at {tj}, with respect to shape parameters only, i.e.: For minimizing the objective (Eq.8) the BFGS quasi-Newton method is applied (Wright and Nocedal, 1999).The derivative of the energy's data term with respect to any parameter a k is expressed as follows (rj = Pj − L(a, tj) is the j − th residual): By splitting the computation into two steps and fixing the tj values, the problem is considerably simplified, because otherwise the terms tj in Eq. 9 would depend on a k , leading to the necessity of calculating the derivative of the product ∂[tj(a) • w(a)]/∂a.Therefore, the derivatives with respect to the entire parameter vector are: As for the angular deviation term Ea, only the axis parameters a1, a3 have a non-zero derivative, given by: We iteratively execute steps 1 and 2 in sequence until convergence is reached.Due to the non-convexity of the optimization problem, local optima may exist.To remedy that, we per-form multiple restarts of the optimization with randomly initialized starting line hypotheses, and pick the lowest-energy solution.This yields the final fitted lines representing the detected single tree stems for all clusters.

Materials
Experiments were conducted in the Hochficht forest close to Holzschlag in Oberösterreich region which is located in Austria.The study area is a type of mountainous forest with a high structural complexity, dominated by Norway spruce (Picea abies), European beech (Fagus sylvatica) and Silver fir (Abies alba).Two sample plots with the approximate area size of 6844.7 m 2 and 17907 m 2 and a mean tree density of 272 trees/ha were selected in the mixed forest to construct the experiment.The high density ALS data was acquired in leaf-on condition by the FMM GmbH Company with the VUX-1 scanner integrated in the VP-1 pod in September 2015 with an average point density of 300 points/m 2 .We assumed that the data is given in a georeferenced coordinate system.The flying altitude between 150-200 m resulted in an average footprint size of 88 mm.Fig. 5 shows a sample scene in the 3D point clouds associated with the visible single tree stems.

Training classifier
We used parts of test plots and two additional plots to train the classifier in the point and segment levels.In the available 3D point clouds, a significant percentage of the stems is not at all represented (specially for the deciduous trees).For every single tree with visible stem, the points and segments were manually marked and assigned either the stem or the non-stem class by visual interpretation, resulting in training sets for binary classification.The total number of the marked points and segments were 20000 and 564, respectively.This represented about 28% of the total number of the stem points and about 3% of the generated segments in plots 1 and 2.

Reference data
The schematic class labels which groups individual segments into tree stems were later obtained based on the relative segment positions and orientations using visualization.In some cases a single tree stem was represented in the point cloud, but it was missing from the reference data, due to the lack of evidence in the 3D point clouds.The total number of the labeled stems were 196.

Choice of parameters
The various control parameter values that we used in our experiment for each level is summarized in the Tabel1.1.Control parameters for the single tree stem detection method; np refers to the number of points inside the clusters.

Evaluation
In the current experiment we use the "recall" and "precision" measures to characterize the detection performance between detected and reference tree stems.The "recall" is defined as the ratio of the reference stem numbers which have at least one associated detected stems to the total number of reference stems.The "precision" expresses the count of detected stems that were successfully connected to reference stems as a fraction of the total number of detected tree stems.We considered the detected and reference stems as matched if the average projected distance between them was not more than 30 cm.This value was derived based on the maximum DBH of trees in the target area.

RESULTS AND DISCUSSION
The procedure for stem detection was applied to the both plots.The output of the single tree stem detection consists of a number of point sets which correspond to the individual stems.The method takes the advantage of the increased point density, which makes more laser reflections available underneath the canopy for regions of test plots dominated by conifers, due to the smaller footprint size compared to the standard ALS.In contrast, the deciduous tree stems are missing in the point clouds due to the dense canopy cover in leaf-on state, and no benefit was achieved despite the lower footprint size.In case that only sparse understory is below the tree base height, stem points are successfully detected by the expressed classifier training and stem line fitting method.Fig. 6 shows the stem detection results for a sample plot (mixed with deciduous and coniferous trees) in three main levels of point, segment and object.The sample plot contains deciduous and coniferous trees.In Fig. 6a the classification results at the point level based on the high and low point probability on the tree stems are presented.The "positive" and "negative" groups of segments are classified using shape context, angular deviation from the world Z and point probability statistics in the Fig. 6b.Finally, at the object level stems with the fitted ODR line (orthogonal distance regression with the 1 norm) after merging are indicated in the Fig. 6c.In the current figure, we used the minimum stem point probability threshold of 0.6 to remove low probability points from the analysis.The detection performance of the proposed method is presented by Fig. 7.Note that in the current test plots, due to the point density and forest characteristics (particularly deciduous trees) up to 30% of the visible tree stems could not be detected.Specifically, in the lower canopy layer limited number of tree stems can be found since most of them are covered by taller tree stems and understory vegetation.Therefore, the majority of the detected stems are located in the upper and intermediate layer of the forest.
Here, we focus on the detection accuracy of tree stems that are derived from the 3D point clouds.In the plot 1, the good trade off results between recall and false positive rate are 0.84 and 0.25, respectively.Also, for the plot 2, the results are 0.63 for the recall and 0.13 for the false positive rate is achieved.By increasing the precision rate, the recall rate is decreased.The average rate of false detected tree stems in the plot 1 and 2 amount to 0.25 and 0.21, respectively.However, no improvement is achieved in the lower layer and deciduous trees since (i) laser hits at the stem of small trees rarely happen, (ii) the stem points are missing for the trees with compact canopy.The restrictions of the approach are that only trees in visible stem in the 3D point cloud visualization can be detected.This problem could be alleviated by acquiring an even denser point cloud e.g. by using UAV-based laser scanning, where more stem hits are to be expected.As mentioned before, the method fails in the regions with high concentration of deciduous trees where stem hits are rare and stems points cannot be clearly clustered.In case of deciduous trees, it is not clear if an increase of the nominal point density of a ALS data with reducing flight height can provide more visible stems.Perhaps conducting the data acquisition in leaf-off state could remedy this problem.The current point density, did not allow to reconstruct the cylindrical shape of stems.Therefore the stem diameter estimation requires higher point density which is left for the future work.From the application point of view, the processing of large areas requires implementing a tiling scheme due to the memory usage of the spatial index structures as well as the large number of generated segment candidates.Moreover, the transferability of 3D descriptors is not perfect, in the sense that in our experiment we had to train and test our method on the same area.Another limitation is related to the optimization problem of the line fitting step.Due to the nonconvexity, which does not guarantee global optimality, it might happen that the fitted line converage to local optimum.

CONCLUSIONS
The study presents a novel method for detecting stems of single trees based on the high density ALS Our results demonstrate that the classification precision is achieved to 0.86 and 0.85, respectively for two sample plot 1 and 2, with recall values of 0.7 and 0.67.In future work, we would like to utilize the stems obtained from the proposed method to enhance the segmentation and delineation of single trees by providing prior knowledge about tree locations.Additionally, higher point density could lead to more precise reconstruction of the stems.

Figure 3 .
Figure 3.The aggregate distance d between two segments; the angular deviation dA is shown with green arrow and the spatial distance between point centroids dC with red arrow.

Figure 5 .
Figure 5.A point cloud visualization of sample forest scene with multiple visible stems (point clouds of scene colored by height over DTM).

Figure 6 .
Figure 6.The detection results for a sample plot:(a), (b) and (c) correspond to point, segment and object levels, respectively (see Sections 2.1-2.3).At point level in (a), the red color shows low and blue high probability.The solid green bars in (b) indicate tree stems classified as positive and red bars refer to unmatched tree stems with references.The points which do not belong to the detected stems are removed from analysis and colored as cyan.Orange ellipses outline examples of the false alarms.The fitted magenta lines in (c) represent the reference tree stems which overlap with colored detected stems (ODR with the 1 norm).

Figure 7 .
Figure 7. Single tree stem detection performance with probability threshold of 0.4 to 0.6 for the test plots.
The values were assigned experimentally based on the forest characteristics.