APPLICATION OF A PATTERN RECOGNITION ALGORITHM FOR SINGLE TREE DETECTION FROM LiDAR DATA

In the present study, we applied the Particle Swarming Optimization (PSO) procedure to parametrize two Local Maxima (LM) algorithms and a pattern recognition model based on raster and point-cloud datasets in order to extract treetops of coniferous forests from high resolution LiDAR-data of different forest structures (monoplane, biplane and multi-layer) in the Alps region. The approach based on the pattern recognition model uses the geomorphon algorithm applied to the DSM to detect the treetops. The geomorphon model gave good results in terms of matching rates (Rmat: 0.8) with intermediate values of commission and omission rates (Rcom: 0.22, Rom: 0.2). Therefore, it could be a valid alternative to the LM-algorithms when only raster products (DSM – CHM) are available. The geomorphon pattern recognition model has been proved to be a powerful method in order to properly detect treetops of coniferous stands with complex forest structures. This model allows to obtain high detection rates and estimation accuracy of forest volume, also in comparison to the most recent available literature data. The models are developed in Java under Free and Open Source license and are integrated in the JGrassTools library, which is now available as SpatialToolbox of the GIS gvSIG.


INTRODUCTION
Starting from nineties the use of Light Detection and Ranging (LiDAR) technology for forest inventories and management has gained more and more importance because of its suitability to capture high-resolution topographic data on large areas (basin or regional level) due to its capability to penetrate canopy and therefore obtain information about the vegetation, the objects on the surface and the ground (Dubayah and Drake, 2000;Lim et al., 2003a).The use of remote sensing technologies is particularly interesting in mountainous regions where the progressive abandonment of the areas increased the problems of accessibility for mapping and therefore for harvesting.A proper characterization of the 3D-structure of a forest involves the extraction from the LiDAR data of the most important structural and biometrics forest parameters, such as tree height and diameter distribution, canopy cover and forest biomass.There are two possible approaches to extract these structural attributes from a discrete return LiDAR point cloud: the individual tree-and the plot-based methods (Wulder et al., 2008).The main difference between these two approaches is that the former focuses on the metrics of each single tree of a given forest, while the second aims at developing a number of site-specific empirical statistical relations to estimate forest attributes without any interest for the biometrics of the single trees forming the forest.Considering the available literature in this regard (Popescu and Wynne, 2004;Naesset and Nelson, 2007), it is clear that to make the most of the LiDAR data the single-tree approach should be preferable.However, for both approaches the validation procedure consists of comparing field and LiDARderived measurements at single-tree or plot scale.Many algorithms were suggested for the single-tree extraction method based both on raster and on raw point cloud data with different performances as function of the forest structure and point density.Some of these methods were tested and compared in other research studies under different forest conditions (Kaartinen and Hyyppä, 2008;Kaartinen et al., 2012;Vauhkonen et al., 2011), focusing on the performance of the models and on the effects on the detection rates of the given study areas and forest types.Applications in Alpine environment are also documented (Eysn et al., 2015).In this study, we present the results of three different models for tree top extraction integrated in the automatic procedure based on the Particle Swarming Optimization (PSO) method to extract single trees on different forest structures and types in a test area in Val Aurina (Bolzano -Italy).The models are two local maxima and a pattern recognition approach for single-tree extraction algorithms applied on both raster and point cloud LiDAR data.In comparison with previous international studies (Kaartinen and Hyyppä, 2008;Kaartinen et al., 2012;Vauhkonen et al., 2011, Eysn et al., 2015) we used fewer algorithms, but we integrated them all in a PSO procedure which automatize both the calibration and matching processes.The use of PSO helps in finding the best combination of parameters for each model and site and at the same time assures the use of the same boundary conditions for all the simulations.In order to focus on those trees that are of real interest to a forest inventory, we decided to work only on trees with diameter greater than 17.5 cm.The aim is to validate the pattern recognition approach for single-tree extraction focusing on the performance of the three models for raster and pointcloud datasets on the different forest types.

