Individual Tree Crown Segmentation in Two-Layered Dense Mixed Forests from UAV LiDAR Data

In forests with dense mixed canopies, laser scanning is often the only effective technique to acquire forest inventory attributes, rather than structure-from-motion optical methods. This study investigates the potential of laser scanner data collected with a low-cost unmanned aerial vehicle laser scanner (UAV-LS), for individual tree crown (ITC) delineation to derive forest biometric parameters, over two-layered dense mixed forest stands in central Italy. A raster-based local maxima region growing algorithm (itcLiDAR) and a point cloud-based algorithm (li2012) were applied to isolate individual tree crowns, compute height and crown area, estimate the diameter at breast height (DBH) and the above ground biomass (AGB) of individual trees. To maximize the level of detection rate, the ITC algorithm parameters were tuned varying 1350 setting combinations and matching the segmented trees with field measured trees. For each setting, the delineation accuracy was assessed by computing the detection rate, the omission and commission errors over three forest plots. Segmentation using itcLiDAR showed detection rates between 40% and 57%, while ITC delineation was successful at segmenting trees with DBH larger than 10 cm (detection rate ~78%), while failed to detect trees with smaller DBH (detection rate ~37%). The performance of li2012 was quite lower with the higher detection rate equal to 27%. Errors and goodness-of-fit between field-surveyed and flight-derived biometric parameters (AGB and tree height) were species-dependent, with higher error and lower r2 for shorter species that constitute the lowermost layer of the forest. Overall, while the application of UAV-LS to delineate tree crowns and estimate biometric parameters is satisfactory, its accuracy is affected by the presence of a multilayered and multispecies canopy that will require specific approaches and algorithms to better deal with the added complexity.


Introduction
Laser scanning is a well-established and consolidated technology used extensively for environmental monitoring and mapping. In forest monitoring and inventories, airborne laser scanning (ALS) started to be used since 1990s [1] with focus on tree height measurements (e.g., [2]), individual tree identification (e.g., [3][4][5]), canopy structure assessment (e.g., [6]), forest volume estimation (e.g., [7][8][9]) and stand basal area estimation (e.g., [10,11]). The applications of LiDAR (Light Detection and Ranging) technology have rapidly developed and today ALS systems are routinely used over extensive areas to collect stand level and regional wall-to-wall forest mensuration data for forest management [12,13]. In Drones 2020, 4, 10 3 of 20 technology that facilitates mobile surveys in GPS-denied below-canopy forest environments mounting a UTM-30LX LiDAR sensor (Hokuyo Automatic Co., LTD, Japan). Ref. [34] used UAV-based LiDAR data acquired by a VUX-1UAV laser scanner (Riegl Laser Measurement Systems GmbH, Austria) over a forest area to extract detailed forest and ground information finding that UAV-based laser scanning is providing both high-quality forest structural and terrain elevation information. Ref. [35] estimated the DBH from a UAV point cloud acquired through a VUX-SYS LiDAR sensor (Riegl Laser Measurement Systems GmbH, Austria) based on modeling the relevant part of the tree stem with a cylinder. Ref. [36] proposed a new concept of autonomous forest field investigation to collect data above and inside the forest canopy by integrating an autonomous driving UAV with a Puck LITE laser scanner (Velodyne LiDAR, Inc., USA). Ref. [37] used an integrated UAV-LS system to collect and process LiDAR data for biodiversity studies (i.e., 3D habitat mapping) in three forest ecosystems using a Puck VLP-16 laser scanner (Velodyne LiDAR, Inc., USA). Ref. [38] tested a LiDAR-hyperspectral image fusion method in treated and control forests with varying tree density and canopy cover as well as in an ecotone environment, integrating an HDL-32E LiDAR sensor (Velodyne LiDAR, Inc., USA). Ref. [39] estimated tree volume through quantitative structure modeling applied on UAV-LS data acquired with a RIEGL Ri-COPTER with VUX-1UAV (Riegl Laser Measurement Systems GmbH, Horn, Austria). Ref. [40] explored a carbon model applied to CHM produced using a cloud point collected by UAV-LS over a deciduous forest using a Velodyne VLP-16 (Velodyne LiDAR, Inc., USA) LiDAR sensor. Ref. [41] estimated the DBH of urban trees by a linear regression model based on variables extracted from a CHM obtained from the UAV-based laser scanning system developed by [32].
Most recent studies were aimed to assess whether UAV-LS point clouds are qualified for individual tree mapping and modeling in various forest stands [42], how their accuracies are by comparison with terrestrial laser scanning [43], whether DBH and tree position estimation algorithms developed for ALS are suitable for UAV-LS data [44] and finally, how the level of processing automation impacts the tree mapping and modeling accuracy [15]. While the accuracy of UAV point clouds is overall high [32], the evaluation of further data processing algorithms for higher level products delivery in forestry is necessary for a comprehensive assessment of UAV-LS especially in complex forest structures. For example, Ref. [44] evaluated five existing segmentation algorithms (i.e., raster-based region growing, local maxima centroidal Voronoi tessellation, point-cloud level region growing, marker controlled watershed and continuously adaptive mean shift) to determine the most suitable method for individual tree detection across a species-diverse forest.
The aim of this study is to assess the capability of the low-cost UAV-LS system developed in [14] to obtain structural and dendrometric parameters in a complex layered forest. Specifically, ITC were delineated by means of the benchmarked [45] tree detection methods developed by [46] and by [47] on ALS data. The DBH and aboveground biomass (AGB) were estimated from crown area and tree height of segmented trees by allometric equations specific for Nearctic-Temperate mixed forests-Gymnosperm [48]. The main research questions of this study are: (i) is the performance of ITC segmentation based on UAV-LS comparable to ALS? (ii) which is the impact of ITC algorithm parameters setting on the number of detected trees, DBH, height and aboveground biomass estimates? (iii) is the accuracy of estimations satisfactory for inventory purposes especially in layered complex forest?

