A Hybrid Method for Individual Tree Detection in Broadleaf Forests Based on UAV-LiDAR Data and Multistage 3D Structure Analysis

: Individual tree detection and segmentation in broadleaf forests have always been great challenges due to the overlapping crowns, irregular crown shapes, and multiple peaks in large crowns. Unmanned aerial vehicle (UAV)-borne light detection and ranging (LiDAR) is a powerful tool for acquiring high-density point clouds that can be used for both trunk detection and crown segmentation. A hybrid method that combines trunk detection and crown segmentation is proposed to detect individual trees in broadleaf forests based on UAV-LiDAR data. A trunk point distribution indicator-based approach is first applied to detect potential trunk positions. The treetops extracted from a canopy height model (CHM) and the crown segments obtained by applying a marker-controlled watershed segmentation to the CHM are used to identify potentially false trunk positions. Finally, the three-dimensional structures of trunks and branches are analyzed at each potentially false trunk position to distinguish between true and false trunk positions. The method was evaluated on three plots in subtropical urban broadleaf forests with varying proportions of evergreen trees. The F-score in three plots ranged from 0.723 to 0.829, which are higher values than the F-scores derived by a treetop detection method (0.518–0.588) and a point cloud-based individual tree segmentation method (0.479–0.514). The influences of the CHM resolution (0.25 and 0.1 m) and the data acquisition season (leaf-off and leaf-on) on the final individual tree detection result were also evaluated. The results indicated that using the CHM with a 0.25 m resolution resulted in under-segmentation of crowns and higher F-scores. The data acquisition season had a small influence on the individual tree detection result when using the hybrid method. The proposed hybrid method needs to specify parameters based on prior knowledge of the forest. In addition, the hybrid method was evaluated in small-scale urban broadleaf forests. Further research should evaluate the hybrid method in natural forests over large areas, which differ in forest structures compared to urban forests.


Introduction
Precise tree inventory plays a critical role in sustainable forest planting, restoration, and management [1].Remote sensing techniques provide powerful tools for collecting tree information, such as location, tree height, crown width, and biomass [2].Airborne light detection and ranging (LiDAR), as an active remote sensing technique that acquires data from above the tree canopy, can capture the three-dimensional (3D) structure information of forests and has been widely used for collecting precise tree information [3,4].The basis of precise tree inventory involves detecting individual trees and delineating their crown boundaries.Due to the complex structures of natural or urban forests, individual tree detection and segmentation from airborne LiDAR data have always been a great challenge [5,6].The individual tree detection and segmentation methods can be divided into three categories: raster-based, point cloud-based, and hybrid methods.Raster-based methods usually use a certain layer of LiDAR returns to generate a two-dimensional (2D) raster, such as a canopy height model (CHM), and identify or segment individual trees from the raster by applying a local maximum filter [7,8], hill-climbing method [9], or watershed segmentation [10,11].Despite the simplicity and efficiency of raster-based methods, the generation of a raster neglects the 3D information of point clouds, and the trees partially or totally covered by the overstory trees cannot be detected by the raster-based methods [12,13].
To overcome the limitations of raster-based methods, point cloud-based methods, which can take advantage of the 3D information of point clouds, were proposed.Considering the spatial distribution of laser points within the cone-shaped crowns of coniferous trees, a distance-based clustering method was proposed by Li et al. [14] to segment individual trees directly from point clouds.This method requires the specification of distance thresholds, which are difficult to specify in heterogeneous forests [15].Under the assumption of a higher point density at the crown center, density-based algorithms, such as mean shift, were used to cluster individual tree points [16,17].Shao et al. proposed a two-phase tree detection approach [1].They applied horizontal mean shift clustering in the horizontal phase to locate possible tree positions and verified each tree cluster in the vertical phase by applying thresholds to the gap in the vertical profile and the vertical length ratio of the tree candidate clusters.However, the methods based on the point density assumption are usually ineffective when the point density in the conjunction areas of overlapping crowns is also high [18].Other point cloud-based methods include the layer stacking method [19] and the spectral clustering method [12].
The abovementioned point cloud-based methods focus on the segmentation of tree crowns, while some other methods adopt a bottom-up strategy, i.e., detecting tree trunks first and then segmenting individual trees based on the detected tree trunks.Such a strategy requires adequate point density in the understory layer to separate the trunks from each other [20].Due to the limitation in the point density of trunks, the bottom-up methods were initially developed for terrestrial LiDAR data, which are collected using a bottom-up scanning mode.The density-based spatial clustering of applications with noise (DBSCAN) algorithm is usually used to cluster trunk points from a horizontal slice of normalized LiDAR point cloud at a certain height [21].
After rapid development in recent years, the technology of unmanned aerial vehicle (UAV)-borne LiDAR has been widely used for individual tree detection and segmentation due to its ability to obtain high-density point clouds.The high point density (>100 points/m 2 or even >1000 points/m 2 ) allows for adequate laser pulses to penetrate the canopy and reflect from the trunk.Therefore, bottom-up methods were developed for the UAV-LiDAR data to detect tree trunks under the assumption that the LiDAR points of tree trunks should be distributed along the vertical direction.Shendryk et al. [22] detected trunks using conditional Euclidean distance clustering and delineated tree crowns through random walk segmentation.Jaskierniak et al. [23] applied a novel watershed clustering within horizontally sliced point clouds and merged the slice-specific clusters into branches, canopy clumps, and trunks through a principal component analysis.Deng et al. [24] proposed a trunk point distribution indicator to detect potential tree trunk positions and applied random sample consistency (RANSAC)-based 3D line fitting to differentiate tree trunks from understory vegetation.Despite the high accuracies achieved by this type of method, the trunk segmentation is subjected to the occlusion problem caused by the dense canopy and tends to be affected by the relatively high understory vegetation around the trunks [24].
Some studies have developed hybrid methods that combine the advantages of the raster-based and point cloud-based methods.The hybrid methods usually segment a CHM first to derive the initial result, followed by a point cloud-based method to refine the initial result [15].Duncanson et al. [11] proposed a two-stage approach that applied a watershed segmentation on a CHM and then generated height histograms within each segment to identify the laser returns from the lower canopy layer so that the understory trees could be detected.Huo et al. [25] took the local maxima extracted from a smoothed CHM as treetop seeds and distinguished between true and false treetop seeds based on the trees' symmetrical structure represented by symmetry curves generated from the point cloud around each seed.
Another type of hybrid method integrates trunk detection with crown segmentation to improve the accuracy of individual tree detection.Reitberger et al. [26] proposed a tree segmentation approach that first applied watershed segmentation on a CHM, then detected tree trunks in each segment using a hierarchical clustering scheme, and finally segmented individual trees using a normalized cut segmentation method.Focusing on poplar plantations, Pu et al. [27] detected trunk seed points from leaf-off LiDAR data using three methods (local maximum, local minimum, and mean shift algorithm), and applied seed-based segmentation methods to segment individual trees from leaf-on LIDAR data.The segmentation result was refined in a horizontal plane based on three tree-crown features.Xu et al. [28] applied Laplacian smoothing and watershed segmentation directly to the point cloud to derive tree clusters and created two segmentations based on treetops and tree bottoms (trunks) within each cluster.Finally, individual trees were identified by intersecting the two segmentations.
In addition to traditional object detection and segmentation methods, deep-learning methods have been proposed to detect trees from point clouds in two-dimensional [29] or three-dimensional space [30].Deep-learning methods provide powerful tools to extract the structure features of trees for individual tree detection.However, training deep-learning models usually needs thousands of (or more) labeled samples.The production of such a sample dataset requires a huge workload.Other deep-learning method challenges include a significant loss of 3D information during the conversion from 3D point clouds into 2D images, and the difficulty in adapting the trained models to forests with diverse species, shapes, and sizes [31].
Although various methods have been developed to identify individual trees from raster data or point clouds, separating overlapping crowns and detecting understory trees remain great challenges [11].In particular, in broadleaf forests, the irregular crown shapes and interlacing branches of neighboring trees bring about more difficulties [22,32,33].Furthermore, the crowns of broadleaf trees usually have multiple peaks and, hence, a tree may be divided into more than one segment, resulting in over-segmentation [34].Therefore, to improve the accuracy of individual tree detection in broadleaf forests, multiple aspects of information should be fully utilized, such as crown morphology and the 3D structure of trunks and branches.The question of how to integrate multiple aspects of information in a hybrid method to improve the accuracy of individual tree detection should be addressed.In the broadleaf forests containing both evergreen and deciduous trees, the dense canopy of evergreen trees obstructs the laser pulses even in the leaf-off season, and the potential of using UAV-LiDAR data to detect individual tree trunks should be explored.Furthermore, trunk detection is usually based on the LiDAR data collected in the leaf-off season (hereafter referred to as leaf-off data) because more laser pulses penetrate the canopy in the leaf-off season [35].On the contrary, the crown segmentation is commonly performed on the LiDAR data acquired in the leaf-on seasons (hereafter referred to as leaf-on data) [27].Since it takes more time to acquire multi-season data, it is necessary to analyze the influence of the data acquisition season on the tree detection result and evaluate the possibility of using the hybrid method to segment crowns and detect individual tree trunks based on single-season data.
The first objective of this work is to develop a hybrid method for detecting individual trees in broadleaf forests based on UAV-LiDAR data.The method combines trunk detection and crown segmentation and takes full advantage of the 3D structure information of branches and trunks.The second objective is to explore the potential of using the hybrid method to detect trees in broadleaf forests containing a relatively high proportion of evergreen trees, whose dense crowns obstruct laser pulses even in the leaf-off season.The third objective is to analyze the influence of data acquisition season (leaf-off and leaf-on) on individual tree detection.The performance of the hybrid method is evaluated on three sample plots in subtropical urban broadleaf forests with varying proportions of evergreen trees and compared with other raster-based and point cloud-based methods.The effect of data acquisition season on individual tree detection is analyzed by utilizing leaf-off and leaf-on UAV-LiDAR data separately for crown segmentation in our hybrid method and the methods for comparison.The remainder of the paper is organized as follows.In Section 2, the study area and the collection of UAV-LiDAR data and reference data are introduced.The proposed hybrid method is detailed in Section 3. The experimental results are given in Section 4 and discussed in Section 5. Finally, we draw some conclusions in Section 6.