Study area
The study area is part of the Aurina Valley located in northeastern Italy in the Autonomous Province of Bolzano.The total area covered by the LiDAR survey is approximately 10 km 2 with the altitude ranging from 1059 to 2406 m a.s.l. and an average elevation of 1729 m a.s.l.Coniferous mixed forests are the dominant forest type of this area with Norway spruce (Picea abies) (60%), larch (Larix decidua) and stone pine (Pinus cembra) as dominant species.Most of the area is of public property and managed by the local community of Aurina Valley, with a limited number of private properties.The site was selected in 2012 for its high variability in forest structure, which represents a very frequent situation in Alpine environment and a complex case study.For this area, high precision LiDAR as well as detailed field surveys data were available.

LiDAR survey
The LiDAR data were collected in one single flight in September 2012 at an average flying height of 500 m a.g.l..An Optech GEMINI Airborne Laser Terrain Mapper (ALTM) was used by Helica s.r.l. to collect the data.Maximum 4 returns were recorded from each pulse.The final LiDAR point cloud had an average density of 10 points/m 2 .LiDAR-data were collected, processed and delivered by Helica s.r.l (Udine -Italy) in Universal Transverse Mercator (UTM) coordinate system, in WGS84 -ETRF2000 datum (Franceschi et al., under revision).The DTM and DSM at a resolution (pixel size) of 0.5 m were elaborated by the surveying company using TerraSolid TerraScan with the TIN-based filtering method.

Field sample survey
For the study areas, detailed field data are available.Specifically, position, diameter, species and height of each tree with diameter at breast height (DBH) bigger than 5 cm were recorded.The diameter at breast height (DBH) was measured using a diameter tape and the height of the trees using a vertex hypsometer.Twelve circular plots were chosen to represent three different vertical forest structures (i.e.monoplane, biplane and multiplane), that were preliminary selected based on the analysis of the available data of the forest structures supplied by the local forestry department and then reassigned based on the distribution of the heights of the measured trees in each plot.Figure 1 shows a schematic representation of the distribution of the heights and DBH for the trees with diameter greater than 17.5 cm in all the plots of the selected study area.