LasUAV System
The system used for LiDAR data acquisition, called LasUAV according to [14] integrates a multi echo laser scanner (LUX 4L, ibeo Automotive Systems GmbH, Germany) able to record up to 3 returns per pulse from targets at a maximum range of 200 m. It scans at a frequency of 12.5 Hz, at a horizontal field of view between −30 •  which uses up to 72 tracking channels from GPS/QZSS L1, GLONASS G1, BeiDou B1, Galileo E1, space-based augmentation system (SBAS). A second GNSS receiver (EVK-M8QCAM, u-blox Ltd., Switzerland) was used to synchronize the laser, INS and GNSS data streams. Data were stored in a mini PC with 64 GB capacity (Zbox Pi321 pico mini PC, Zotac International Ltd., Taiwan). The system weight was 1.5 kg. During operation, a local WLAN generated by a mobile router was used to connect the LasUAV computer, the RTK GNSS receiver and a laptop ground station, to remotely operate the sensors, perform diagnostic checks and control the acquisition. A Matlab (MathWorks, Inc., Natick, MA, USA) custom software was developed to generate the LiDAR point cloud in LAS (LASer) file format in post-processing.
LasUAV was mounted in a 12 brushless motors hexa-coaxial UAV rotor, equipped with a control and navigation system capable of autonomous flight via uploaded waypoints. The net endurance was around 15 min with this payload, contained within a case which was fixed in a rigid isolated frame ( Figure 1). A more detailed and thorough illustration of LasUAV components and point cloud generation performance is reported in [14].
Drones 2020, 4, x FOR PEER REVIEW 4 of 20 B1, Galileo E1, space-based augmentation system (SBAS). A second GNSS receiver (EVK-M8QCAM, u-blox Ltd., Switzerland) was used to synchronize the laser, INS and GNSS data streams. Data were stored in a mini PC with 64 GB capacity (Zbox Pi321 pico mini PC, Zotac International Ltd., Taiwan). The system weight was 1.5 kg. During operation, a local WLAN generated by a mobile router was used to connect the LasUAV computer, the RTK GNSS receiver and a laptop ground station, to remotely operate the sensors, perform diagnostic checks and control the acquisition. A Matlab (MathWorks, Inc., Natick, MA, USA) custom software was developed to generate the LiDAR point cloud in LAS (LASer) file format in post-processing. LasUAV was mounted in a 12 brushless motors hexa-coaxial UAV rotor, equipped with a control and navigation system capable of autonomous flight via uploaded waypoints. The net endurance was around 15 min with this payload, contained within a case which was fixed in a rigid isolated frame ( Figure 1). A more detailed and thorough illustration of LasUAV components and point cloud generation performance is reported in [14].

Study Area
The study area was located on Monte Morello (934 m ASL) at N43°52'55.92" and E11°13'53.04" in the NW side of Florence (Italy). In the historic past, this mountain area was affected by an intensive browsing and grazing that altered the forest vegetation [49] and triggered hydrogeological instability. During the last century, the area was therefore subject to reforestation concluded in the early 1980s. Conifer species, in particular European black pine (Pinus nigra J. F. Arnold), Brutia pine (Pinus brutia Ten.), Mediterranean cypress (Cupressus sempervirens L.) and Atlas cedar (Cedrus atlantica (Endl.) Carrière), were chosen for their capability to colonize bare-earth with occasional broadleaves such as Turkey oak (Quercus cerris L.) and leaving already existing native trees such as downy oak (Quercus pubescens Willd.). Nowadays, the stands consist of two-layered mixed even-aged high forest. The upper layer was composed of dominant trees where coniferous species and oaks were prevailing; the bottom layer consisted of overtopped trees where the presence of manna ash (Fraxinus ornus L.), field maple (Acer campestre L.) and hawthorn (Crataegus monogyna Jacq.) was significant in terms of number of trees, but less significant in terms of biomass. Table 1 summarizes the characteristics of the study area.