Research Materials 2.1. Study Area
The study area is located on the campus of Zhejiang A&F University (119 • 22 ′ 43 ′′ E, 30 • 02 ′ 01 ′′ N) in Lin'an District, Hangzhou, southeastern China (Figure 1).This area has a subtropical monsoon climate and is covered by patches of broadleaf forests.In the study area, three sample plots with varying sizes were established.Their boundaries were delineated in UAV-RGB imagery and the details of the plots are shown in Table 1.
branches and trunks.The second objective is to explore the potential of using the hybrid method to detect trees in broadleaf forests containing a relatively high proportion of evergreen trees, whose dense crowns obstruct laser pulses even in the leaf-off season.The third objective is to analyze the influence of data acquisition season (leaf-off and leaf-on) on individual tree detection.The performance of the hybrid method is evaluated on three sample plots in subtropical urban broadleaf forests with varying proportions of evergreen trees and compared with other raster-based and point cloud-based methods.The effect of data acquisition season on individual tree detection is analyzed by utilizing leaf-off and leaf-on UAV-LiDAR data separately for crown segmentation in our hybrid method and the methods for comparison.The remainder of the paper is organized as follows.In Section 2, the study area and the collection of UAV-LiDAR data and reference data are introduced.The proposed hybrid method is detailed in Section 3. The experimental results are given in Section 4 and discussed in Section 5. Finally, we draw some conclusions in Section 6.

Study Area
The study area is located on the campus of Zhejiang A&F University (119°22′43″ E, 30°02′01″ N) in Lin'an District, Hangzhou, southeastern China (Figure 1).This area has a subtropical monsoon climate and is covered by patches of broadleaf forests.In the study area, three sample plots with varying sizes were established.Their boundaries were delineated in UAV-RGB imagery and the details of the plots are shown in Table 1.Plot P1 is mostly covered by tall broadleaf trees, amongst which some trees, such as Cinnamomum camphora (L.)J. Presl, have relatively large crowns.Their branches intertwine with surrounding trees, and their crown boundaries are difficult to delineate.A few trees are partially or completely covered by neighboring large trees.They are difficult to identify from CHM or RGB images.Some species of trees, such as osmanthus, have dense crowns, relatively low tree heights (3-5 m), and multiple trunks growing together.Their crowns have low bases (1-2 m) and even nearly touch the ground.The proportion of evergreen trees in plot P1 is 56.7%.Plot P2 has a similar tree species composition to plot P1 but contains more (65.2%)evergreen trees.In plot P2, the deciduous trees are mainly willows which are concentrated along the eastern boundary of the plot.Plot P3 is on a gentle slope of approximately 9 degrees.In plot P3, more than half (56.6%) of the trees are deciduous trees.The main deciduous tree species is Choerospondias axillaris (Roxb.)B. L. Burtt & A. W. Hill, which is distributed in the southwest part of the plot.Plot P3 contains a larger proportion of low trees, like osmanthus, than both plots P1 and P2.In all three plots, the understory vegetation is mainly composed of low bushes and herbs.

UAV-LiDAR Data Collection
The UAV-LiDAR data were collected in both leaf-off (February 2023) and leaf-on seasons (May 2023) using a DJI Zenmuse L1 system (DJI Co., Shenzhen, Beijing).The L1 system is mounted on the M300 RTK platform and integrates a Livox LiDAR module, a high-accuracy IMU, and a camera.It can capture point clouds at 240,000 points per second, offering 5 cm vertical and 10 cm horizontal accuracy.The LiDAR data were captured using the three-return mode with a flight altitude of 90 m, a flight speed of 3.5 m per second, and an overlap between strips of 40%.
The DJI Terra v3.1.4software was used to turn the raw data into point clouds.Then the point clouds in the three plots were extracted from the original point cloud in ArcMap 10.6 based on the plot boundaries (Figure 2).The LiDAR point density in each plot ranged from 892 to 1090 points/m 2 .In CloudCompare v2.13.1 software, the noises were removed and each point cloud was classified into ground and non-ground points using the cloth simulation filter [36].The ground points were then used to generate a digital terrain model (DTM) with a 0.5 m resolution, and the non-ground points were normalized using the DTM by subtracting the terrain elevation from the elevation values of each point to remove the topographic effect.

Reference Data Collection
To evaluate the accuracy of individual tree detection by using our method, the reference tree trunk positions were measured using a total station in combination with a realtime kinematic (RTK) receiver.Due to the weak satellite signal under the dense canopy, the RTK was mainly used to establish control points, on which the total station was set up.For each plot, we measured the tree trunk positions in an area larger than the plot, but the trees with trunks outside the plot boundary were not used for the accuracy evaluation.The DBH of each tree was measured using a tape and the trees with DBH less than 6 cm were not included in the procedure of accuracy assessment.We did not measure the heights of all trees but only the minimum and maximum heights in each plot.This is due to the high degree of canopy overlap and the difficulty in identifying treetops.The