Positioning:
The position of the center of each plot together with the position of the trees within each plot were estimated using the application for digital field mapping Geopaparazzi, an Android app for engineering and geology surveys.Through this tool it is possible to store georeferenced notes (text, images and personalized forms) and log GPS tracks (http://geopaparazzi.github.io/geopaparazzi,Brovelli et al., 2016, Leidig andTeeuw, 2015).In Geopaparazzi we considered both the GPS signal and the Canopy Height Model (CHM) loaded as background layer.Since the CHM image was used as reference, the tree position was referred to the top of the trees.Using this application, the relative planimetric accuracy of positions of each tree is estimated to be less than 1 m.

Tree volume:
In order to calculate the stem volume from DBH and height of the single tree, we applied a sitespecific allometric relationship developed by Floreancig (2014) using the data collected in all the 12 plots of the study area.

Estimating single tree position and height from laser data
A single-tree top extraction method based on the identification of Local Maxima (LM) (Hyyppa et al., 2001;Persson et al., 2002) applied to both raster and point cloud data and a pattern recognition raster based model (Geomorphon), were the approaches used in the present study.LM are points with higher elevation respect to the others in a window of a given size and shape.Geomorphon, instead, characterizes the surface in the neighboring of a central cell through the use of the line of sight principle (Yokoyama et al., 2002, Jasiewics andStepinski, 2012).These techniques are based on the assumption that the trees have a tip and this is much more reliable with conifers rather than other species.In order to optimize the parametrization of the used algorithms (e.g. the size of the window or the thresholds on the zenith and nadir angle) a Particle Swarming Optimization method was used to compare the results of each simulation with the data measured in the field using a complete automated matching procedure (Franceschi et al., under revision).Specifically, the data analysis and optimization method consisted of five steps: 1. run the algorithm on a set of randomly defined parameters, the result is a set of tree positions and heights from the input dataset; 2. compare positions and heights of the LiDARextracted tree tops with those of trees mapped in the field using an automatic matching procedure, which considers both the distance and the difference in height between the LiDAR-derived and the fieldmeasured trees; 3. asses the obtained results based on a target function that takes into account the number of true positives (TP: LiDAR-extracted trees that match for position and height the field-measures ones), false positives (FP: LiDAR-extracted trees that do not exist in the field) and false negatives (FN: trees measured in the field that cannot be paired to a LiDAR-extracted tree tops); 4. after the first simulation, based on simple mathematical functions the parameters of the algorithms are iteratively changed to move toward the best possible solution of the target function; 5. the optimizer iterates until either the target solution has been reached or the number of maximum imposed iterations has been got.

2.4.1
Local maxima extractor for raster: The LM algorithm, when applied to raster data, scanned the raster canopy height model (CHM) by using a moving window and identifying local maxima as cells in the center of the moving window with the highest height value in comparison with all the cells covered by the window.All the details regarding this model and the way it has been integrated in the PSO process are described in Franceschi et al. (under revision).

2.4.2
Local maxima extractor for point-cloud data: The LM algorithm compared the height of each point against all the points included in a circular window centered on the point itself.If the central point resulted to be the highest point, then it was automatically considered as a tree top candidate and added to a temporary list of extracted trees.All the details regarding this model and the way it has been integrated in the PSO process are described in Franceschi et al. (under revision).

2.4.3
The pattern recognition extractor for raster: The third method that we implemented in the PS workflow is based on the geomorphon algorithm (Jasiewics and Stepinski, 2012) applied to the raster of CHM.Geomorphons are a finite number of fundamental micro-structures of landscape which can be extracted from a reference surface or Digital Elevation Model (DEM).The geomorphon concepts base on those of pattern recognition for image classification and has been originally developed by Jasiewics and Stepinski (2011) for the morphological classification of landforms.In this application for single-tree extraction the reference surface is the CHM.The most important difference between this method and the LM methods consists in the evaluation of the extension of the area considered to identify the shape of the canopy.The LM algorithm considers the elevation of a fixed-size and fixedshape neighborhood window regardless of the specific local conditions, while in the geomorphon model the size and shape of the moving window self adapt to the local topography using the principles of the line-of-sight (Lee, 1991;Nagy, 1994;Yokoyama et al., 2002) and ensuring that each shape is identified as its appropriate scale.The advantage of using the line-of-sight instead of the grid based concept for analyzing the neighborhood of a central cell is that using the line-of-sight it is possible to identify the shapes of a surface on a wider range of scales in comparison with the grid based neighborhood.The original model proposed by Yokoyama et al. (2002) considers the surface relief and horizontal distance related to vertical angles, the so-called zenith and nadir angles, along the eight principal neighboring directions.The zenith angle, in particular, is the angle between the horizontal plane and the line connecting the central cell with the point located on the profile at the line-of-sight distance or at the maximum allowed distance.This angle is negative if the central cell has an elevation higher than the one of the connected point on the profile, positive otherwise.The nadir angle is the angle between the nadir and the hypothetical line-of-sight obtained by reflecting the surface profile in the current direction to the horizontal plane.The developed model calculates the map of the 10 most common geomorphon structures from the input CHM and implements the single tree extraction algorithm by isolating the combination of these structures that better reproduce the presence of a tree.We first compared the output geomorphon raster map with the field observations in each plot and thus, isolated the most common combinations, which usually identify a treetop and its surrounding area.After some tests we found that the best combination for all vegetation type and trees species consists in peaks which are surrounded by: i) no more than 2 not valid cells (otherwise it is considered a point on the boundary of a crown), ii) no other peaks, iii) not a coexistence of pits, valleys, spurs and hollows.This means that, usually, treetops extracted from the geomorphon peaks do not have in their surrounding other sharp morphologies, on the contrary, there should be slopes, ridges or shoulders which indicate a continuity of the shape that generates the central peak.A simplified symbolic 3D representation of the main geomorphon structures is shown in Figure 3.
As for the LM model, very often the branches create the morphological conditions for the identification of a local maximum, therefore also in this case we introduced a filter based on the difference in height between the height of the candidate treetop and that of all the neighboring cells contained in the moving window.The filter requires to set this vertical threshold difference between neighboring cells representative of the difference between the final apex of a branch in respect to the lower vegetation or ground layers.This value was added as parameter of the model to be calibrated.Therefore, by the PSO procedure the following three parameters were optimized: the maximum lookup distance to consider for extracting the line-of-sight, the flatness threshold as the difference between the nadir and the zenith angles and the maximum difference in height between the cells of the moving window for identifying the branches.The result of the point-cloud-based treetop extractor was a set of tree positions and heights and a raster map containing the reclassification of the input map (CHM) in the 10 main morphological structures.