Study Area
The study area was located on Monte Morello (934 m ASL) at N43 • 52'55.92" and E11 • 13'53.04" in the NW side of Florence (Italy). In the historic past, this mountain area was affected by an intensive browsing and grazing that altered the forest vegetation [49] and triggered hydrogeological instability. During the last century, the area was therefore subject to reforestation concluded in the early 1980s. Conifer species, in particular European black pine (Pinus nigra J. F. Arnold), Brutia pine (Pinus brutia Ten.), Mediterranean cypress (Cupressus sempervirens L.) and Atlas cedar (Cedrus atlantica (Endl.) Carrière), were chosen for their capability to colonize bare-earth with occasional broadleaves such as Turkey oak (Quercus cerris L.) and leaving already existing native trees such as downy oak (Quercus pubescens Willd.). Nowadays, the stands consist of two-layered mixed even-aged high forest. The upper layer was composed of dominant trees where coniferous species and oaks were prevailing; the bottom layer consisted of overtopped trees where the presence of manna ash (Fraxinus ornus L.), field maple (Acer campestre L.) and hawthorn (Crataegus monogyna Jacq.) was significant in terms of number of trees, but less significant in terms of biomass. Table 1 summarizes the characteristics of the study area.  Figure 2). Three plots, positioned in the northeast side of the study area, were selected for this study.  Figure 2). Three plots, positioned in the northeast side of the study area, were selected for this study. Field surveys were carried out from June to September 2016, in leaf-on period. For each tree with DBH bigger or equal to 3.5 cm, the species, the health state (unharmed, broken or dead) and the DBH were recorded. The polar coordinates of each tree were collected using a TruPulse Laser Rangefinder 200/B (Laser Technology, Inc., Centennial, USA) and then converted into Cartesian coordinate with origin in the plot center. For a subsample of 10 trees per plot, the height was also measured with the rangefinder-hypsometer Haglöf Vertex IV (Haglöf Sweden AB, Sverige, Sweden). These trees were selected to distribute height measures along the DBH range for each plot and to measure only trees with very visible tops. To estimate the height of the trees for which only DBH was measured in the field, a bivariate model was fitted with the function aov() of R [50] using log(DBH) as a continuous predictor, tree species as a categorical predictor and measured tree heights as the dependent variable The model then estimates tree height (H) as: Field surveys were carried out from June to September 2016, in leaf-on period. For each tree with DBH bigger or equal to 3.5 cm, the species, the health state (unharmed, broken or dead) and the DBH were recorded. The polar coordinates of each tree were collected using a TruPulse Laser Rangefinder 200/B (Laser Technology, Inc., Centennial, USA) and then converted into Cartesian coordinate with origin in the plot center. For a subsample of 10 trees per plot, the height was also measured with the rangefinder-hypsometer Haglöf Vertex IV (Haglöf Sweden AB, Sverige, Sweden). These trees were selected to distribute height measures along the DBH range for each plot and to measure only trees with very visible tops. To estimate the height of the trees for which only DBH was measured in the field, a bivariate model was fitted with the function aov() of R [50] using log(DBH) as a continuous predictor, tree species as a categorical predictor and measured tree heights as the dependent variable The model then estimates tree height (H) as: where a and b are the model intercept and slope, and DBH is the stem diameter at 1.3 m height.
[Species] is the categorical predictor on which the model builds the contrast matrix and is composed by the following 8 categories: (i) Acer campestre, (ii) Cupressus sempervirens, (iii) Fraxinus ornus, (iv) Pinus nigra, Volume for all live, unbroken trees was computed using the equations developed by [51] for the second Italian national forest inventory which allows the prediction of the total above ground biomass using DBH and total tree height as independent variables for twenty-five of the most important forest species growing in Italy.
A univariate analysis of variance (ANOVA), followed by a multiple comparison of tree height means, was used to identify statistically significant differences between the biometric parameters of each species, giving insight on species-specific effects on the algorithms' output. Both the ANOVA and the multiple comparison test were part of the Statistics and Machine Learning Toolbox of MATLAB R2017a (functions anova1 and multcompare, MathWorks, Inc., Natick, MA, USA).
The dendrometric features of each plot are reported in Table 2. In plot 6.1, 24.2% of total trees were present in the canopy layer and the remaining 75.8% in the understory layer. In plot 6.2, 43.8% of total trees were present in the canopy layer and the remaining 56.3% in the understory layer. In plot 7.2, 64.8% of total trees were present in the canopy layer and the remaining 35.2% in the understory layer.

LiDAR Data Acquisition and Point Cloud Pre-processing
LiDAR data acquisition was carried out on 28 July 2017. The flight altitude was on average 50 m above ground level, the flying speed was 2.5 m/s and the flight sampling time was 10 min. A total of 8 flight lines, which resulted in a flight line spacing of about 20 m, were necessary to cover the 155 m × 102 m area where plots 6.1, 6.2 and 7.2 were located. Based on a ±30 • scanning angle, the average point density was about 193 points/m 2 . A summary of the characteristics of the UAV flight and LiDAR data collected is reported in Table 3. The whole point cloud gathered by the LiDAR is shown in Figure 3a, while a slice of the cloud highlighting the two-layered structure of the canopy is shown in Figure 3b.  After full point cloud production, three subsets were extracted as cylinders with a radius of 16 m centered at each plot. The choice of the radius larger than the field plot size (13 m) was made to avoid "edge effects" that would have produced the loss of data points belonging to trees with stem inside the plot but crown partially outside. On the other hand, trees with the stem outside and crown partially inside plot boundaries were excluded after their segmentation based on their location outside the plot border.

Individual Tree Crown Delineation Approach and Trees Matching
This approach was translated by [52] in the itcLiDAR() function contained in the itcSegment package which was used in this study. In the itcLiDAR() function, in addition to the x,y,z coordinates of each LiDAR return and the coordinate reference system, the input parameters defined in Table 4 were provided. After full point cloud production, three subsets were extracted as cylinders with a radius of 16 m centered at each plot. The choice of the radius larger than the field plot size (13 m) was made to avoid "edge effects" that would have produced the loss of data points belonging to trees with stem inside the plot but crown partially outside. On the other hand, trees with the stem outside and crown partially inside plot boundaries were excluded after their segmentation based on their location outside the plot border.