Reference Data Collection
To evaluate the accuracy of individual tree detection by using our method, the reference tree trunk positions were measured using a total station in combination with a real-time kinematic (RTK) receiver.Due to the weak satellite signal under the dense canopy, the RTK was mainly used to establish control points, on which the total station was set up.For each plot, we measured the tree trunk positions in an area larger than the plot, but the trees with trunks outside the plot boundary were not used for the accuracy evaluation.The DBH of each tree was measured using a tape and the trees with DBH less than 6 cm were not included in the procedure of accuracy assessment.We did not measure the heights of all trees but only the minimum and maximum heights in each plot.This is due to the high degree of canopy overlap and the difficulty in identifying treetops.The maximum tree height was extracted from the normalized point cloud in each plot.As for the minimum tree height, we selected several low trees and measured their heights using a total station in the field.The height of the lowest tree was taken as the minimum height.The information on the reference data is shown in Table 1.

Methodology
The workflow of our hybrid individual tree detection method (hereafter referred to as the hybrid method) is visualized in Figure 3.The leaf-off UAV-LiDAR data were used to detect potential trunk positions using a trunk point distribution indicator (TPDI)-based approach (Section 3.1).The detected trunk positions were taken as an initial trunk detection result.Then a CHM was generated from the normalized point cloud, smoothed for identification of potential treetops, and segmented using a marker-controlled watershed segmentation (MCWS) algorithm (Section 3.2).The leaf-off LiDAR data acquired in February and the leaf-on data acquired in May were separately used to generate the CHM in order to evaluate the effect of data acquisition season.The potential trunk positions were overlaid on the crown segments so that the segments containing more than one trunk position were identified.Within such segments, the false trunk positions were filtered out by analyzing 3D structure features, i.e., the distribution of LiDAR points in the 3D space (Section 3.3).The methods for comparison with the hybrid method are given in Section 3.4.The accuracy assessment method is explained in Section 3.5.
Forests 2024, 15, x FOR PEER REVIEW 7 of 22 maximum tree height was extracted from the normalized point cloud in each plot.As for the minimum tree height, we selected several low trees and measured their heights using a total station in the field.The height of the lowest tree was taken as the minimum height.
The information on the reference data is shown in Table 1.

Methodology
The workflow of our hybrid individual tree detection method (hereafter referred to as the hybrid method) is visualized in Figure 3.The leaf-off UAV-LiDAR data were used to detect potential trunk positions using a trunk point distribution indicator (TPDI)-based approach (Section 3.1).The detected trunk positions were taken as an initial trunk detection result.Then a CHM was generated from the normalized point cloud, smoothed for identification of potential treetops, and segmented using a marker-controlled watershed segmentation (MCWS) algorithm (Section 3.2).The leaf-off LiDAR data acquired in February and the leaf-on data acquired in May were separately used to generate the CHM in order to evaluate the effect of data acquisition season.The potential trunk positions were overlaid on the crown segments so that the segments containing more than one trunk position were identified.Within such segments, the false trunk positions were filtered out by analyzing 3D structure features, i.e., the distribution of LiDAR points in the 3D space (Section 3.3).The methods for comparison with the hybrid method are given in Section 3.4.The accuracy assessment method is explained in Section 3.5.

Potential Trunk Detection
The TPDI-based approach was proposed by Deng et al. [24] for individual tree detection and segmentation.The principle of this approach is that the LiDAR points on a tree