2.4.4
Particle Swarming Optimization: The Particle Swarm Optimization (PSO) (Eberhart and Kenedy, 1995) is an iterative computational method used to find an optimal solution given a measured set of data.The fitness of this optimization process is based on the results of a fitness function, whose mathematical structure should be chosen according to the aim of the analysis (e.g.forest volume assessment).In our application the fitness function is based on the comparison between the lidar-extracted and the field-mapped trees (Franceschi et al., under revision).This comparison is based on an automatic tree matching algorithm that identifies the number of true positives (TP), false positives (FP) and false negatives (FN).This tree matching procedure is based on the assumptions that the horizontal distance between the field-mapped tree and the lidarextracted treetop cannot be more than 2 m and the difference in height between field-mapped and lidar-extracted tree should be less than a percentage of the measured height of the tree (Monnet et al., 2010).At the end of the matching procedure a list of extracted trees matching field-mapped trees (TP), a list of extracted trees that have no matching tree and are assumed to be non-existing (FP) and a list of the remaining unmatched mapped trees (FN) is produced (Franceschi et al., under revision).

2.4.5
Fitness function: PSO searches the best solution using a fitness function which calculates the index of goodness based on the comparison between LiDAR extracted and existing trees using the statistical concepts of Information Retrieval, precision and recall (Makhoul et al., 1999).In this study, for all the algorithms, we implemented three different fitness functions for: minimizing the number of the trees extracted by the algorithms but not mapped in the field (FP), minimizing the existing trees not found by the software (FN) and the classic combination between precision and recall with harmonic mean (F-score).The optimization of the fitness function is performed by maximizing its value in the range between 0 and 1. Starting from the harmonic mean of precision and recall, traditionally called F-measure or balanced F-score, where FP and FN are evenly weighted, we used the general Fβ measure (with nonnegative real values of β) to give β-times more important to precision or recall and consequently to FN or FP.The two values of β used in this study are β = 2 for minimizing the FN and β = 0.01 for minimizing the FP (Franceschi et al., under revision).The model converges when the fitness function is maximized (approximately 1), which means that all extracted trees have a match in the mapped dataset, while no false positives or false negative have been generated.Even if the ideal situation of reaching 1 does not occur, and the optimization stops when it reaches the maximum number of iterations, it is assumed that the best available solution for the current set of parameters has been reached.This is only partially true, since the particle swarming method itself doesn't assure to reach an optimal solution, but it has proved to be acceptable during the testing.

Data analysis
For each combination of LM algorithm and fitness function, the following statistical parameters were calculated for each plot: • Total number of LIDAR-extracted trees (Next); The height of each single LIDAR-extracted tree was calculated as the elevation of the LiDAR data in the point where a LM was detected and the volume of the aboveground biomass (AGB) was estimated for both field-mapped and the LiDAR-extracted trees through a field-derived species-independent allometric function relating the tree volume to the tree height and DBH (diameter at breast height) as follows: (4) where a = constant DBH = tree diameter at 1.3 m of height H = total tree height.The variable a is derived from the non linear regression of the field observations.Since the volume of both field-measured and LiDAR-extracted trees was estimated using the same equation we did not introduce an additional source of discrepancy between the two estimates.However, since only position and height of trees were extracted from the LiDAR data, the DBHs were estimated using the following hypsometric relationship derived from the site specific field measurements: (5) where b = 0.0096 c = 1.298 d = 0.00