Individual Tree Crown Delineation Approach and Trees Matching
This approach was translated by [52] in the itcLiDAR() function contained in the itcSegment package which was used in this study. In the itcLiDAR() function, in addition to the x,y,z coordinates of each LiDAR return and the coordinate reference system, the input parameters defined in Table 4 were provided. Table 4. Input parameters of itcLiDAR() function and their definition. In the definition CHM stands for canopy height model. From the itcLiDAR() function, x,y coordinates, height and crown area of each tree were obtained. The dbh() function of itcSegment was then used to predict the DBH of each delineated tree crown and the agb() function, based on the equations developed by [48], was used to predict AGB of each delineated tree crown. To ensure that the itcLiDAR() function performed optimally for the dense mixed forest under investigation, a tuning of input parameters was made (Table 4). On each plot, a single parameter was varied in a wide interval while maintaining the others (i.e., resolution and cw) fixed according to the ecological meaning of the parameter and the level of knowledge of forest stands characteristics (Table 5). This procedure resulted in a set of 1350 parameter combinations. For each combination, the delineated trees were matched to the field surveyed trees using the [53] method: if only one field-measured tree was included inside a delineated crown, then that tree was associated with the corresponding segmented tree; otherwise if more than one field-measured tree was included in a segmented crown, the field-measured tree with the height closer to the segmented tree height was chosen.
For each plot and for each combination of itcLiDAR() function parameters the ITC delineation accuracy (accuracy index, AI) was computed as: where OE is the omission error, i.e., existing trees not detected; CE is the commission error, i.e., delineated crowns that do not exist in reality. Detection rate (DET) was computed as the ratio between the number of matched trees and the number of trees measured in field. The best set of itcLiDAR() function parameters was the one that maximized AI. At tree level, for each delineated tree crown, DBH and height measured in the field and volume computed from the second Italian national forest inventory equations, were compared to the corresponding variables resulting from the UAV point cloud processing with the itcSegment() package. A comparison between the DBH frequency distribution from both datasets was also performed.
At plot level, for the entire set of matched trees a comparison between average tree density, DBH, height and AGB was made. The AGB estimate uncertainty was calculated as the percentage difference between AGB estimated from the UAV-LS data and from second Italian national forest inventory equations.
For comparison, the approach of individual tree crowns delineation developed by [47] was accomplished. The algorithm, developed for small footprint discrete return ALS point cloud, segments individual trees in a sequence by taking advantage of the relative spacing between trees. Starting from a treetop, the algorithm identifies and "grows" a target tree by including nearby points and exclude points of other trees based on their relative spacing. To overcome the difficulty of classifying points at lower level due to the diminished spacing between trees, points were classified sequentially, from the highest to the lowest. Points with spacing larger than a specified threshold were excluded from the target tree; points with spacing smaller than the threshold were classified based on a minimum spacing rule [47]. This segmentation method is contained in the li2012 algorithm in the lastrees() function of lidR package [54]. A parameter tuning procedure similar to the one used for tuning the itcLiDAR() function was applied, varying the value of li2012 algorithm parameters in the ranges reported in Table 6.

Calibration of itcLiDAR() Function Parameters
For each of the three plots, the optimal combination of itcLiDAR() function parameters obtained with the iterative tuning procedure is reported in Table 7. The value of maximum size in pixels of the moving window used to detect the local maxima of the CHM (MaxSearchFilSize) was the same in the three plots. The minimum size was equal in plots 6.1 and 6.2, while changed in plot 7.2 that was characterized by a lower tree density and a higher mean DBH. This means that in stands with mature trees, the value of the minimum size of the moving window was close to the maximum, due to a lower crown size variability. Similarly, both the maximum crown diameter of a detected tree (maxDIST) and the first threshold of the minimum vertical distance between the height of each neighbor pixel extracted from the CHM and the height of the local maximum (TRESHSeed) were the same in plot 6.1 and 6.2 and different in plot 7.2. The second threshold of the minimum vertical distance between the height of each neighbor pixel extracted from the CHM and the height of the local maximum (TRESHCrown) and the minimum value of the crown diameter of a detected tree (minDIST) had different values in each plot after the tuning process (Table 7).

ITC Detection and Dendrometric Variables Estimation with the itcLiDAR() Function
The known position of individual trees over imposed to the crowns delineated by the itcLiDAR() function is reported in Figure 4.

Calibration of itcLiDAR() Function Parameters
For each of the three plots, the optimal combination of itcLiDAR() function parameters obtained with the iterative tuning procedure is reported in Table 7. The value of maximum size in pixels of the moving window used to detect the local maxima of the CHM (MaxSearchFilSize) was the same in the three plots. The minimum size was equal in plots 6.1 and 6.2, while changed in plot 7.2 that was characterized by a lower tree density and a higher mean DBH. This means that in stands with mature trees, the value of the minimum size of the moving window was close to the maximum, due to a lower crown size variability. Similarly, both the maximum crown diameter of a detected tree (maxDIST) and the first threshold of the minimum vertical distance between the height of each neighbor pixel extracted from the CHM and the height of the local maximum (TRESHSeed) were the same in plot 6.1 and 6.2 and different in plot 7.2. The second threshold of the minimum vertical distance between the height of each neighbor pixel extracted from the CHM and the height of the local maximum (TRESHCrown) and the minimum value of the crown diameter of a detected tree (minDIST) had different values in each plot after the tuning process (Table 7).