Potential Trunk Detection
The TPDI-based approach was proposed by Deng et al. [24] for individual tree detection and segmentation.The principle of this approach is that the LiDAR points on a tree trunk should be evenly distributed along the vertical direction.The steps of the TPDI-based approach are as follows (Figure 4).trunk should be evenly distributed along the vertical direction.The steps of the TPDIbased approach are as follows (Figure 4).Firstly, a layer of LiDAR points within a certain height interval was extracted from the normalized non-ground points.The lower and upper bounds of the layer were defined depending on the crown base heights of the trees in the experimental plots.Because many trees in the plots have a relatively low crown base, the lower bound of 1 m and the upper bound of 4 m were employed.If the upper bound was assigned a larger value, many branch points would be included.The points within the layer were then divided into a 3D grid.Each grid cell is called a voxel.The TPDI was calculated based on the voxels as follows: , = { 0,  ,, = 0 or  ,, = 0 1,  ,, ≠ 0 and  ,, ≠ 0 where i, j, and k are the voxel indices, TPDIi,j represents the TPDI value on the location of the i-th row and j-th column, m is the total number of voxels on the location (i, j), and ni,j,k is the point count in the i,j,k-th voxel.The resulting raster has a TPDI value on each cell in the i-th row and j-th column.The voxel sizes, represented by the horizontal size (dH) and vertical height (dV), were defined based on the LiDAR point density and the maximum tree DBH in the plots.As indicated in the study of Deng et al. [24], if dH is small (e.g., 0.1 m), a trunk may span many cells in the TPDI raster, and the TPDI value is lowered.Additionally, if the trunk points are sparse, a relatively large dV should be adopted.Therefore, the same voxel size (dH = 0.25 m and dV = 0.5 m) used by Deng et al. was adopted in this study.The dH value of 0.25 m is approximately half the maximum DBH (Table 1) and is greater than the DBH of most trees in each plot.
To reduce the effects of the voxel size and the layer height, the TPDI raster was normalized using the following formula: where nTPDIi,j represents the normalized TPDI value for the i,j-th raster cell, TPDImin and TPDImax are the minimum and maximum TPDI values, respectively.Firstly, a layer of LiDAR points within a certain height interval was extracted from the normalized non-ground points.The lower and upper bounds of the layer were defined depending on the crown base heights of the trees in the experimental plots.Because many trees in the plots have a relatively low crown base, the lower bound of 1 m and the upper bound of 4 m were employed.If the upper bound was assigned a larger value, many branch points would be included.The points within the layer were then divided into a 3D grid.Each grid cell is called a voxel.The TPDI was calculated based on the voxels as follows: where i, j, and k are the voxel indices, TPDI i,j represents the TPDI value on the location of the i-th row and j-th column, m is the total number of voxels on the location (i, j), and n i,j,k is the point count in the i,j,k-th voxel.The resulting raster has a TPDI value on each cell in the i-th row and j-th column.The voxel sizes, represented by the horizontal size (d H ) and vertical height (d V ), were defined based on the LiDAR point density and the maximum tree DBH in the plots.As indicated in the study of Deng et al. [24], if d H is small (e.g., 0.1 m), a trunk may span many cells in the TPDI raster, and the TPDI value is lowered.Additionally, if the trunk points are sparse, a relatively large d V should be adopted.Therefore, the same voxel size (d H = 0.25 m and d V = 0.5 m) used by Deng et al. was adopted in this study.The d H value of 0.25 m is approximately half the maximum DBH (Table 1) and is greater than the DBH of most trees in each plot.
To reduce the effects of the voxel size and the layer height, the TPDI raster was normalized using the following formula: where nTPDI i,j represents the normalized TPDI value for the i,j-th raster cell, TPDI min and TPDI max are the minimum and maximum TPDI values, respectively.A local maximum filter in the form of a sliding window with a fixed size (w T ) was applied to the nTPDI raster to identify potential trunk positions.The value of w T was specified based on the minimum distance between tree trunks (D min ) and the maximum DBH (DBH max ) in the plots, following: We adopted a window size smaller than 2 × D min because the local maximum of a trunk may deviate from the trunk center, and thus the distance between the local maxima of two neighboring trunks may be smaller than D min .
Subsequently, a vertical cylinder with a constant radius (r = D min /2) was centered at each potential trunk position and RANSAC-based 3D line fitting was applied to the LiDAR points inside the cylinder.The RANSAC method randomly selected two points to compute the parameters of the 3D line.The points lying within a distance (δ = DBH max /2) were taken as inlier points.This process was repeated for N = 500 iterations and the optimal 3D line fit was determined to be the one with the largest ratio of inlier points among all iterations.The inclination angle of the 3D line was compared with a threshold.If the inclination angle of the 3D line on a trunk position was larger than the threshold, the position was regarded as a false trunk position.Considering the positional error of LiDAR data and the inclined trunks, a loose threshold of 45 • was adopted.Such a threshold value is appropriate for most cases and does not need to be changed when applied to other forests.

Crown Segmentation
To segment individual crowns, a CHM was generated from the normalized point cloud by interpolating the aboveground heights of the highest laser points in each grid cell using the triangulation-based natural neighbor interpolation method.The CHM was then smoothed using a Gaussian filter to reduce the effect of canopy surface irregularities [37].A sliding window was applied to the smoothed CHM to identify local maxima, which were regarded as potential treetops.Several parameters should be defined for the Gaussian filter and the sliding window.The size of the sliding window (w C ) for the treetop detection was defined in the same way as that for the identification of potential trunk positions, i.e., w C = 2 × ⌈(D min − DBH max /2)/CS⌉ + 1, where CS represents the cell size of the CHM.The Gaussian filter was also applied within a window moving over each CHM cell.The width of the moving window was defined as ⌊w C /2⌋, and the standard deviation of the Gaussian filter was 0.25 m.
Subsequently, individual crowns were segmented using an MCWS algorithm similar to the one proposed by Chen et al. [38].The operations of opening and closing by reconstruction were first performed on the CHM image.Regional maxima (connected components with a constant value) were then extracted from the processed CHM and were taken as markers for the watershed algorithm to delineate crown segments.Since the morphological operations can smooth the CHM, the MCWS algorithm was applied to the original CHM rather than the CHM smoothed by the Gaussian filter.We did not take the potential treetops as markers because the crown of a large broadleaf tree may have multiple peaks and over-segmentation tends to occur when using the potential treetops as markers.In contrast, the MCWS algorithm based on morphological operations can avoid the over-segmentation problem to a certain degree [38].The detected potential treetops were used in the following refinement process.

Refinement of Initial Trunk Detection Result
The potential trunk positions and the crown segmentation result were stacked to identify the crown segments containing multiple potential trunk positions.Such segments were processed one after another to filter out false trunk positions.The steps are given as follows.
Step 1. Identification of potentially false trunk positions.
For a crown segment containing multiple potential trunk positions, the LiDAR points within a certain height interval (1-4 m) were extracted and a 2D DBSCAN algorithm was used to cluster the points.The DBSCAN algorithm defines three types of points: core points, border points, and noise points [39].Starting from a random point P, the points in its neighborhood with a radius of ε 1 were searched.If the number of neighboring points is not less than the minimum number of neighboring points (minPts 1 ), P is defined as a core point.A border point connects to a core point but does not have enough neighboring points, while the remaining points are the noise points.The connected core points and border points form a cluster [20].In this step, the distances between points were calculated using the planimetric coordinates of the LiDAR points.The value of ε 1 influences the clustering result, and more clusters are generated when applying a smaller ε 1 .Since the aim is to identify potentially false trunk positions, ε 1 is related to the thickness of a trunk, and it was set to DBH max /2.The value of minPts 1 has a minor influence and was set to 3.After clustering, the potential trunk positions were superimposed on the clusters.If a cluster enclosed multiple trunk positions (Figure 5), these trunk positions were possibly associated with bushes or low branches and leaves and should be further evaluated in the following steps.
Forests 2024, 15, x FOR PEER REVIEW 10 of 22 Step 1. Identification of potentially false trunk positions.For a crown segment containing multiple potential trunk positions, the LiDAR points within a certain height interval (1-4 m) were extracted and a 2D DBSCAN algorithm was used to cluster the points.The DBSCAN algorithm defines three types of points: core points, border points, and noise points [39].Starting from a random point P, the points in its neighborhood with a radius of ε1 were searched.If the number of neighboring points is not less than the minimum number of neighboring points (minPts1), P is defined as a core point.A border point connects to a core point but does not have enough neighboring points, while the remaining points are the noise points.The connected core points and border points form a cluster [20].In this step, the distances between points were calculated using the planimetric coordinates of the LiDAR points.The value of ε1 influences the clustering result, and more clusters are generated when applying a smaller ε1.Since the aim is to identify potentially false trunk positions, ε1 is related to the thickness of a trunk, and it was set to DBHmax/2.The value of minPts1 has a minor influence and was set to 3.After clustering, the potential trunk positions were superimposed on the clusters.If a cluster enclosed multiple trunk positions (Figure 5), these trunk positions were possibly associated with bushes or low branches and leaves and should be further evaluated in the following steps.Step 2. Analysis of 3D structure on potentially false trunk positions.For a potentially false trunk position, the treetops extracted in Section 3.2 were searched around the trunk position with a small radius.The trunk position was considered correct with a high probability if at least one treetop was found near the trunk position, and the subsequent validation was no longer carried out.Considering the 10 cm horizontal accuracy of the LiDAR data and the trunk inclination, the radius was set to DBHmax + 0.2.If no treetop was found, a cylinder with a radius of  min − DBH max /2 was centered on the trunk position, within which the 3D structure of the LiDAR points was analyzed as follows: (1) Remove crown points based on the analysis of a vertical profile.The LiDAR points were binned with a vertical resolution of 0.1 m, generating a height histogram representing the vertical profile.Each bin in the histogram has a value of the percentage of points, i.e., ni/N, where ni is the point count within this bin and N is the total point number.The largest peak in the histogram at the height, Hp, corresponds to a value represented by maxP.The height of the highest LiDAR point is denoted by maxH.If Hp was above threequarters of maxH, the approximate crown base height (CBH) was estimated as the height corresponding to the highest bin with a percentage of points below 0.1 × maxP (Figure 6a).Step 2. Analysis of 3D structure on potentially false trunk positions.For a potentially false trunk position, the treetops extracted in Section 3.2 were searched around the trunk position with a small radius.The trunk position was considered correct with a high probability if at least one treetop was found near the trunk position, and the subsequent validation was no longer carried out.Considering the 10 cm horizontal accuracy of the LiDAR data and the trunk inclination, the radius was set to DBH max + 0.2.If no treetop was found, a cylinder with a radius of D min − DBH max /2 was centered on the trunk position, within which the 3D structure of the LiDAR points was analyzed as follows: (1) Remove crown points based on the analysis of a vertical profile.The LiDAR points were binned with a vertical resolution of 0.1 m, generating a height histogram representing the vertical profile.Each bin in the histogram has a value of the percentage of points, i.e., n i /N, where n i is the point count within this bin and N is the total point number.The largest peak in the histogram at the height, H p , corresponds to a value represented by maxP.The height of the highest LiDAR point is denoted by maxH.If H p was above threequarters of maxH, the approximate crown base height (CBH) was estimated as the height corresponding to the highest bin with a percentage of points below 0.1 × maxP (Figure 6a).If the estimated CBH was lower than half the maxH, which meant a large percentage of LiDAR points were below the crown, the CBH was re-estimated as H p − maxH − H p .The LiDAR points above the CBH were regarded as crown points and removed.If the height corresponding to the largest peak was below three-quarters of maxH, there could be either a low tree covered by the upper canopy or large branches and leaves at the potential trunk position (Figure 6b).The highest void space with a span larger than 0.5 m was searched and the LiDAR points above the void space were removed.Since the aim is to reduce the effect of the crown points in the following steps, we only need to remove most crown points and a precise CBH estimation is unnecessary.If the estimated CBH was lower than half the maxH, which meant a large percentage of LiDAR points were below the crown, the CBH was re-estimated as  p − ( −  p ).The LiDAR points above the CBH were regarded as crown points and removed.If the height corresponding to the largest peak was below three-quarters of maxH, there could be either a low tree covered by the upper canopy or large branches and leaves at the potential trunk position (Figure 6b).The highest void space with a span larger than 0.5 m was searched and the LiDAR points above the void space were removed.Since the aim is to reduce the effect of the crown points in the following steps, we only need to remove most crown points and a precise CBH estimation is unnecessary.(2) Remove the points of large branches.Trunk, branch, and foliage points differ by the relative position of their neighboring points [40].To distinguish the true trunk positions from the false ones caused by branches and leaves, the LiDAR points remaining after the crown point removal were clustered using the DBSCAN algorithm.First, a Principal Component Analysis (PCA) was performed to calculate the maximum extension direction corresponding to the first component for the LiDAR points.All the LiDAR points were then projected to the vertical plane containing the vector of the maximum extension direction.Next, the projected points were clustered using a 2D DBSCAN algorithm with ε2 (2) Remove the points of large branches.Trunk, branch, and foliage points differ by the relative position of their neighboring points [40].To distinguish the true trunk positions from the false ones caused by branches and leaves, the LiDAR points remaining after the crown point removal were clustered using the DBSCAN algorithm.First, a Principal Component Analysis (PCA) was performed to calculate the maximum extension direction corresponding to the first component for the LiDAR points.All the LiDAR points were then projected to the vertical plane containing the vector of the maximum extension direction.Next, the projected points were clustered using a 2D DBSCAN algorithm with ε 2 = 0.1 m and minPts 2 = 3, where ε 2 and minPts 2 are the neighborhood radius and minimum number of neighboring points.A relatively small ε 2 was employed because using a large ε 2 could result in large-size clusters containing multiple branches and even the trunk.A 3D DBSCAN with ε 3 = DBH max /2 and minPts 3 = 3 was applied to each cluster to remove isolated points and small clusters with a size ≤ minPts 3 .We did not directly apply a 3D DBSCAN to the LiDAR point cloud because the points of branches were usually incomplete due to the occlusion of the canopy.We found that a 3D DBSCAN clustering may break a large branch into small clusters, while a 2D DBSCAN in the projection plane may generate more complete branch clusters when the same parameters were adopted.After DBSCAN clustering, the maximum extension direction of each cluster was computed through PCA.The angle between the maximum extension direction and the vertical direction was then calculated.If the angle exceeded the loose threshold of 45 • , the cluster was regarded as a large branch and removed.Figure 7 displays two examples of DBSCAN clustering and removal of the points of large branches.
= 0.1 m and minPts2 = 3, where ε2 and minPts2 are the neighborhood radius and minimum number of neighboring points.A relatively small ε2 was employed because using a large ε2 could result in large-size clusters containing multiple branches and even the trunk.A 3D DBSCAN with ε3 = DBHmax/2 and minPts3 = 3 was applied to each cluster to remove isolated points and small clusters with a size ≤ minPts3.We did not directly apply a 3D DBSCAN to the LiDAR point cloud because the points of branches were usually incomplete due to the occlusion of the canopy.We found that a 3D DBSCAN clustering may break a large branch into small clusters, while a 2D DBSCAN in the projection plane may generate more complete branch clusters when the same parameters were adopted.After DBSCAN clustering, the maximum extension direction of each cluster was computed through PCA.The angle between the maximum extension direction and the vertical direction was then calculated.If the angle exceeded the loose threshold of 45°, the cluster was regarded as a large branch and removed.Figure 7 displays two examples of DBSCAN clustering and removal of the points of large branches.(3) Distinguish between true and false trunk positions based on RANSAC-based 3D line fitting.After removing the points of large branches, we assume that the remaining points on a true trunk position should be mostly trunk points (Figure 7b).If a 3D line is (3) Distinguish between true and false trunk positions based on RANSAC-based 3D line fitting.After removing the points of large branches, we assume that the remaining points on a true trunk position should be mostly trunk points (Figure 7b).If a 3D line is fitted to the points, it should have a small inclination angle.In contrast, the remaining points on a false trunk position should be distributed disorderly or have a large inclination angle (Figure 7a).Considering that some trees have inclined trunks, we divided the point cloud into two layers and applied the RANSAC-based 3D line fitting separately in the two layers.The first layer was the same as the layer for the detection of potential trunk positions (Section 3.1).The second layer spanned from the upper bound of the first layer to the height of the highest point.After fitting 3D lines in the two layers, the horizontal distance between the upper endpoint of the first line and the lower endpoint of the second line was calculated.If the distance did not exceed a threshold T D and at least one 3D line had an inclination angle <45 • , the trunk position was regarded as a true one.The definition of T D is associated with the maximum DBH (DBH max ).Considering the 10 cm horizontal accuracy of the UAV-LiDAR data, T D was specified as DBH max + 0.2.

Methods for Comparison
We selected two individual tree detection and segmentation methods for comparison.The first method detects treetops based on a variable-size window which is adaptively adjusted depending on the relationship between tree height and crown size [41].We performed the method in Fusion v4.41 software.A mean filter was applied to smooth the CHM before detecting treetops.The main parameters include the CHM resolution and the size of the moving window for the mean filter.To select the optimal parameters, we separately applied the moving window with sizes of 3, 5, and 7 cells to the CHM with a 0.25 m resolution and the moving window with sizes of 11, 13, 15, 17, and 19 cells to the CHM with a 0.1 m resolution.In addition, the method was applied to the leaf-off and leaf-on LiDAR data, respectively, to evaluate the effect of data acquisition season.
The second method was a region-growing method (Li2012) working at the point cloud level [14].The method segments individual trees in sequence from the point cloud by taking advantage of the relative spacing between trees.We performed the Li2012 method from the R package lidR [42].The parameters that need to be specified include two distance thresholds (dt1 and dt2), a searching radius (R), and the minimum height (h min ) of detected trees.The distance thresholds are associated with both the minimum distance between tree trunks (D min ) and the crown diameter.The dt1 and dt2 values of 1, 1.5, 2, 2.5, and 3 were separately applied to select the optimal value.The value of R is also related to D min and was set to 1.5 m.The threshold h min was set to 3.5 m.The leaf-off and leaf-on LiDAR data were separately utilized to segment individual trees.The planimetric coordinates of the highest point within each tree segment were taken as the tree position and used for accuracy calculation.

Accuracy Assessment
The individual tree detection accuracy was assessed by comparing the final tree trunk positions with the reference data.Within a searching distance from a reference tree position, the nearest detected trunk was searched and regarded as a true positive (TP).The reference and the detected trunk position were labeled as matching.If no detected trunk position can be found within the searching distance, the reference tree position is regarded as a false negative (FN).The unmatched detected trunk positions were regarded as false positives (FPs).Recall (R), precision (P), and F-score (F) were calculated for each plot using the following equations [24]: R = N TP /(N TP + N FN ) (4) where N TP , N FN , and N FP are the number of TPs, FNs, and FPs, respectively.The searching distance was defined based on the maximum DBH and the minimum distance between trees in the plots.Since some detected trunk positions deviated from the reference positions due to the inclined trunks, a searching distance of 1 m was utilized to calculate the accuracy of our hybrid method.Because some trees were located on the plot boundaries and only a part of the trunk and crown fell within the boundaries, the accuracy calculation only considered the detected and reference trees within a smaller extent that moved inward 0.5 m from each boundary.The 1.0 m searching distance is a relatively strict threshold, as our hybrid method directly extracted trunk positions, rather than extracting treetops to represent tree positions.Due to the deviation between a treetop and a trunk position when the trunk is inclined, a larger searching distance (1.5 m) was used to calculate the accuracy when comparing our hybrid method with the treetop detection and the tree segmentation methods.

Detection of Individual Trees Using Both Leaf-Off and Leaf-On LiDAR Data
Since the tree trunk detection requires adequate LiDAR points of the tree trunks, the detection of potential trunk positions and the refinement of initial trunk detection results were performed on the leaf-off LiDAR data acquired in February.The leaf-on LiDAR data acquired in May was used to identify potential tree tops and segment individual crowns.In our hybrid method, most parameters and thresholds were specified based on the minimum distance between tree trunks (D min ) and the maximum DBH (DBH max ).According to the field-measured data (Section 2.3), D min was set to 1.5 m and DBH max was 0.5 m.As the CHM resolution can influence the crown segmentation, the MCWS was performed on the CHMs with a resolution of 0.25 m and 0.1 m, respectively.The final individual tree detection results in the three plots are shown in Table 2.
Table 2. Accuracies of individual tree detection using the hybrid method based on both leaf-off and leaf-on LiDAR data.The searching distance for accuracy assessment is 1.0 m.

Plot
No As indicated in Table 2, the highest and the lowest F-scores were derived in plots P1 and P3, respectively.Compared to the other two plots, plot P3 contains more trees with low heights, low crown bases, and multiple trunks growing together, such as osmanthus.The trunks of such trees had barely LiDAR points and were undetectable in the point cloud.On the other hand, their low crowns which nearly touched the ground resulted in many false trunk positions.In plot P2, the deciduous trees are mainly concentrated along the eastern boundary of the plot and the rest of the trees are almost evergreen trees.Due to the dense canopy, some trees had sparse LiDAR points on their trunks and, hence, were difficult to detect.
Table 2 also shows the effect of the CHM resolution on the crown segmentation and the final individual tree detection result.When the CHM resolution of 0.25 m was adopted, the number of crown segments in all three plots was less than the number of reference trees, indicating under-segmentation.As the CHM resolution varied from 0.25 to 0.1 m, many more crown segments were generated and over-segmentation occurred.Despite the great change in the number of crown segments, the variation in the F-score value was small.The largest variation in F-score (0.025) can be found in plot P1.In this plot, the number of crown segments grew from 69 to 230, whereas the number of detected individual trees increased by only 7. In all three plots, the F-scores derived by using the 0.25 m resolution CHM were higher than those derived by using the 0.1 m resolution CHM.This is because in the case of over-segmentation, many crown segments contained none or only one potential trunk position and, hence, the false trunk positions were not validated by the refinement procedure.

Detection of Individual Trees Using Leaf-Off LiDAR Data Alone
To analyze the effect of the season for data acquisition on crown segmentation and individual tree detection, the leaf-off LiDAR data were applied to both trunk detection and crown segmentation.The final individual tree detection results derived by using leaf-off LiDAR data alone are shown in Table 3.
Table 3. Accuracies of individual tree detection using the hybrid method based on leaf-off LiDAR data alone.The searching distance for accuracy assessment is 1.0 m.

Plot
No Similar to Table 2, Table 3 also shows a significant increase in crown segment number and a decrease in F-score in all three plots when the CHM resolution varied from 0.25 to 0.1 m, indicating that false trunk positions tended to be retained in the final result in the case of over-segmentation.The largest variation in F-score (0.030) can be found in plot P3.A total of 19 trunk positions, including 4 true positions and 15 false positions, were removed when using a 0.25 m resolution CHM, but retained in the final result when using a 0.1 m resolution CHM on plot P3.

Comparison with Existing Methods
A treetop detection method implemented in Fusion software and a tree segmentation method (Li2012) were selected for comparison with our hybrid method.The two methods were applied to the leaf-off and leaf-on LiDAR data, respectively.By trial-and-error of different parameter values (Section 3.4) in each plot using either leaf-off or leaf-on data, the results with the highest F-scores are shown in Table 4.The optimal LiDAR data were the leaf-off data for both methods in all three plots.The results of our hybrid method shown in the table were derived by using the leaf-off LiDAR data for trunk detection and leaf-on data for crown segmentation.The CHM resolution of 0.25 m was adopted.
As shown in Table 4, both the treetop detection method and the Li2012 method derived much lower accuracies than the hybrid method, especially in plots P1 and P2.The low accuracy is due to the overlapping canopy and the multi-peaks of large crowns of broadleaf trees.The overlapping canopy tends to result in under-segmentation, while the multi-peaks of large-size crowns may result in over-segmentation.The low accuracy may also be due to the large deviations between treetops and trunk positions.This means that using treetops to represent tree positions may result in large errors in tree positions.
We also calculated the accuracies of the initial trunk detection result (Section 3.1) and the tree detection result without crown segmentation, i.e., all potential trunk positions were validated using the steps given in Section 3.3.As indicated in Table 4, the initial trunk detection results had high recall and low precision values.Both true and false trunk positions were filtered out by the refinement processing, but more false trunk positions were filtered out than true trunk positions.Therefore, by carrying out the refinement steps, the F-score was increased by 0.062, 0.066, and 0.115 in plots P1, P2, and P3, respectively.In contrast, the tree detection without crown segmentation derived low recall and high precision values, indicating that many true trunk positions were removed after the refine-ment.This also implies that the 3D structures of some correctly and incorrectly detected trunks were indistinguishable in the point clouds, possibly due to the sparse trunk points, abundant lateral branches, and inclined trunks.Nevertheless, integrating trunk detection with crown segmentation can effectively improve individual tree detection accuracy.

Challenges in This Study
The aim of this study is to detect the trunk positions of individual trees in broadleaf forests.The study was conducted in an urban area covered by patches of broadleaf forests.The first challenge is the relatively sparse trunk points compared to the terrestrial LiDAR data, meaning that the methods developed for the terrestrial LiDAR data may be unsuitable for this task.The second challenge involves abundant lateral branches and low crown bases of some trees, which make it hard to identify individual trunks from the point cloud.Figure 8 displays partial LiDAR points within a certain height interval together with the reference trees in plot P1.Within the 1-2 m height interval, some overstory and understory trees have no corresponding LiDAR points.Within the 2-4 m height interval, almost all the trees have corresponding LiDAR points but some large clusters of LiDAR points enclose multiple trees.The commonly used clustering methods, such as DBSCAN clustering [21] or clustering based on horizontal distances [26], may be ineffective in such a case.The tree trunks enclosed in a large-size cluster as shown in Figure 8 could not be separated from each other using the clustering methods.
In our hybrid method, the true and false trunk positions were distinguished through analysis of the 3D structure of the point cloud around a trunk position.Some similar studies also proposed methods for 3D structure analysis.For example, Huo et al. [25] distinguished between true and false treetops based on the trees' symmetrical structure represented by symmetry curves generated from the point cloud.However, in our study area, some trees have unsymmetrical lateral branches or inclined trunks.Furthermore, due to the relatively low point density underneath the canopy, the difference in 3D structure between the trunks and the large branches or low crowns is not distinct.In this case, it is difficult to distinguish between true and false trunk positions based on trees' symmetrical structure.
tree trunks enclosed in a large-size cluster as shown in Figure 8 could not be separated from each other using the clustering methods.In our hybrid method, the true and false trunk positions were distinguished through analysis of the 3D structure of the point cloud around a trunk position.Some similar studies also proposed methods for 3D structure analysis.For example, Huo et al. [25] distinguished between true and false treetops based on the trees' symmetrical structure represented by symmetry curves generated from the point cloud.However, in our study area, some trees have unsymmetrical lateral branches or inclined trunks.Furthermore, due to the relatively low point density underneath the canopy, the difference in 3D structure between the trunks and the large branches or low crowns is not distinct.In this case, it is difficult to distinguish between true and false trunk positions based on trees' symmetrical structure.

Performance of Our Hybrid Method
In the proposed hybrid method, we identified potential trunk positions based on a TPDI indicator and RANSAC-based 3D line fitting, rather than using clustering methods to directly segment individual trunks.The treetops identified from the CHM and the crown segments derived by MCWS were then used to find potentially false trunk positions.The 3D structure of the point cloud at each potentially false trunk position was analyzed to distinguish between true and false trunk positions.It is based on the observation that the LiDAR points of trunks are distributed along the vertical or near vertical direction, while the distribution of branch and foliage points is either random or approaches the horizontal direction.Therefore, the crown points were first removed based on an estimated CBH and the points of large branches were then removed based on their maximum extension directions.On a true trunk position, the remaining points should be trunk points that were vertically distributed.Our hybrid method achieved relatively high accuracies in recall, precision, and F-score, especially compared to the treetop detection and the individual tree segmentation methods (Table 4).
We did not validate each potential trunk position but found potentially false trunk positions based on crown segments for validation.According to the comparison with the refinement result without crown segmentation (Table 4), adding the constraint of crown segments greatly improved the F-scores.The improvement in F-score ranged from 0.077 to 0.133.This is because the 3D structures of some correctly and incorrectly detected trunks were indistinguishable in the point clouds.Taking the crown segments as a

Performance of Our Hybrid Method
In the proposed hybrid method, we identified potential trunk positions based on a TPDI indicator and RANSAC-based 3D line fitting, rather than using clustering methods to directly segment individual trunks.The treetops identified from the CHM and the crown segments derived by MCWS were then used to find potentially false trunk positions.The 3D structure of the point cloud at each potentially false trunk position was analyzed to distinguish between true and false trunk positions.It is based on the observation that the LiDAR points of trunks are distributed along the vertical or near vertical direction, while the distribution of branch and foliage points is either random or approaches the horizontal direction.Therefore, the crown points were first removed based on an estimated CBH and the points of large branches were then removed based on their maximum extension directions.On a true trunk position, the remaining points should be trunk points that were vertically distributed.Our hybrid method achieved relatively high accuracies in recall, precision, and F-score, especially compared to the treetop detection and the individual tree segmentation methods (Table 4).
We did not validate each potential trunk position but found potentially false trunk positions based on crown segments for validation.According to the comparison with the refinement result without crown segmentation (Table 4), adding the constraint of crown segments greatly improved the F-scores.The improvement in F-score ranged from 0.077 to 0.133.This is because the 3D structures of some correctly and incorrectly detected trunks were indistinguishable in the point clouds.Taking the crown segments as a constraint may reduce the number of true trunk positions being filtered out.On the other hand, adding such a constraint also reduced the computation time.
The crown segmentation result had an influence on the individual tree detection result.In Section 4, we compared the tree detection results derived by using the CHMs of 0.25 and 0.1 m resolution.Using the CHM with a 0.25 m resolution resulted in undersegmentation while using the CHM with a 0.1 m resolution resulted in over-segmentation.In the case of over-segmentation, false trunk positions tended to be retained in the final result, represented by the lower Precision values.In plots P1 and P2, the recall value was almost invariant when the CHM resolution changed, indicating that the true trunk positions were less affected by the CHM resolution.The higher F-scores derived by using the 0.25 m resolution CHM in all three plots indicate that under-segmentation is preferred when applying our hybrid method.

Effect of Data Acquisition Season
In this study, the leaf-off LiDAR data were used for trunk position identification and the subsequent 3D structure analysis.The leaf-on data were not used for trunk position identification due to severe obstruction by the dense canopy during the leaf-on season.The leaf-off and leaf-on data were separately applied to the crown segmentation to evaluate the effect of the data acquisition season.By comparing the results in Tables 2 and 3, it can be seen that the difference in F-score between the results obtained using different datasets (leaf-off and leaf-on) for the crown segmentation was small (ranging from 0.005 to 0.029), indicating a small influence of the data acquisition season on the individual tree detection when using our hybrid method.
The effect of the data acquisition season was also evaluated for the treetop detection and the individual tree segmentation methods.Both methods derived higher F-scores in all three plots when using the leaf-off data.The reason for the higher accuracies derived from the leaf-off data may be that the treetops and the gap between crowns are more distinct in the leaf-off data.
The three plots used for the evaluation of our hybrid method have different proportions of evergreen and deciduous trees.Plots P1 and P2 are dominated by evergreen trees, while more than half of the trees in plot P3 are deciduous trees.Theoretically, a forest with a higher proportion of deciduous trees should have more laser pulses penetrating the canopy layer, and the individual tree detection accuracy should be higher.Since plot P1 has a higher proportion of deciduous trees (43.3%) than plot P2 (34.8%), higher recall and F-score values were derived in plot P1.An exception is plot P3, within which the lowest F-score was obtained.This is because the deciduous trees in plot P3 are concentrated in the southwest part of the plot and there are many undetectable low trees with no trunk points.
Nevertheless, due to the requirement of the leaf-off LiDAR data for the identification of potential trunk positions and analysis of 3D structure, a single leaf-off datum is adequate in practice to reduce the data collection efforts in multiple seasons.

Limitations
The hybrid method proposed in this work needs to specify some parameters.The main parameters include the height interval and the voxel size for the calculation of TPDI, the size of the sliding window for the identification of potential trunk positions, and the neighborhood radius of the 2D or 3D DBSCAN algorithm in the refinement procedure.The sensitivity of the parameter specification for the identification of potential trunk positions can be referred to in the study of Deng et al. [24].In this work, we defined the parameters based on the prior knowledge of the forests in the study area, namely, the minimum distance between tree trunks (D min ) and the maximum DBH (DBH max ).In practice, these two values are difficult to acquire and can only be estimated.An overestimated D min could result in a larger sliding window and thus fewer trees could be detected.
In the hybrid method, the crown segments that did not contain trunk positions were not processed further, even though some of them corresponded to undetected trees.This is because the crown segments corresponding to undetected trees usually had very sparse trunk laser points (<5 points), resulting in the difficulty of determining whether such crown segments corresponded to undetected trees or were caused by over-segmentation.The method needs to be improved to validate the crown segments without trunk positions based on additional evidence, such as the morphology of crown segments.
We evaluated the hybrid method on three sample plots and the forest type is the urban broadleaf forest containing both evergreen and deciduous trees.Although the trees are unevenly distributed in the plots, the understory vegetation is relatively sparse, which is unlike the natural subtropical forests with dense and high understory vegetation.Moreover, the number and sizes of the plots in this study are relatively small.This is because the target is subtropical broadleaf forests with varying proportions of evergreen trees and the objectives include evaluating the influence of data acquisition season.There is a lack of suitable UAV-LiDAR data collected in multiple seasons over a large area.Further research should evaluate the transferability of the proposed hybrid method in natural forests over a large area.Although multiple-season data were used in this study, the comparison results indicate that using single-season (leaf-off) data can produce promising results.Therefore, when applying our hybrid method to a larger scale, single-season data can be used and the main problems are the computation cost and acquisition of the parameters, such as the minimum distance between trees and the maximum DBH.One possible solution is to use LiDAR data to generate estimated values for these parameters, and then adjust them based on forest type and age.
Our hybrid method requires adequate LiDAR points along the tree trunks and, hence, the UAV-LiDAR data with a high point density was utilized in this study.In addition, the leaf-off LiDAR data were used for trunk position identification as more laser pulses penetrated the canopy layer in the leaf-off season.Although our method has taken into account the sparse trunk laser points, the tree trunks with too few points (<5 points) could not be detected in all three plots.Further research should evaluate the applicability of our hybrid method to airborne LiDAR data with lower point densities.This depends on the penetration rate of laser pulses, which is affected by forest type and structure.
In our study area, there are some trees with a relatively low tree height and low crown base height (<2 m).These trees were almost undetectable by using the hybrid method due to the lack of trunk points.The 3D structures of such trees should be further studied to improve the method.In addition, we did not segment individual trees but only detected trunk positions.The first reason is the trunk position is an important parameter in precise tree inventory.The second reason is the reference crown boundaries are difficult to obtain due to the overlapping canopy, and thus the segmentation accuracy could not be assessed.The commonly used segmentation methods, e.g., region growing, can be adapted to segment individual trees by taking the detected trunk positions as seeds.

Conclusions
In this study, we proposed a hybrid method for detecting individual trees in broadleaf forests.Potential trunk positions were first detected from the normalized TPDI raster and were then validated based on the crown segmentation result and the 3D structures of trunks and branches in the point cloud.The method was evaluated on three experimental plots in urban broadleaf forests with varying proportions of evergreen trees and tree species composition.The leaf-off LiDAR data were used for the identification of potential trunk positions, while the leaf-off and leaf-on LiDAR data were separately applied to the crown segmentation, for the analysis of the influence of data acquisition season.In addition, the 0.25 m and 0.1 m resolution CHMs were separately used for the crown segmentation to analyze the influence of crown segmentation on the individual tree detection result.The comparison results indicate that the tree detection results derived by using the 0.25 m resolution CHM had higher accuracies, and the influence of data acquisition season was small.The highest F-score in each plot ranged from 0.695 to 0.820 when using both leaf-off and leaf-on LiDAR data, while the range was from 0.676 to 0.830 when using the leaf-off data alone.The variation in the F-score was related to the species composition and the spatial pattern of the deciduous and evergreen trees in each plot.The performance of the hybrid method was compared with a treetop detection and a point cloud-based individual tree segmentation method.Much lower F-scores were obtained by both methods.In addition, both methods derived higher F-scores when using the leaf-off LiDAR data than using the leaf-on data.
One limitation of the proposed hybrid method is that the specification of parameters depends on the prior knowledge of the forest, mainly the minimum distance between tree trunks and the maximum DBH.The second limitation is that the hybrid method was evaluated on three relatively small plots in urban broadleaf forests, which differ from the natural subtropical forest in terms of forest structure.Further research should be conducted to evaluate the transferability of the proposed hybrid method in natural forests over a large area.In addition, this study used UAV-LiDAR data with a high point density to ensure the

Figure 1 .
Figure 1.Experimental plots in the study area.(a,b) show the locations of the three plots, P1, P2, and P3; (c-e) display the UAV-RGB images of the three plots, respectively.

Figure 1 .
Figure 1.Experimental plots in the study area.(a,b) show the locations of the three plots, P1, P2, and P3; (c-e) display the UAV-RGB images of the three plots, respectively.

Figure 2 .
Figure 2. Elevation-rendering point clouds on three experimental plots acquired in February (left column) and May (right column).(a,b) Point clouds in plot P1; (c,d) point clouds in plot P2; (e,f) point clouds in plot P3.

Figure 2 .
Figure 2. Elevation-rendering point clouds on three experimental plots acquired in February (left column) and May (right column).(a,b) Point clouds in plot P1; (c,d) point clouds in plot P2; (e,f) point clouds in plot P3.

Figure 3 .
Figure 3. Workflow of the proposed hybrid individual tree detection method.nTPDI: normalized trunk point distribution indicator; RANSAC: random sample consistency.

Figure 3 .
Figure 3. Workflow of the proposed hybrid individual tree detection method.nTPDI: normalized trunk point distribution indicator; RANSAC: random sample consistency.

Figure 4 .
Figure 4. Steps of the TPDI-based trunk detection approach.nTPDI: normalized trunk point distribution indicator; RANSAC: random sample consistency.The red and blue colors of the normalized point cloud correspond to high and low points above ground.

Figure 4 .
Figure 4. Steps of the TPDI-based trunk detection approach.nTPDI: normalized trunk point distribution indicator; RANSAC: random sample consistency.The red and blue colors of the normalized point cloud correspond to high and low points above ground.

Figure 5 .
Figure 5.The 2D DBSCAN clustering of the LiDAR points between 1 and 4 m in height in a crown segment.The potential tree trunk positions are superimposed on the clustering result.

Figure 5 .
Figure 5.The DBSCAN clustering of the LiDAR points between 1 and 4 m in height in a crown segment.The potential tree trunk positions are superimposed on the clustering result.

Figure 6 .
Figure 6.Removal of crown points when the largest peak in the vertical profile is (a) above or (b) below three-quarters of the maximum height value on a potential trunk position.

Figure 6 .
Figure 6.Removal of crown points when the largest peak in the vertical profile is (a) above or (b) below three-quarters of the maximum height value on a potential trunk position.

Figure 7 .
Figure 7. Two examples of DBSCAN clustering, removal of the points of large branches, and RAN-SAC-based 3D line fitting for determination of (a) false and (b) true trunk positions.The colors in the middle column represent different point clusters.

Figure 7 .
Figure 7. Two examples of DBSCAN clustering, removal of the points of large branches, and RANSAC-based 3D line fitting for determination of (a) false and (b) true trunk positions.The colors in the middle column represent different point clusters.

Figure 8 .
Figure 8. Partial LiDAR points within (a) 1-2 m, and (b) 2-4 m height intervals in plot P1.The reference tree positions are superimposed on the point cloud.

Figure 8 .
Figure 8. Partial LiDAR points within (a) 1-2 m, and (b) 2-4 m height intervals in plot P1.The reference tree positions are superimposed on the point cloud.

Table 1 .
Detailed information on experimental plots.

Table 4 .
Accuracies of individual tree detection using different methods.The searching distance for accuracy assessment is 1.5 m.TPDI-based tree detection refers to the initial tree detection based on TPDI and RANSAC 3D line fitting without the subsequent refinement procedure.** Refinement without crown segmentation refers to validating each trunk position without the constraint of crown segmentation in our method. *