Detection rates
According to the three different optimization functions, we compared the detection rates of the LM-algorithms based on raster, point-cloud and geomorphon for each vegetation structure.
The matching and extracted results per method and per fitness function (Rmat, Rext) indicate how well each combination performed for each vegetation structure, while the omission (Rom) and commission (Rcom) rates indicate the errors in detecting the trees.
In particular, the highest extraction rate (Rext: 1.24) was obtained by the raster based LM method using the fitness function which minimizes the FN on multilayers structures, while the lowest rate (Rext: 0.63) was registered by the point-cloud-based method by minimizing the FP, in any case on multilayers forests structures.The geomorphon algorithm shows intermediate values of extraction rate with a little overestimation of the total number of trees (about 10%) in particular with biplane and multilayer stands by minimizing the FN.
In general, the two raster-based methods showed the same behaviour with highest average values of extraction rates for biplane forests structures using the harmonic mean or the min_FN as target functions.In contrast, for the matching rates, the method with the highest matching rate was the point-cloud based (Rmat: 0.88) with minimization of FN on biplane forests structures, while the one with the lowest (Rmat: 0.59) was the LM raster-based method using min_FP on monoplane structures.The geomorphon methods again had performance in the middle of the other two and in general very closed to the point-cloud LM method in terms of matching rate.Considering all the structures and the fitness functions, the overall Rmat of the three methods were very similar, but the geomorphon algorithm and raster-based LM algorithms detected more FPs, leading to a higher Rcom with a general overestimation of the total number of trees.Regarding the impact of the vegetation structure, the average matching rates ranged between 0.71 for multilayer and 0.80 for biplane forests, while the commission rates ranged between 0.15 for monoplane and 0.21 for biplane structures and the omission rates between 0.18 for biplane and 0.29 for multilayer coniferous forests.
Considering the incorrect detection, the highest commission rates (Rcom: 0.43) was produced for the raster-based method by minimizing the FN, while the lower commission rates (Rcom: 0.00) were obtained with the PS point-cloud based method by minimizing the FP on monoplane and multilayer forest structures.In general, the method with the best omission rates was the point-cloud-based method when min_FN was applied.On average, the Rmat was higher when min_FN was applied, but Rcom was significantly higher using min_FN in comparison to min_FP function.Therefore, for the test area, the fitness function was the parameter that had the most relevant impact on the correct detection of trees.
On the whole, irrespective of the target function, the pointcloud-based LM-algorithm seemed to be the less sensitive to the type of function that was used for the calibration of the parameters with high values of the matching rate and the lowest values of the commission rates.The geomorphon algotithm showed better results rather than the raster-based LM one and candidated itself as a good choice to be used in case of only raster data are available for the study area.

Tree heights and forest volume
The optimization process was based on an automatic matching procedure, which used the distance and the difference in height between the LiDAR-extracted and the field-measured trees as parameters.Therefore, the heights of the LiDAR-extracted trees necessarily fitted very well those of the field-measured ones.
In the study area the error in the evaluation of the height of the trees was significantly affected only by the forest structure, in particular, it was lower for the monoplane and higher for the multilayers plots, respectively.The method and the fitness function did not prove to be able to affect the tree height estimation accuracy.
In general for the study area the algorithm based on LM pointcloud performed better than the raster-based ones.The rate volume of LiDAR extracted TP trees (RVtp) respect to the total field-measured volume ranged between 0.60 and 0.96 and averaged 0.87 for the LM point-cloud, while it averaged 0.80 and 0.84 for the raster based LM and geomorphon algorithms respectively.However, no statistical difference between the extraction methods, the optimization functions and the forest structure where found in terms of extracted volume.Although, on average the results in terms of total LiDARextracted volume were similar between all the models, the geomorphon algorithm showed to be less sensitive than the LMalgorithms to the selected target function.