ITC Detection and Dendrometric Variables Estimation with the itcLiDAR() Function
The known position of individual trees over imposed to the crowns delineated by the itcLiDAR() function is reported in Figure 4.   Detection rate ranged between 40% to 57%, omission error between 42% and 59% and commission error between 41% and 53% ( Table 8). The highest detection rate and overall accuracy was obtained in plot 7.2, characterized by a lower tree density and with 65% of trees (composed by European black pine, Brutia pine and Mediterranean cypress) in the upper layer and 35% in the understory. This indicates that plot 7.2 was easier to segment due to the shape of conifer crowns that allowed the clear identification of the apex. Quite high levels of under segmentation were obtained in plot 6.1 and 6.2 where the overlap between close trees was relevant, being the plot density very high. Table 8. ITC delineation accuracy for three plots (DET = detection rate, crown correctly detected, OE = omission error, failure to detected a crown that exists, CE = commission error, delineation of a crown that do not exists in reality, AI = accuracy index, 100 -(OE + CE). In plot 6.1 and 6.2, DBH of small trees was estimated successfully, while DBH of medium and large trees was underestimated by the segmentation process. In these two plots, tree density was high because of a large presence of small trees in the understory (75% and 56% in plot 6.1 and 6.2 respectively), but at the same time there were canopy gaps that facilitated the identification of individual crowns ( Figure 4). Instead, in plot 7.2, DBH of small trees was over-estimated while DBH of medium and large trees was under-estimated ( Figure 5).
Drones 2020, 4, x FOR PEER REVIEW 11 of 20 Detection rate ranged between 40% to 57%, omission error between 42% and 59% and commission error between 41% and 53% ( Table 8). The highest detection rate and overall accuracy was obtained in plot 7.2, characterized by a lower tree density and with 65% of trees (composed by European black pine, Brutia pine and Mediterranean cypress) in the upper layer and 35% in the understory. This indicates that plot 7.2 was easier to segment due to the shape of conifer crowns that allowed the clear identification of the apex. Quite high levels of under segmentation were obtained in plot 6.1 and 6.2 where the overlap between close trees was relevant, being the plot density very high. Table 8. ITC delineation accuracy for three plots (DET = detection rate, crown correctly detected, OE = omission error, failure to detected a crown that exists, CE = commission error, delineation of a crown that do not exists in reality, AI = accuracy index, 100 -(OE + CE). In plot 6.1 and 6.2, DBH of small trees was estimated successfully, while DBH of medium and large trees was underestimated by the segmentation process. In these two plots, tree density was high because of a large presence of small trees in the understory (75% and 56% in plot 6.1 and 6.2 respectively), but at the same time there were canopy gaps that facilitated the identification of individual crowns ( Figure 4). Instead, in plot 7.2, DBH of small trees was over-estimated while DBH of medium and large trees was under-estimated ( Figure 5). Figure 5. Estimation of the diameter at breast height (DBH) for the field-measured trees (measured DBH, y-axis) by the segmentation process (estimated DBH, x-axis) in plots 6.1 (a), 6.2 (b) and 7.2 (c). Red lines and associated statistics show linear fits (with the intercept forced to zero) considering the whole dataset, while green lines and statistics show the linear fit after removing the cluster of shorter species (field maple, hawthorn and manna ash). The blue dashed line indicates the 1:1 reference line in all the three plots.

Number of Trees in Field
Individual tree crown delineation was successful for trees with DBH larger than 10 cm (detection rate of 78.3%, 50.0% and 67.6% for plot 6.1, 6.2 and 7.2 respectively), while failed to detect trees with DBH smaller than 10 cm (detection rate of 19.5%, 33.5% and 36.8% for plot 6.1, 6.2 and 7.2 respectively). The algorithm over-segmented trees with DBH between 16 and 24 cm, while undersegmented trees with DBH between 28 and 36 cm ( Figure 6). Figure 5. Estimation of the diameter at breast height (DBH) for the field-measured trees (measured DBH, y-axis) by the segmentation process (estimated DBH, x-axis) in plots 6.1 (a), 6.2 (b) and 7.2 (c). Red lines and associated statistics show linear fits (with the intercept forced to zero) considering the whole dataset, while green lines and statistics show the linear fit after removing the cluster of shorter species (field maple, hawthorn and manna ash). The blue dashed line indicates the 1:1 reference line in all the three plots.
Individual tree crown delineation was successful for trees with DBH larger than 10 cm (detection rate of 78.3%, 50.0% and 67.6% for plot 6.1, 6.2 and 7.2 respectively), while failed to detect trees with DBH smaller than 10 cm (detection rate of 19.5%, 33.5% and 36.8% for plot 6.1, 6.2 and 7.2 respectively). The algorithm over-segmented trees with DBH between 16 and 24 cm, while under-segmented trees with DBH between 28 and 36 cm ( Figure 6).
In all the three plots, tree heights around 10 m and in the range between 15 to 20 m were quite well estimated, while in the range between 10 to 15 m they were overestimated by about 50%. In terms of tree species, height overestimation especially affects species that were significantly shorter (field maple, hawthorn, manna ash, p < 0.05 for the ANOVA test) and that clearly cluster separately from the other species (Figure 7). Linear fitting between measured and derived heights (red lines and statistics,  In all the three plots, tree heights around 10 m and in the range between 15 to 20 m were quite well estimated, while in the range between 10 to 15 m they were overestimated by about 50%. In terms of tree species, height overestimation especially affects species that were significantly shorter (field maple, hawthorn, manna ash, p < 0.05 for the ANOVA test) and that clearly cluster separately from the other species (Figure 7). Linear fitting between measured and derived heights (red lines and statistics, Figure 7a  The correlation between tree species and height overestimation also affects RMSE between measured and estimated height when pooling the data by species rather than by plot (Figure 8). The RMSE shows a strong inverse correlation with the average measured species height. Field maple (with an average height of 9.37±2.31 m) has the highest RMSE (9.75 m), while the Atlas cedar (single exemplar with a measured height of 20.2 m) has the lowest RMSE (0.07 m).  In all the three plots, tree heights around 10 m and in the range between 15 to 20 m were quite well estimated, while in the range between 10 to 15 m they were overestimated by about 50%. In terms of tree species, height overestimation especially affects species that were significantly shorter (field maple, hawthorn, manna ash, p < 0.05 for the ANOVA test) and that clearly cluster separately from the other species (Figure 7). Linear fitting between measured and derived heights (red lines and statistics, Figure 7a  The correlation between tree species and height overestimation also affects RMSE between measured and estimated height when pooling the data by species rather than by plot (Figure 8). The RMSE shows a strong inverse correlation with the average measured species height. Field maple (with an average height of 9.37±2.31 m) has the highest RMSE (9.75 m), while the Atlas cedar (single exemplar with a measured height of 20.2 m) has the lowest RMSE (0.07 m). The correlation between tree species and height overestimation also affects RMSE between measured and estimated height when pooling the data by species rather than by plot (Figure 8). The RMSE shows a strong inverse correlation with the average measured species height. Field maple (with an average height of 9.37±2.31 m) has the highest RMSE (9.75 m), while the Atlas cedar (single exemplar with a measured height of 20.2 m) has the lowest RMSE (0.07 m).
The total above ground biomass estimated for each tree was compared with field estimates in Figure 9. An overall bias was obtained, with the biomass of small trees resulting overestimated and that of large trees underestimated. Table 9 reports a comparison between average value of tree density, DBH, height and AGB at plot level. Considering the average values of above ground biomass measured in the field and estimated by the individual tree delineation process, the biomass estimates uncertainty resulted about 13%, 41% and −9% for plot 6.1, 6.2 and 7.2 respectively. The algorithm, instead, tends to overestimate both mean height and mean DBH (with the exception of DBH in plot 7.2). The total above ground biomass estimated for each tree was compared with field estimates in Figure 9. An overall bias was obtained, with the biomass of small trees resulting overestimated and that of large trees underestimated. Figure 9. Estimation of the above ground biomass (estimated AGB, x-axis) for the field-measured trees (measured AGB, y-axis) in plots 6.1 (a), 6.2 (b) and 7.2 (c). Red lines and associated statistics show linear fits (with the intercept forced to zero) considering the whole dataset, while green lines and statistics show the linear fit after removing the cluster of shorter species (field maple, hawthorn and manna ash). The blue dashed line indicates the 1:1 reference. Table 9 reports a comparison between average value of tree density, DBH, height and AGB at plot level. Considering the average values of above ground biomass measured in the field and estimated by the individual tree delineation process, the biomass estimates uncertainty resulted about 13%, 41% and -9% for plot 6.1, 6.2 and 7.2 respectively. The algorithm, instead, tends to overestimate both mean height and mean DBH (with the exception of DBH in plot 7.2). Table 9. Tree number, mean DBH, mean height, total above ground biomass (AGB) estimation field surveyed (= FS) and derived from itcLiDAR() function (= UAV). The total above ground biomass estimated for each tree was compared with field estimates in Figure 9. An overall bias was obtained, with the biomass of small trees resulting overestimated and that of large trees underestimated. Figure 9. Estimation of the above ground biomass (estimated AGB, x-axis) for the field-measured trees (measured AGB, y-axis) in plots 6.1 (a), 6.2 (b) and 7.2 (c). Red lines and associated statistics show linear fits (with the intercept forced to zero) considering the whole dataset, while green lines and statistics show the linear fit after removing the cluster of shorter species (field maple, hawthorn and manna ash). The blue dashed line indicates the 1:1 reference. Table 9 reports a comparison between average value of tree density, DBH, height and AGB at plot level. Considering the average values of above ground biomass measured in the field and estimated by the individual tree delineation process, the biomass estimates uncertainty resulted about 13%, 41% and -9% for plot 6.1, 6.2 and 7.2 respectively. The algorithm, instead, tends to overestimate both mean height and mean DBH (with the exception of DBH in plot 7.2). Table 9. Tree number, mean DBH, mean height, total above ground biomass (AGB) estimation field surveyed (= FS) and derived from itcLiDAR() function (= UAV).

Plot
FS UAV FS UAV FS UAV FS UAV Figure 9. Estimation of the above ground biomass (estimated AGB, x-axis) for the field-measured trees (measured AGB, y-axis) in plots 6.1 (a), 6.2 (b) and 7.2 (c). Red lines and associated statistics show linear fits (with the intercept forced to zero) considering the whole dataset, while green lines and statistics show the linear fit after removing the cluster of shorter species (field maple, hawthorn and manna ash). The blue dashed line indicates the 1:1 reference. Table 9. Tree number, mean DBH, mean height, total above ground biomass (AGB) estimation field surveyed (= FS) and derived from itcLiDAR() function (= UAV).

Calibration of Function Parameters, ITC Detection and Dendrometric Variables Estimation with the li2012 Algorithm in Lastrees() Function
The optimal set of parameters of the li2012 algorithm after the iterative tuning procedure for each plot is reported in Table 10.  Figure 10 reports the segmented crowns of trees with respect to the position of tree stems surveyed in field, showing a clear tendency of the algorithm to aggregate the crowns of more than one tree and consequently delineate bigger crowns than in reality.

Calibration of Function Parameters, ITC Detection and Dendrometric Variables Estimation with the li2012 Algorithm in Lastrees() Function
The optimal set of parameters of the li2012 algorithm after the iterative tuning procedure for each plot is reported in Table 10.  Figure 10 reports the segmented crowns of trees with respect to the position of tree stems surveyed in field, showing a clear tendency of the algorithm to aggregate the crowns of more than one tree and consequently delineate bigger crowns than in reality. Values of detection rate are quite far from those obtained with the application of itcLiDAR, varying between 25% and 27% (Table 11). In case of plot 6.1 a value of commission error equal to 0% was obtained because no segmented trees were not matched: indeed, as the obtained crowns had a large diameter there were always one or more field measured tree inside the delineated crown.  Values of detection rate are quite far from those obtained with the application of itcLiDAR, varying between 25% and 27% (Table 11). In case of plot 6.1 a value of commission error equal to 0% was obtained because no segmented trees were not matched: indeed, as the obtained crowns had a large diameter there were always one or more field measured tree inside the delineated crown. As the level of accuracy reached by li2012 algorithm in lastrees() was very low, no details in the performance DBH, height and AGB are provided and the following discussion is therefore focused on the best-performing algorithm.

Discussion
Delineation accuracy of ITC segmentation obtained in this study was in line with other studies based on ALS. Ref. [55] evaluated the quality, accuracy and feasibility of nine automatic tree extraction methods based on ALS data acquired in boreal forest stands with point densities of 2, 4 and 8 points/m 2 . The capabilities of detecting trees resulted between 25% and 102%, while the point density had a minor impact on the segmentation rate. The commission error varied between 0% and 26.5%, while the omission error varied between 35% and 68%. As expected, the taller trees were detected with a better accuracy; the RMSE of height estimation for the delineated trees ranged from 60 to 80 cm for the best methods. Ref. [45] benchmarked eight ALS tree detection methods, including the method that was used here [53], across different regions, forest types and structures in the Alps. Considering all the methods, the omission error ranged between 51% and 63%, being higher for multi-layered mixed and coniferous forests and lower for single-layered coniferous forest (37%). The method of [53] that was used here produced an overall omission error of 61%, while the omission error in single layered mixed forests was 25%, in single layered coniferous forests was 40%, in multi layered mixed forests was 75% and in multi layered coniferous forest was 65%. Ref. [56] benchmarked five ALS-based individual tree detection methods with a point cloud density of 2, 4 and 8 points/m 2 acquired in five boreal forests with diverse stand conditions and different developmental stages (i.e., from deciduous-dominated to mixed, to coniferous-dominated and from old to mature to very young stands) considering four crown classes, i.e., dominant, codominant, intermediate and suppressed trees. Their study revealed that omission error varied between 25% and 63% and that all methods tended to overestimate the number of dominant trees and underestimate the number of subordinate trees. Delineation of dominated trees was challenging for all methods.
Similar to [45,55,56], one finding of our study was that the omission error varies between 40% and 57% and that the largest fraction of not detected trees was in the understory, which suggests they were dominated trees. This outcome, that needs to be confirmed by further studies considering larger datasets, confirms what [55] highlighted concerning the effect of point density on the ITC delineation accuracy and retract the outcomes of [56]: high levels of laser point density do not increase the delineation accuracy. In fact, also using LiDAR data with very high density as in this study (i.e., 193 points/m 2 ) the accuracy of tree segmentation was similar to that obtained in the studies that benchmarked individual tree detection methods using ALS data with point density ranging from 2 to 8 points/m 2 .
Ref. [15] used high-resolution LiDAR data captured from a small multirotor UAV to extract individual trees in 6 circular plots in a four-year-old Tasmanian blue gum regular plantation (density ranging from 680 to 1560 trees/ha) located in Southern Tasmania (Australia) using an adaptation of the method outlined by [47]. By means of this approach, 91% of field-measured trees were correctly matched to the UAV segmented trees across all plots. In their successive study, Ref. [57], using the same data set, tested five tree detection routines to assess the tree delineation accuracy. The percentage of field-measured trees correctly matched with UAV-LS delineated trees resulted between 92% and 97%. The highest rate of omission occurred in small trees, while commission errors for most algorithms was caused by over-segmentation of large trees. The study also highlighted that improvement in omission rate occurs using higher point densities (~50 points/m 2 ).
Ref. [36] delineated trees from UAV-borne laser scanner data acquired over a 64 m × 64 m study area located in Evo, Southern Finland, with a density of about 490 trees/ha. A two steps approach for individual tree detection was used, based on the application of a local maximum filtering algorithm on the CHM and on the application of a watershed algorithm to successively delineate the tree crowns. The average detection rate resulted 84.2%, being 100% for isolated and dominant trees, but the relatively low value of tree density of the study area, facilitating the detection process, should be considered. When utilizing direct measurement in the point cloud, Ref. [36] obtained a relative RMSE equal to 6.47% in height estimation and 27.46% in DBH estimation.
Ref. [38] delineated individual tree crows from UAV-LS data acquired in an ecotone area in Northern Arizona (USA) which includes a shrubland-grassland meadow adjacent to a ponderosa pine forest with a density of 220 trees/ha. They obtained a determination coefficient between the number of trees delineated from the LiDAR data and field-based tree density equal to 0.77 and a determination coefficient for tree height within dense canopy cover equal to 0.79, where only the tallest tree canopy heights were successfully estimated and the shorter trees under the tall canopies were missed by the laser scanner.
Ref. [58] applied a methodology to UAV-LS data based on fitting and modeling the 3D-shape of trees in a hazel grove using RANSAC [59] obtaining a detection rate equal to 86%.
Ref. [60] assessed the performance of four algorithms (three based on CHM and one on point cloud) for individual tree detection using UAV-LS data in pure ginkgo planted forest in China with density ranging from 400 to 700 trees per hectare, DBH ranging from 16 cm to 20 cm and Lorey's height from 12 m to 11 m. The individual tree detection algorithm based on point cloud was the best performing in plot with low density but was the worst as the tree density increased.
Unlike our study, that was conducted in multi-layer forest stands with a density higher than 1000 trees/ha and with a high level of crown cover, the results of the studies reported above were obtained in isolated trees, in stands composed by only dominated trees, in plantation where trees were planted with a regular pattern, and overall, in stands with a tree density approximately half of that of our analysis. And as a matter of fact, the levels of accuracy obtained in our study were slightly lower. The optimal setting of itcLiDAR() function has determined a detection rate ranging between 40% and 57% in the three plots, while the main source of omission error was the failure to separate two or three medium trees that have similar heights and were located close to each other. The tuning of the ITC algorithm parameters resulted quite important in reaching the highest determination rates that ranged from 23.4% to 40.6% in plot 6.1, 21.2% and 40.9% in plot 6.2 and 32.1% and 57.1% in plot 7.2 across the parameter optimization space. This was an important outcome: knowing the possible range of variability of each parameter was the starting point for the individual tree segmentation process to reach the optimal determination rates. Hence, the application of a parameter tuning procedure is suggested to optimize the segmentation process for a specific set of forest characteristics. Our findings agree with [44], which reports that a certain degree of calibration is necessary to obtain the best results in tree detection and delineation, where size and shape of the 3D kernel play a crucial role. Another reason that explains the lower values of detection rate obtained here with respect to other studies is the approach used to match the delineated trees with the field surveyed trees. The approach of [53] considers the inclusion of a field-measured tree inside the delineated crown and difference of height between the field measured tree and the delineated tree. Other methods consider also the planar distance between delineated and reference trees (i.e., [61]) or the maximum values of height differences (i.e., [62]). Depending on the method used for tree matching, the results can be quite different. Vertical forest stratification in combination with the matching process were probably the factors that have determined the overestimation of canopy height. Intermediate trees, that occupied a position underneath the dominant and co-dominants trees, were erroneously coupled with the higher trees. Indeed, in some cases crowns of dominant trees, due to the specific shape of species present in the stands, were segmented as two distinct trees, so in this way a delineated tree was coupled with the dominant tree and the other one with the codominant or intermediate tree, leading to an over estimation of tree height and consequently of tree biomass.
The effect of stratification on the estimation of tree heights was clearly visible when removing shorter trees from the analysis. The coefficient of determination for estimated versus measured tree height, in fact, greatly increased when significantly shorter trees were removed, reaching even higher values than [38] in plot 6.1 (0.90 vs. 0.79). In fact, when considering individually the single species, the effect of stratification was relevant: taller species that stand out clearly from the dense mixed canopy have lower RMSEs between measured and estimated heights, while shorter species were less easily distinguished and have associated higher errors. We should also consider the error that occur during the ground-based tree height measurements which was related to the tree top shape, determined by the apical dominance: the higher the apical dominance, the greater the ability to detect the tree top in field survey. This implies that in species with a strong apical dominance (e.g., Atlantic cedar and Mediterranean cypress) the uncertainty in tree height field measurement was lower than in species with moderate or low apical dominance (e.g., European black pine and Turkey oak). This aspect reflects in the RMSE that was lower in species with greater apical dominance (Figure 8).
This study was based on a sample size of 186 trees, which could prevent the upscaling of the results obtained here to larger areas. However, comparable sample sizes were used in several studies focusing on small areas without loss of generality [58]. Ref. [43] employed UAV-LS data to build 200 tree models by a quantitative structure approach and compared them with terrestrial laser scanning obtaining reliable volume estimations for trunks and branches ≥ 30 cm. Similarly, Ref. [44] used multiple segmentation algorithms over about 120 trees, while [58] performed 3D volume reconstruction with UAV-LS on 181 hazel trees. Overall, these studies increase the confidence in the validity of the approach proposed here, that definitely calls for further and deeper analysis based on forthcoming larger datasets.

Conclusions
This study explored the application of two different algorithms, i.e., itcLiDAR() and li2012, for individual tree crown segmentation in a two-layered dense mixed forest. The parameters of both algorithms were tuned for this experimental site through an optimization procedure. The raster-based local maxima region growing algorithm itcLiDAR() performed better than the point cloud-based li2012 and the performance of its estimates resulted comparable to those obtained in ITC segmentation based on ALS data. The parameters calibration for individual tree crown detection was necessary and allowed to increase the performance of the segmentation process. The accuracy obtained with the method based on CHM maxima was reasonable for forest inventory, but still not satisfactory. In fact, such method proved to be very efficient for forests composed by conifers with a well-defined apex (i.e., European spruce, fir, European larch) and with uniform height implying they were not multilayered. In presence of broadleaved forests, forests composed by trees having round, broad, spreading and irregular crowns-or in the case of mixed multilayered forests-the segmentation procedure becomes challenging. The exploration of algorithms-other than those based on generating a CHM and searching for local maxima in the top of the canopy indicating the position of the tree-is needed. These should allow achieving higher detection rates and accuracy on dendrometric variables estimated across a wide range of species and in multilayered forests where most of trees were hidden under or between taller and larger trees. Given the complex and heterogeneous structure of many forests it is important that future research is conducted in this direction, assuming that the main advantage of individual tree detection stays in the capacity to provide true stem distribution series, enabling accurate estimates of timber assortments.