DISCUSSION AND CONCLUSION
Methods and algorithms to extract treetops from LiDAR data are usually calibrated on the subjective choice of a number of parameters, which are defined considering the characteristics and the spatial distribution of trees of the studied forests.However, this approach has often produced variable and inconsistent results according to the structure, species composition and development stage of the forest.In the present study, we applied the particle swarming optimization procedure to parametrize two LM algorithms and a pattern recognition model in order to extract treetops of coniferous stands with different forests structures from high resolution LiDAR-data.
The obtained results were compared between each other.The results in this study showed remarkable differences between the algorithms more in terms of tree detection rates than in height and volume estimation accuracies.Similar results were obtained by Vauhkonen et al., (2011) in which they compared six different algorithms on coniferous and broadleaved forests in Europe and a tropical plantation in Brasil.There were two main differences between the algorithms developed and tested in this study.First, one algorithm used point cloud data whereas the others require a raster model of the surface (CHM or DSM-DTM).Second, two algorithms used the local maxima approach and the other a pattern recognition approach to identify the treetops (as maxima).
The new developed method based on the geomorphon algorithm seems to be very promising as it showed detection rates of treetops (Rmat: 0.77 -0.86) very similar to those obtained with the LM method based on point-cloud (Rmat: 0.76 -0.88) and higher than those obtained with LM method based on raster data (Rmat: 0.72 -0.82).Furthermore, even if the estimation of the volume of the above ground biomass, were similar between all the models, the geomorphon algorithm showed to be less sensitive than the LM-algorithms to the selected target function.
On overall the use of the PSO iterative parametrization process to calibrate the parameters of the models seems to have big influence on the quality of the results in terms of treetop detection and volume estimation.This confirms that a rough parametrization of the selected algorithm can produce a strong reduction of the estimation accuracy even when the forest structure is well known.
A high accuracy can be reached when a high matching rate is combined with a low commission rate.In this context, the pointcloud-based method combined with the harmonic mean as target function was found to be the best compromise, as it allowed us to get high matching rates associated with the lowest commission and omission rates (Rmat: 0.81, Rcom: 0.03).On the contrary, the raster-based LM algorithm (Rcom: 0.13 -0.27) and the geomorphon model (Rmat: 0.86, Rcom: 0.24) tended to overperform due to high commission rate.The geomorphon model in general for all the forest structures showed higher values of matching rates respect to the other two LM algorithms (Rmat: 0.80) with similar values of omission and commission rates to those of the LM raster based (Rcom: 0.22, Rom: 0.20) leading to a general estimation of the number of the trees very closed to the real value with an average overestimation of only 1% on all the forest structures.
The new developed model based on geomorphon (Jasiewics and Stepinski, 2012) gave good results in terms of matching with intermediate values of commission and omission rates.Therefore, it could be a valid alternative to the LM-algorithms when only raster products (DSM -CHM) are available.
In conclusion, this study validated and compared the use of a pattern recognition algorithm to extract treetops form high resolution LiDAR data.The geomorphon model has proved to be able to properly detect treetops of coniferous stands with complex forest structures.When integrated in an automatic optimization procedure (PSO) it allowed us to obtain high detection rates and estimation accuracy also of the forest volume, in comparison to the most recent available literature data.The differences in performance between the methods were found to be higher for tree detection than for height and volume estimation.
Further studies are needed to test and validate the method on a larger datase and in particular on broadleaves forests.Finally, the effect of different point density and flight parameters on the detection accuracy of the proposed method should be performed to understand its robustness and sensitivity to the most common variables of LiDAR dataset.

Figure 1 :
Figure 1: Distribution of the heights and DBH for the trees with DBH greater than 17.5 cm in all the plots of the study areas.The colors represent the forest structures: green = monoplane structure, orange = biplane structure, blue = multilayer structure.
Figure 2 illustrates the concepts of line-ofsight distance and zenith and nadir angles.Therefore, evaluation of the Local Ternary Pattern (LTP) in each direction depends of three main parameters, the lookup distance and the two vertical angles (Jasiewics and Stepinski, 2012).Finally, there are 498 possible output patterns resulting from the application of the basic model for Local Binary Pattern recognition, therefore a reclassification is needed to convert these combinations to the 10 most common geomorphons types: flat, peak, ridge, shoulder, spur, slope, pit, valley, footslope and hollow (Jasiewics and Stepinski, 2012).

Figure 2 :
Figure 2: Main concepts of the geomorpon algorithm: zenith and nadir angles.

Figure 3 :
Figure 3: Simplified symbolic 3D representation of the main geomorphon structures used to identify the treetops.