Automatic Detection of Single Trees in Airborne Laser Scanning Data through Gradient Orientation Clustering

Currently, existing methods for single-tree detection based on airborne laser scanning (ALS) data usually require some thresholds and parameters to be set manually. Manually setting threshold or parameters is laborious and time-consuming, and for dense forests, the high commission and omission rate make most existing single-tree detection techniques inefficient. As a solution to these problems, this paper proposed an automatic single-tree detection method in ALS data through gradient orientation clustering (GOC). In this method, the rasterized Canopy Height Model (CHM) was derived from ALS data using surface interpolation. Then, potential trees were assumed as approximate conical shapes and extracted based on the GOC. Finally, trees were identified from the potential trees based on the compactness of the crown shape. This method used the gradient orientation information of rasterized CHM, thus increasing the generalization of single-tree detection method. In order to verify the validity and practicability of the proposed method, twelve 1256 m2 circular study plots of different forest types were selected from the benchmark dataset (NEWFOR), and the results from nine different methods were presented and compared for these study plots. Among nine methods, the proposed method had the highest root mean square matching score (RMS_M = 43). Moreover, the proposed method had excellent detection (M > 47) in both single-layer coniferous and single-layered mixed stands.


Introduction
Remote sensing (RS) is a synthesized, efficient, and rapid means of obtaining earth surface information, which can help forestry management departments achieve high-precision forest resource surveys, timber production estimation, and other applications [1][2][3].Light detection and ranging (LiDAR) is an RS method that uses light in the form of a pulsed laser to measure ranges (i.e., variable distances) to the Earth.These light pulses, combined with other data recorded by the airborne system, generate precise, three-dimensional information about the shape of the Earth and its surface characteristics.Airborne laser scanning (ALS) is a widely used LiDAR technology in forestry, which has a unique advantage in 3D terrain and forest height detection.It can quickly and accurately obtain a digital surface model (DSM) and forest height information for forested land [4,5].According to the statistics, single-tree detection algorithms that use active RS data, such as ALS, have grown rapidly in recent years, accounting for 80% of the total increase in RS area, at the same time, those algorithms using passive data or a fusion of passive and active data accounts for 8% and 12% respectively [6].To date, there are many methods for single-tree detection based on ALS data acquired for forest mapping, precision forestry management, and forest timber evaluation [7][8][9][10][11][12].
Currently, some researchers have proposed automatic or semi-automatic single-tree detection methods.According to the data processing, the single-tree detection method based on ALS data can be divided into three categories: (1) raster-based methods; (2) point cloud-based methods; and (3) methods combing raster, point cloud, and a prior information [13].Our research focuses on the raster-based method which is based on rasterized canopy height model (CHM).This type of method needs to rasterize the point cloud data and detect single trees in the rasterized data.In this method, height normalised ALS points are classified into ground points and above-ground points.Then, the points are rasterised to produce a digital terrain model (DTM, ground points only) and digital surface model (DSM, maximum of all points) by following the generation of a triangulated irregular network.As a result, the CHM was built based on the DTM and the DSM.For example, Hyyppa et al. [14] firstly used high-pulse-rate laser scanners to extract the information at tree level and accurate 3D tree height maps, which is used for forest inventories at stand level.This method calculated the stem volume and basal area for each tree by using the tree height and crown boundary information obtained with the laser scanner for each individual tree.The precision obtained by the new method is better than that in conventional standwise forest inventories.Persson et al. [15] generated digital canopy model with different scales to estimate the height and crown diameter of identified trees by using a parabolic surface.This study demonstrates the usefulness of airborne laser scanning data for detecting and measuring individual trees.Popescu et al. [16] proposed a Local Maximum (LM) technique with a variable window size, which is demonstrated more appropriate than the static LM filter or the grid approach that arbitrarily selects only the highest values in a grid with a user defined size.But, this study was based on rather limited samples with no stratification based on tree species or ages.The LiDAR data were not fully exploited by neglecting secondary returns from within the tree canopy.Chen et al. [17] proposed a marker-controlled watershed segmentation for LiDAR CHM.This method identifies the treetops using variable window sizes to traverse the CHM and then uses these treetops as markers to avoid over segmentation of the watershed algorithm [18,19].This method is effective for isolating trees, and Chen's software TIFFS (toolbox for LiDAR data filtering and forest studies) is well-established in forest remote sensing.The disadvantage of this algorithm is that the single-tree detection results depend on the seed point selection, and over-segmentation often occurs.Dalponte et al. [20] proposed a method to exploit the rasterized CHM and ALS point cloud with normalized heights.This method uses a moving window approach to select seed points if they are the highest point inside the window and higher than 2.5 m.Then, it uses the region growing algorithm to detect the single trees.Mohan et al. [21] proposed an individual tree detection system for open canopy mixed conifer forests.This system evaluates the applicability of low consumer grade cameras attached to unmanned aerial vehicles (UAVs) and a structure-from-motion (SfM) algorithm for automatic individual tree detection using a local-maxima-based algorithm on UAV-derived CHM.Zhao et al. [22] presented a method of region-based hierarchical cross-section analysis for individual tree crown delineation using ALS data.This method uses the features of a mountain-like topographic surface and adopts the vertical structure of crowns to detect treetops and delineate crown boundaries.However, it is difficult to detect and delineate some flat crowns without a dominant protrusion on CHM.Kaartinen et al. [23] conducts lots of experiment to show that the extraction method is the main factor impacting achieved accuracy.However, sometimes the success of tree detection is dependent on tree density and clustering [24].
In this paper, we provide a new idea for single tree detection, clustering by gradient orientation information in CHM, which is named gradient orientation clustering (GOC).We group together the raster cells in CHM whose gradient orientation point to the same one as a potential tree, then use a density threshold that describes the compactness of the shape to select from potential trees as the results.Finally, we compared different model parameter settings and the performance of our method with eight other methods on a public dataset.Moreover, we analyzed the algorithm performance for different forest types.

Study Areas
The study areas in this paper were obtained from the NEWFOR (NEW technologies for a better mountain FORest timber mobilization) project [25,26].The distribution of the study areas is shown in Figure 1.

Study Areas
The study areas in this paper were obtained from the NEWFOR (NEW technologies for a better mountain FORest timber mobilization) project [25,26].The distribution of the study areas is shown in Figure 1.NEWFOR aimed at enhancing the wood supply chain within the Alpine Space and provided a benchmark dataset within different forest types for LiDAR-based single-tree detection.There are eighteen plots across five countries in the benchmark dataset.As some of the plots were non-public, twelve plots from the benchmark dataset were chosen as our study areas.Some important forest parameters, such as the region and forest type of experimental plots are presented in Table 1.The more detailed information, such as field measurements and plot description can be found in NEWFOR project.NEWFOR aimed at enhancing the wood supply chain within the Alpine Space and provided a benchmark dataset within different forest types for LiDAR-based single-tree detection.There are eighteen plots across five countries in the benchmark dataset.As some of the plots were non-public, twelve plots from the benchmark dataset were chosen as our study areas.Some important forest parameters, such as the region and forest type of experimental plots are presented in Table 1.The more detailed information, such as field measurements and plot description can be found in NEWFOR project.

ALS Data and Reference Data
For all study plots, ALS data, DTM, as well as reference data were available in NEWFOR project, and the detailed data used in this article can be downloaded from its website [25].Since the dataset originates from different institutions, the size and shape of these experimental plots are different.Most of the twelve plots are 40 m in diameter, and some plots (i.e., 3,10,11,12) are 50 m in diameter.For the sake of comparative analysis, this study normalized twelve plots to the same size of 40 m diameter.The 3D view of plot 10 ALS data is shown as an example in Figure 2, where the colour represents the height of point.Red indicates high, and green indicates low.

ALS Data and Reference Data
For all study plots, ALS data, DTM, as well as reference data were available in NEWFOR project, and the detailed data used in this article can be downloaded from its website [25].Since the dataset originates from different institutions, the size and shape of these experimental plots are different.Most of the twelve plots are 40 m in diameter, and some plots (i.e., 3,10,11,12) are 50 m in diameter.For the sake of comparative analysis, this study normalized twelve plots to the same size of 40 m diameter.The 3D view of plot 10 ALS data is shown as an example in Figure 2, where the colour represents the height of point.Red indicates high, and green indicates low.The region and forest type of experimental plots are presented in Table 1, and more details can be seen in NEWFOR project.Among these plots, 5 plots were multi-layered mixed forests (ML/M), 2 plots were multi-layered coniferous forests (ML/C), 2 plots were single-layered mixed forests (SL/M), and 3 plots are single-layered coniferous forests (SL/C).The field measurements for plots 1, 2, 4, 5, 6, 7, and 8 occurred after the ALS data, while other plots (i.e., 3,9,10,11,12) had field measurements before the ALS data.

Methods
Existing single-tree detection algorithms require selection of the local maxima and thresholds for window dimensions.For example, the watershed method extracts the similar conical geometry of trees.However, this method must select the local maxima as seed points to avoid over-segmentation.The seed point selection depends on the manual setting of the parameters that significantly affects the experimental results.
In this study, we proposed a novel and automatic method to detect single trees from ALS data through gradient orientation clustering.Like other single-tree detection algorithms (such as watershed [19,27,28], valley-following [29,30], and active contour model [31,32]), our method considers tree canopies as a mountain-like 3D topographic surface in the rasterized CHM, as shown in Figure 3a.The gradient orientation of each point was calculated from the CHM.The gradient orientation of most points moves towards the local maxima in rasterized CHM, as shown in Figure 3b.If it starts from any point in rasterized CHM and moves along the gradient orientation of each point, it eventually reaches the highest point, which is the local maxima point.All points that reach the same highest point are clustered into the same class.Therefore, the points in the same cones are clustered together.The region and forest type of experimental plots are presented in Table 1, and more details can be seen in NEWFOR project.Among these plots, 5 plots were multi-layered mixed forests (ML/M), 2 plots were multi-layered coniferous forests (ML/C), 2 plots were single-layered mixed forests (SL/M), and 3 plots are single-layered coniferous forests (SL/C).The field measurements for plots 1, 2, 4, 5, 6, 7, and 8 occurred after the ALS data, while other plots (i.e., 3,9,10,11,12) had field measurements before the ALS data.

Methods
Existing single-tree detection algorithms require selection of the local maxima and thresholds for window dimensions.For example, the watershed method extracts the similar conical geometry of trees.However, this method must select the local maxima as seed points to avoid over-segmentation.The seed point selection depends on the manual setting of the parameters that significantly affects the experimental results.
In this study, we proposed a novel and automatic method to detect single trees from ALS data through gradient orientation clustering.Like other single-tree detection algorithms (such as watershed [19,27,28], valley-following [29,30], and active contour model [31,32]), our method considers tree canopies as a mountain-like 3D topographic surface in the rasterized CHM, as shown in Figure 3a.The gradient orientation of each point was calculated from the CHM.The gradient orientation of most points moves towards the local maxima in rasterized CHM, as shown in Figure 3b.If it starts from any point in rasterized CHM and moves along the gradient orientation of each point, it eventually reaches the highest point, which is the local maxima point.All points that reach the same highest point are clustered into the same class.Therefore, the points in the same cones are clustered together.The process of the proposed method consists of three steps: (1) calculate the DSM based on ALS data by retaining the highest altitude value of the points in each pixel, and then, the CHM is derived by subtracting the DTM from the DSM; (2) obtain the potential single-trees from gradient orientation clustering; and (3) select the final individual trees based on the density feature and determine its coordinates.

Obtain the CHM from ALS Data
The point cloud in ALS data was irregularly distributed.However, the ALS height values can be summarized into regular grids.Figure 4    The value will be set to 0 if there is no ALS-information available in a raster cell, which we called black holes.Usually, there are some black holes in the rasterized DSM, which are caused by the abnormal changes of the height values.This unnatural phenomenon may cause potential errors in the treetop detection, but we can solve it through interpolation.Interpolation algorithm is a commonly used algorithm in image processing [33,34].Usually, there are linear, cubic and nearestneighbor interpolation methods, and they are integrated into some image processing softwares, such as Matlab (The MathWorks Inc., Natick, MA, USA).What we used is the cubic interpolation method.The process of the proposed method consists of three steps: (1) calculate the DSM based on ALS data by retaining the highest altitude value of the points in each pixel, and then, the CHM is derived by subtracting the DTM from the DSM; (2) obtain the potential single-trees from gradient orientation clustering; and (3) select the final individual trees based on the density feature and determine its coordinates.

Obtain the CHM from ALS Data
The point cloud in ALS data was irregularly distributed.However, the ALS height values can be summarized into regular grids.Figure 4 is an example of the process of rasterization ALS point cloud data to DSM.In this example, ∆x = 5 dm and ∆y = 5 dm are the length and width of the grid, respectively, and r = 6 dm is the search radius of the rasterized pixels.There may be more than one point in the same search range.We set the value of the rasterized pixels to the maximum height values within the search range.The process of the proposed method consists of three steps: (1) calculate the DSM based on ALS data by retaining the highest altitude value of the points in each pixel, and then, the CHM is derived by subtracting the DTM from the DSM; (2) obtain the potential single-trees from gradient orientation clustering; and (3) select the final individual trees based on the density feature and determine its coordinates.

Obtain the CHM from ALS Data
The point cloud in ALS data was irregularly distributed.However, the ALS height values can be summarized into regular grids.Figure 4    The value will be set to 0 if there is no ALS-information available in a raster cell, which we called black holes.Usually, there are some black holes in the rasterized DSM, which are caused by the abnormal changes of the height values.This unnatural phenomenon may cause potential errors in the treetop detection, but we can solve it through interpolation.Interpolation algorithm is a commonly used algorithm in image processing [33,34].Usually, there are linear, cubic and nearestneighbor interpolation methods, and they are integrated into some image processing softwares, such as Matlab (The MathWorks Inc., Natick, MA, USA).What we used is the cubic interpolation method.The value will be set to 0 if there is no ALS-information available in a raster cell, which we called black holes.Usually, there are some black holes in the rasterized DSM, which are caused by the abnormal changes of the height values.This unnatural phenomenon may cause potential errors in the treetop detection, but we can solve it through interpolation.Interpolation algorithm is a commonly used algorithm in image processing [33,34].Usually, there are linear, cubic and nearest-neighbor interpolation methods, and they are integrated into some image processing softwares, such as Matlab (The MathWorks Inc., Natick, MA, USA).What we used is the cubic interpolation method.
As shown in Figure 5, the CHM can be directly subtracted from DSM to DTM.The value of each raster cell in the CHM is the absolute height above ground.The values of some raster cells in the calculated CHM may be less than zero, and those raster cells will be replaced using nearest-neighbour interpolation [35].Then, a 3 × 3 Gaussian filtering with σ = 0.5 was used for image smoothing and denoising of the CHM.Section 4.4 will discuss the selection of Gaussian filtering in detail.As shown in Figure 5, the CHM can be directly subtracted from DSM to DTM.The value of each raster cell in the CHM is the absolute height above ground.The values of some raster cells in the calculated CHM may be less than zero, and those raster cells will be replaced using nearest-neighbour interpolation [35].Then, a 3  3 Gaussian filtering with σ = 0.5 was used for image smoothing and denoising of the CHM.Section 4.4 will discuss the selection of Gaussian filtering in detail.

Gradient Orientation Clustering
The x-gradient-component G x (x, y) of point (x, y) in the CHM can be calculated using the Sobel operator [36]: The y-gradient-component G y (x, y) is The gradient orientation θ(x, y) is The pseudocode of 4-N GOC algorithm (Algorithm 1) is shown as follows.
for y = 1 → N // Traversing every point in the matrix 3.
(x c , y c ) ← (x, y) // Select a current point from the matrix 4.
for (i, j) in map 12.
label(i, j) ← index // Set every point in map In the 4-N GOC algorithm, label is a label matrix with the size of CHM, which is initialized to 0. These points in label have a same value belong to the same class, (e.g., label(x, y) = k indicates that point (x, y) belongs to class k).index is initialized to 1, and when it increments by 1 indicates that one cluster is completed.map is a point set that stores the points belong to the same class in clustering, and (i, j) represents a point in the map.(x c , y c ) is the point being processed currently.∠θ(α, β) is the angle between the gradient orientation θ(x c , y c ) and the orientation from (x c , y c ) to (α, β), where (α, β) represents one neighbour point of (x c , y c ). (x , y ) is the neighbour point of (x c , y c ) that has the minimum ∠θ(α, β).
The gradient orientation clustering (GOC) algorithm is a procedure to find a path from one point to the last point of the cluster.According to the neighbour size, GOC can be divided into 4-N (4-Neighbour) and 8-N (8-Neighbour) GOC algorithms.Compared with Pseudocode of 4-N GOC, Pseudocode of 8-N GOC is only different in line 9 where Figure 6 is an example of the process of GOC.During the clustering, both 4-N and 8-N GOC algorithms, the point that is nearest to the current point at the gradient direction is added to the clustering path.
Figure 6 is an example of the process of GOC.During the clustering, both 4-N and 8-N GOC algorithms, the point that is nearest to the current point at the gradient direction is added to the clustering path.

Single-Tree Selection
In CHM, the gradient orientation of the shrubs and edges of canopies are chaotic.Thus, after the process of gradient orientation clustering, there is some noise and obviously erroneous detection objects in the image, such as rough borders and small plaques.This noise and erroneous detection objects can be eliminated by performing an image morphology opening and closing operation.We had performed some experiments for the morphology opening and closing operation with different structural elements.Finally, the 5  5 square structuring element is selected, as the 3  3 square has a poor smoothing effect, and the 7  7 square may remove some trees correctly detected and distort the image.Although most objects removed are noise and erroneous in the morphology operation, it is inevitable to remove some small trees.For example, if the resolution of CHM is 0.5 m, we will not extract the crown less than 2 m 2 , because the object less than 3  3 cells was removed.
The shape of the tree canopy should approximate a circle.As shown in Figure 7, the object of region 2 is more such as a tree than region 1.The features of the image object that do not satisfy the tree canopy shape must be removed.By defining the density value to describe the compactness of the regions, the noise region (with an abnormal density value) can be removed.The density d is expressed as follows [37]: where n is the number of image pixels for the object; X and Y are the x and y coordinates of all image pixels that constitute the object; and Var is the variance operation.Density d describes the degree of compactness of an object.An image object commonly resembles a square with high density.
We defined the density threshold Td as follows: where res is the image resolution,

Single-Tree Selection
In CHM, the gradient orientation of the shrubs and edges of canopies are chaotic.Thus, after the process of gradient orientation clustering, there is some noise and obviously erroneous detection objects in the image, such as rough borders and small plaques.This noise and erroneous detection objects can be eliminated by performing an image morphology opening and closing operation.We had performed some experiments for the morphology opening and closing operation with different structural elements.Finally, the 5 × 5 square structuring element is selected, as the 3 × 3 square has a poor smoothing effect, and the 7 × 7 square may remove some trees correctly detected and distort the image.Although most objects removed are noise and erroneous in the morphology operation, it is inevitable to remove some small trees.For example, if the resolution of CHM is 0.5 m, we will not extract the crown less than 2 m 2 , because the object less than 3 × 3 cells was removed.
The shape of the tree canopy should approximate a circle.As shown in Figure 7, the object of region 2 is more such as a tree than region 1.The features of the image object that do not satisfy the tree canopy shape must be removed.By defining the density value to describe the compactness of the regions, the noise region (with an abnormal density value) can be removed.The density d is expressed as follows [37]: where n is the number of image pixels for the object; X and Y are the x and y coordinates of all image pixels that constitute the object; and Var is the variance operation.Density d describes the degree of compactness of an object.An image object commonly resembles a square with high density., the image object is considered a tree.For example, the density of region 2 in Figure 7a is 1.92, and the density of region 1 in Figure 7a is 1.12.Therefore, the object in region 2 is more likely a tree.
The final step is to determine the centre points of the trees.In this paper, the centre of the minimum circumscribed circle around the image object is considered the centre of a tree.The altitude of the centre point is considered the tree height.The radius of the minimum circumscribed circle is considered the crown radius of a tree.The effect of density features is detailly discussed in Section 4.1.

Parameters for Accuracy Evaluation
To prove the effectiveness of this method, we analysed the test results of this study with the reference data.The detected tree and reference tree were considered a match if the horizontal difference ∆ and vertical difference ∆ were within a certain range.The matching process is as follows.(1) Candidate tree selection: for a test tree, the reference trees are added to the candidate set if the horizontal difference ∆ and vertical difference ∆ are within a certain threshold We defined the density threshold T d as follows: where res is the image resolution, k 1 and k 2 are two constants, and T d ∈ (1, 2).A higher image resolution corresponds to a larger tree density in the image.Therefore, the appropriate k 1 and k 2 must be selected to obtain the threshold at different resolutions.k 1 is the initial value of a threshold, and . Use k 2 to adjust the threshold for different resolutions, where k 2 ∈ [0, 1].When the density of the image object d > T d , the image object is considered a tree.For example, the density of region 2 in Figure 7a is 1.92, and the density of region 1 in Figure 7a is 1.12.Therefore, the object in region 2 is more likely a tree.The final step is to determine the centre points of the trees.In this paper, the centre of the minimum circumscribed circle around the image object is considered the centre of a tree.The altitude of the centre point is considered the tree height.The radius of the minimum circumscribed circle is considered the crown radius of a tree.The effect of density features is detailly discussed in Section 4.1.

Parameters for Accuracy Evaluation
To prove the effectiveness of this method, we analysed the test results of this study with the reference data.The detected tree and reference tree were considered a match if the horizontal difference ∆H 2H and vertical difference ∆V 2V were within a certain range.The matching process is as follows.
(1) Candidate tree selection: for a test tree, the reference trees are added to the candidate set if the horizontal difference ∆H 2H and vertical difference ∆V 2V are within threshold range.If the tree is higher than 15 m, we set ∆H 2H < 5 m and ∆V 2V < 2 m, else we set ∆H 2H < 4 m and ∆V 2V < 1.5 m; (2) Choose the best candidate tree: find the nearest reference tree to the test tree from the candidate set as the best candidate tree; (3) Candidate testing: the matching problem is not a one-way problem.A test tree needs to find the best matching reference tree.The reference tree also needs to find the best matching test tree.These two trees are considered as a successful match only when the candidate tree of the test results and the candidate tree of reference data are two candidate trees for each other.
The evaluation parameters and their calculation methods are defined as follows [26]: These evaluation parameters were used to measure the matching results for each experimental dataset.The extraction rate should be within a reasonable range.If the extraction rate is too high, (e.g., R extr = 200%, that indicates at least half of the extracted trees cannot match reference trees, or many extracted trees are incorrect).A higher matching rate indicates a better single-tree detection method.A higher commission rate and omission rate indicate a worse single-tree detection method.Since each above-mentioned indicator can only show one aspect of algorithm performance, we measure the performance of the algorithm more comprehensive and directly by using a matching score.The idea of using a matching score comes from Larsen's work [38].In Larsen's work the matching score is M Laesen = 100 × N matching /(N matching + N com + N om ).In order to make it easier to compare with the previous studies, we modified the matching score as below.

•
M (Matching score): Although the expressions of these two formulas are different, a high matching score under our criteria can also indicate that the algorithm performance is better.In addition, root mean square (RMS) was used to evaluate the accuracy of model predictions.The evaluation parameters of the RMS were defined as follows [26]: The RMS of X is defined as Formula ( 7), where X can be H mean , V mean , R extr , R mat , R com , R om and M, i represents the i-th plot.

Experiment & Discussion
This section aims to verify the effectiveness of our algorithm by conducting some related experiments.Firstly, it shows the selection of appropriate parameters, such as 4-N or 8-N GOC, different resolutions, different Gaussian filters, the values of k 1 , k 2 and T d .Next, it shows the results of 4-N GOC and 8-N GOC for each plot, and evaluates the detection results by using the evaluation criteria, such as extraction rate, matching rate, commission rate, omission rate, mean horizontal distance and mean height distance.Then, to show the advantages of our method, it presents the experiment results of RMS_R extr , RMS_R mat , RMS_R com , RMS_R om , RMS_H, and RMS_V for 12 plots to compare the detection results of our method with other existing methods.Finally, it presents the applicability of our algorithm for different forest types.

Analysis of Model Parameters
This section mainly presents the selection of experimental parameters.It first shows the effects of using different Gaussian filters.Then, it shows the effects of density threshold.Finally, it presents the effects of 4-N and 8-N and the effects of different resolutions of rasterized CHM.
We have tried several filter-specifications such as 3 × 3, 5 × 5 and 7 × 7 filtering to smoothen the CHM.With the same window size, there is almost no difference in the performances for different σ values, so in our experiment we chose σ = 0.5.
Table 2 shows matching scores of all 12 plots and there RMS_M of no filter, 3 × 3, 5 × 5 and 7 × 7 Gaussian filters.Moreover, Figure 8 shows 0.5 m CHM and the matching results of no filter, 3 × 3, 5 × 5 and 7 × 7 Gaussian filters in plot 10 using 4-N GOC method as an example.As we can see, the matching score M of no filter, 3 × 3, 5 × 5 and 7 × 7 Gaussian filters in plot 10 are 39, 45 45, 45.For plot 10, the detection result for the CHM using filter was better than the CHM without filter, and the single tree detection using different gaussian filters almost had the same results.Moreover, RMS_M of all these twelve plots showed that the detection result for the CHM using filter will have better performance than the CHM without filter.For these filters, the 3 × 3 and 5 × 5 filter had a better performance.For further experiments, we also have implemented the related algorithms.For example, we compared the matching results using our method with LM + filtering [39].Figure 11 shows the matching results for 4-N GOC and the LM + filtering method.Figure 10a shows the matching results of 4-N GOC with 0.5 m resolution in plot 10 using 4-N GOC method.The extraction rate is 170%, the matching rate is 76%, the commission rate is 56%, the  For the impact of density threshold T, we also took plot 10 as an example.Table 3 shows the matching results of plot 10 with T d of 0, 1.2, 1.4, and 1.6 respectively.The extracted number decreases quickly when the threshold increases, and the matched number decreases rarely.Simply said, the density threshold T d can help reduce the commission rate.Moreover, R om does not increase with the same rate as R com decreases.Therefore, the matching score M by setting the density threshold T d will have better performance.When the threshold is within a certain range (1.2 to 1.6), the impact on the result is not significant.In addition, we used two detection methods to present the experimental results: 4-N and 8-N gradient orientation clustering.The resolution of the CHM significantly affects the experimental results.If the density of the ALS data is not sufficient to support a high-resolution CHM, the CHM will have many null points.Although they can be filled using the interpolation method, this problem inevitably increases the noise in the CHM.If the CHM resolution is too low, the information in the CHM to compute the gradient orientation is inadequate and inaccurate.Then, the detection rate decreases; in particular, some small trees cannot be detected.To compare the experimental results, we generated three different resolution (0.3 m, 0.5 m, 0.7 m) CHM in this paper.Table 4 shows different parameters and the experimental results for all datasets.The three groups of k 1 and k 2 have notably similar experiment results, which implies that the parameters have a minimal effect on the single-tree matching result.The best matching result is obtained when k 1 is 1.55 and k 2 is 0.5.Therefore, in our experiments, the values of k 1 k 2 can be set as 1.55 and 0.5, respectively.We do not need manually adjust the values of these parameters according to the spatial resolutions of remote-sensing images, which enable the automation of single-tree detection.
Table 5 shows the 4-N GOC and 8-N GOC matching results for CHM at 0.3 m, 0.5 m, and 0.7 m resolutions over the entire dataset.The 4-N GOC single-tree detection method achieved relatively satisfactory results at 0.5 m resolution.If the resolution of CHM is too high or too low, the results will not be reasonable.If the CHM resolution is too high, the point density is sparse, which results in more noise in the CHM.If the CHM resolution is too low, too much gradient information is lost in the CHM.As to the reason that 4-N GOC performs better than 8-N GOC in our algorithm, the remote sensing image is not clear enough, so the gradient information is partly incorrect.In 8-N GOC, one raster cells will have eight directions to point to, and there are less raster cells point to a same one compared with 4-N GOC.Therefore, the clustering of 8-N is easy to converge prematurely, that results in an error clustering, which can be verified from the commission rate.As we can see in Table 5, the commission rate of 8-N GOC is higher than that of 4-N GOC no matter what the extraction rate is.Therefore, we finally adopted 4-N GOC in our comparative experiment.

Results for Each Plot of the Dataset
This section gives the results of each plot and the evaluation results of relevant six parameters (extraction rate, matching rate, commission rate, omission rate, mean horizontal distance and mean height distance) of 4-N GOC and 8-N GOC for each plot with the 0.3 m, 0.5 m and 0.7 m CHM.
This study conducted the proposed single-tree detection algorithm in twelve plots of approximately 20 metres in diameter, which belonged to multi-layered mixed forests (ML/M), multi-layered coniferous forests (ML/C), single-layered mixed forests (SL/M), and single-layered coniferous forests (SL/C), as shown in Table 1. Figure 9 illustrates the results of the extracted individual trees by the 4-N GOC method with 0.5 m CHM.As to the reason that 4-N GOC performs better than 8-N GOC in our algorithm, the remote sensing image is not clear enough, so the gradient information is partly incorrect.In 8-N GOC, one raster cells will have eight directions to point to, and there are less raster cells point to a same one compared with 4-N GOC.Therefore, the clustering of 8-N is easy to converge prematurely, that results in an error clustering, which can be verified from the commission rate.As we can see in Table 5, the commission rate of 8-N GOC is higher than that of 4-N GOC no matter what the extraction rate is.Therefore, we finally adopted 4-N GOC in our comparative experiment.

Results for Each Plot of the Dataset
This section gives the results of each plot and the evaluation results of relevant six parameters (extraction rate, matching rate, commission rate, omission rate, mean horizontal distance and mean height distance) of 4-N GOC and 8-N GOC for each plot with the 0.3 m, 0.5 m and 0.7 m CHM.
This study conducted the proposed single-tree detection algorithm in twelve plots of approximately 20 metres in diameter, which belonged to multi-layered mixed forests (ML/M), multilayered coniferous forests (ML/C), single-layered mixed forests (SL/M), and single-layered coniferous forests (SL/C), as shown in Table 1. Figure 9 illustrates the results of the extracted individual trees by the 4-N GOC method with 0.5 m CHM.In this study, there were 1286 reference trees, and a total of 1189 trees were extracted using the proposed algorithm in these twelve plots.The forest type of plots 1 and 3 were ML/C, and plots 2, 4, and 6 were SL/C.They were composed of coniferous trees, and the gradient of a coniferous tree usually points to only one treetop, so these plots had good extracted rates.The main problem of these In this study, there were 1286 reference trees, and a total of 1189 trees were extracted using the proposed algorithm in these twelve plots.The forest type of plots 1 and 3 were ML/C, and plots 2, 4, and 6 were SL/C.They were composed of coniferous trees, and the gradient of a coniferous tree usually points to only one treetop, so these plots had good extracted rates.The main problem of these plots was the omission rate.The forest type of plots 5, 7, 8, 11, and 12 were ML/M, and that of plots 9 and 10 was SL/M.They are mixed forests of conifer and broadleaf trees.As the gradient of a broadleaf tree may point to multiple treetops, the commission rate of these plots was slightly high.Overall, the proposed algorithm extracted most individual treetops in four types of forests.The matching score was M = 39 for ML/C plots, M = 51 for SL/C plots and M = 39 for ML/M plots, and M = 47 for SL/M plots.However, there were some omitted and misrecognized trees in each plot.We will discuss the matching results for different forest types in Section 4.4.
In Figure 10a, the extraction rate is higher when the resolution of the CHM resolution is higher, but the extraction rate being too high is not necessarily a good thing.It will also result in a high commission rate.The extraction rate should be in a suitable range (e.g., 50% to 150%).Then, we counted the plot number of 4-N and 8-N GOC's extraction rate in a suitable range in total of 0.3 m, 0.5 m and 0.7 m are 11, 17, and 19, respectively.
A higher resolution corresponds to a higher matching rate because the extraction rate is high at high resolution, where more trees are accurately extracted at the same time, and the matching rate increases.At low resolution, the extraction rate is low and the number of accurately extracted trees decreases.Therefore, the matching rate in the single-tree detection is low.A higher matching rate is better, but we think that a result with greater than 60% the matching rate is a relatively good result.As in Figure 10b, 4-N GOC and 8-N GOC in 0.5 m CHM achieved the most plots with an excellent result.There were 10 plots of 4-N GOC with excellent results at 0.5 m resolution.These 10 plots include all single-layered forests (plots 2, 4, 6, 9 and 10), which indicates that the method achieves good results for single-layered forest types.
The commission is the misrecognized trees in the detection results; thus, the commission rate increase with the extraction rate increase.As seen in Figure 10c, the commission rates in 0.3 m CHM are higher than in 0.5 m and 0.7 m CHM.Because the density of ALS points in a pixel of the rasterized CHM of 0.3 m resolution is low, the 0.3 m CHM contains more noise; therefore, there are more misrecognized trees in the detection results.The omission is leakage in the reference trees, so the omission rate increases with the extraction rate decrease.As seen in Figure 10d, the omission rates in 0.7 m CHM are higher than that in 0.3 m and 0.5 m CHM.Too much gradient information is lost when generating CHM with 0.7 m resolution with ALS point cloud data, and the trees density feature in low resolution CHM is difficult to recognize.Figure 10e,f show the mean horizontal and vertical distance between the matched trees and reference trees.The mean horizontal distance of 0.3 m, 0.5 m, 0.7 m CHM is 2.0 m, 1.9 m and 1.9 m respectively.The mean vertical distance of 0.3 m, 0.5 m, 0.7 m CHM is 1.1 m, 1.0 m and 0.9 m, respectively.The mean horizontal and vertical distance difference between CHM of the three resolutions is small.Through the above histogram of six evaluation parameters, we find that the proposed 4-N GOC and 8-N GOC achieve best result in CHM of 0.5 m resolution.plots was the omission rate.The forest type of plots 5, 7, 8, 11, and 12 were ML/M, and that of plots 9 and 10 was SL/M.They are mixed forests of conifer and broadleaf trees.As the gradient of a broadleaf tree may point to multiple treetops, the commission rate of these plots was slightly high.Overall, the proposed algorithm extracted most individual treetops in four types of forests.The matching score was M = 39 for ML/C plots, M = 51 for SL/C plots and M = 39 for ML/M plots, and M = 47 for SL/M plots.However, there were some omitted and misrecognized trees in each plot.We will discuss the matching results for different forest types in Section 4.5.
In Figure 10a, the extraction rate is higher when the resolution of the CHM resolution is higher, but the extraction rate being too high is not necessarily a good thing.It will also result in a high commission rate.The extraction rate should be in a suitable range (e.g., 50% to 150%).Then, we counted the plot number of 4-N and 8-N GOC's extraction rate in a suitable range in total of 0.3 m, 0.5 m and 0.7 m are 11, 17, and 19, respectively.
A higher resolution corresponds to a higher matching rate because the extraction rate is high at high resolution, where more trees are accurately extracted at the same time, and the matching rate increases.At low resolution, the extraction rate is low and the number of accurately extracted trees decreases.Therefore, the matching rate in the single-tree detection is low.A higher matching rate is better, but we think that a result with greater than 60% the matching rate is a relatively good result.As in Figure 10b, 4-N GOC and 8-N GOC in 0.5 m CHM achieved the most plots with an excellent result.There were 10 plots of 4-N GOC with excellent results at 0.5 m resolution.These 10 plots include all single-layered forests (plots 2, 4, 6, 9 and 10), which indicates that the method achieves good results for single-layered forest types.
The commission is the misrecognized trees in the detection results; thus, the commission rate increase with the extraction rate increase.As seen in Figure 10c, the commission rates in 0.3 m CHM are higher than in 0.5 m and 0.7 m CHM.Because the density of ALS points in a pixel of the rasterized CHM of 0.3 m resolution is low, the 0.3 m CHM contains more noise; therefore, there are more misrecognized trees in the detection results.The omission is leakage in the reference trees, so the omission rate increases with the extraction rate decrease.As seen in Figure 10d, the omission rates in 0.7 m CHM are higher than that in 0.3 m and 0.5 m CHM.Too much gradient information is lost when generating CHM with 0.7 m resolution with ALS point cloud data, and the trees density feature in low resolution CHM is difficult to recognize.Figure 10e,f show the mean horizontal and vertical distance between the matched trees and reference trees.The mean horizontal distance of 0.3 m, 0.5 m, 0.7 m CHM is 2.0 m, 1.9 m and 1.9 m respectively.The mean vertical distance of 0.3 m, 0.5 m, 0.7 m CHM is 1.1 m, 1.0 m and 0.9 m, respectively.The mean horizontal and vertical distance difference between CHM of the three resolutions is small.Through the above histogram of six evaluation parameters, we find that the proposed 4-N GOC and 8-N GOC achieve best result in CHM of 0.5 m resolution.

Comparison with Previous Studies
To verify the effectiveness of the proposed GOC method, we conducted a comparative experiment with eight previous single-tree detection methods.The experimental results (i.e., RMS_Rextr, RMS_Rmat, RMS_Rcom, RMS_Rom, RMS_H, RMS_V) of these eight methods are provided in the study of Eysn et al. [26].Although these indicators are very specific, it is difficult to directly know the effect of the algorithm.Therefore, we calculated the root mean square of matching scores (RMS_M), which can directly compare the performance of these algorithms.In order to make it easier to compare with the other algorithms compared by Eysn et al., we modified Larsen's formula of matching score.Although the expressions of these two formulas are different, the ideas of computing matching score are basically consistent.Moreover, our method can produce similar results when Nref = Ntest.
Table 6 shows the matching results of various single-tree detection methods.The ninth method is our method with 4-N GOC at 0.5 m resolution.The single-tree detection method with 4-N GOC can achieve a 63% matching rate at 124% extraction rate, which is much higher than the other detection methods.The commission rate in our method is 42%, which is lower than those of methods 3, 4, 5, and 6.The omission rate in our method is 42%, which is lower than any other method.RMS_H is 1.9 m, which is slightly higher than those of other methods; RMS_V is 1.0 m, which can reach the average level of all other methods.For the RMS_M, our single-tree detection method of 4-N GOC achieves a maximum RMS_M over the entire dataset.Therefore, the 4-N GOC is an excellent singletree detection method to detect single trees from the ALS data.

Comparison with Previous Studies
To verify the effectiveness of the proposed GOC method, we conducted a comparative experiment with eight previous single-tree detection methods.The experimental results (i.e., RMS_R extr , RMS_R mat , RMS_R com , RMS_R om , RMS_H, RMS_V) of these eight methods are provided in the study of Eysn et al. [26].Although these indicators are very specific, it is difficult to directly know the effect of the algorithm.Therefore, we calculated the root mean square of matching scores (RMS_M), which can directly compare the performance of these algorithms.In order to make it easier to compare with the other algorithms compared by Eysn et al., we modified Larsen's formula of matching score.Although the expressions of these two formulas are different, the ideas of computing matching score are basically consistent.Moreover, our method can produce similar results when N ref = N test .
Table 6 shows the matching results of various single-tree detection methods.The ninth method is our method with 4-N GOC at 0.5 m resolution.The single-tree detection method with 4-N GOC can achieve a 63% matching rate at 124% extraction rate, which is much higher than the other detection methods.The commission rate in our method is 42%, which is lower than those of methods 3, 4, 5, and 6.The omission rate in our method is 42%, which is lower than any other method.RMS_H is 1.9 m, which is slightly higher than those of other methods; RMS_V is 1.0 m, which can reach the average level of all other methods.For the RMS_M, our single-tree detection method of 4-N GOC achieves a maximum RMS_M over the entire dataset.Therefore, the 4-N GOC is an excellent single-tree detection method to detect single trees from the ALS data.For further experiments, we also have implemented the related algorithms.For example, we compared the matching results using our method with LM + filtering [39].Figure 11 shows the matching results for 4-N GOC and the LM + filtering method.For further experiments, we also have implemented the related algorithms.For example, we compared the matching results using our method with LM + filtering [39].Figure 11 shows the matching results for 4-N GOC and the LM + filtering method.Figure 10a shows the matching results of 4-N GOC with 0.5 m resolution in plot 10 using 4-N GOC method.The extraction rate is 170%, the matching rate is 76%, the commission rate is 56%, the Figure 11a shows the matching results of 4-N GOC with 0.5 m resolution in plot 10 using 4-N GOC method.The extraction rate is 170%, the matching rate is 76%, the commission rate is 56%, the omission rate is 24%, and the matching score is 49. Figure 11b,c use the LM + filtering.In Figure 11b, the window size is 5, the extraction rate is 91%, the matching rate is 59%, the commission rate is 35%, the omission rate is 41%, and the matching score is 44.In Figure 11c, the window size is 3, the extraction rate is 329%, the matching rate is 81%, the commission rate is 75%, the omission rate is 19%, and the matching score is 46.The results show that the 4-N GOC method can provide superior results, whereas the LM + filtering method is strongly affected by the window size.

Analysis of Matching Results for Different Forest Types
The existing single-tree detection methods based on rasterized CHM are commonly used to detect trees in a single-layered forest or dominant trees in a multi-layered forest.It is difficult to effectively detect suppressed trees in a multi-layered forest.Though GOC cannot detect suppressed trees in a multi-layered forest either, it can appropriately capture the dominant trees at least.
To explore the effect of the GOC algorithm on different types of forests, we applied our method to detect single trees in both single-layered forests and multi-layered forests.Table 7 shows the matching results for different types of forests using our 4-N GOC method at 0.5 m resolution.According to the matching results, there are larger values of RMS_R extr , RMS_R mat and RMS_R com in single-layered mixed forests and single-layered coniferous forests.The values of RMS_M are 47 and 51 for the single-layered mixed forests and single-layered coniferous forests, respectively.Therefore, our single-tree detection method is more suitable for non-overlapping trees in a single-layered forest.
Because CHM only contains the top layer information of the tree canopy, that leads to better performance of the algorithm in single-layer forests than multi-layer forests.The matching results of multi-layered mixed forests and multi-layered coniferous forests have smaller values of RMS_R extr , RMS_R mat and RMS_R com than those of the single-layered forest.The values of RMS_M for the multi-layered mixed forest and multi-layered coniferous forest is 39.Although the matching results of the multi-layered forest are not as satisfactory as those of single-layered forests, they are still better than most the other existing methods in Table 6.Therefore, GOC method can also be used to automatically and effectively identify dominant trees in multi-layered forests.SL/M: single-layered mixed forests; SL/C: single-layered coniferous forests; ML/M: multi-layered mixed forests; ML/C: multi-layered coniferous forests.

Conclusions
This paper proposed an automatic method for single-tree detection in airborne laser scanning data through gradient orientation clustering.The proposed method combines the height value of CHM and gradient orientation information to cluster points belonging to the same crown as a conical shape in a canopy to detect single-trees.In the experiments, we used the proposed method with two different neighbour sizes (4-N and 8-N), and we presented the experimental results in six different aspects, the extraction rate, matching rate, commission rate, omission rate, mean horizontal distance and mean height distance.In our method, the CHM resolution significantly affects the experimental results.If the ALS points density is not sufficient to support a high-resolution CHM, the CHM may have some pits.Although they can be filled by using the linear interpolation method, this problem inevitably distorts the CHM.If the CHM resolution is too low, the gradient orientation information of the CHM is inadequate and inaccurate.This result in the extraction rate and matching rate obviously decreasing and many small trees cannot be detected.Experiment shows that the overall matching scores of the proposed GOC method in four types of forests are higher than those of the other eight methods, and the root mean square matching score of GOC method (RMS_M is 43) at least 7% higher than the other eight methods.Although the proposed method can detect dominant single trees in multi-layered forests, this method has more satisfactory results in single-layered forests.
As to the types of forest communities, such as the broad-leaved forest, coniferous forest, and mixed coniferous broad-leaved forest, as long as the tree crown has a conical shape, our method is applicable.However, the single-tree detection method still has room for improvement, particularly for the detection of multi-layer forest, as the information for different layers of the multi-layered forest is insufficient in the CHM.Because of GOC method needs to calculate the gradient orientation of the trees and small trees' gradient orientation have deviation, it is not suitable for small trees' detection.For example, GOC cannot correctly detect the crown whose projection area is below 2 m 2 .In the future, our method will be combined with a point cloud-based method to detect small trees and find a way to detect the suppressed trees in multi-layered forests.

Figure 1 .
Figure 1.Study areas distribution.Red triangles represent different study areas that this paper used, and yellow triangles represent the remaining areas of NEWFOR (NEW technologies for a better mountain FORest timber mobilization) project.
(a) Top view (b) Side view
is an example of the process of rasterization ALS point cloud data to DSM.In this example, the length and width of the grid, respectively, and 6 dm r  is the search radius of the rasterized pixels.There may be more than one point in the same search range.We set the value of the rasterized pixels to the maximum height values within the search range.

Figure 4 .
Figure 4.The process of rasterizing ALS point cloud data to the digital surface model (DSM).The blue point represent the ALS point.The black rectangles represent the rasterized pixels and the red cross represents the value of a rasterized pixel.
is an example of the process of rasterization ALS point cloud data to DSM.In this example, the length and width of the grid, respectively, and 6 dm r  is the search radius of the rasterized pixels.There may be more than one point in the same search range.We set the value of the rasterized pixels to the maximum height values within the search range.

Figure 4 .
Figure 4.The process of rasterizing ALS point cloud data to the digital surface model (DSM).The blue point represent the ALS point.The black rectangles represent the rasterized pixels and the red cross represents the value of a rasterized pixel.

Figure 4 .
Figure 4.The process of rasterizing ALS point cloud data to the digital surface model (DSM).The blue point represent the ALS point.The black rectangles represent the rasterized pixels and the red cross represents the value of a rasterized pixel.

Figure 5 .
Figure 5. Computing canopy height within each pixel based on the difference between the surface model and the terrain model.DSM = digital surface model; DTM = digital terrain model; CHM = canopy height model.

Figure 5 .
Figure 5. Computing canopy height within each pixel based on the difference between the surface model and the terrain model.DSM = digital surface model; DTM = digital terrain model; CHM = canopy height model.

Figure 6 .
Figure 6.Example of gradient orientation clustering (GOC).The blue arrow represents the gradient orientation of a point and the red squares represent the clustering path.

1 k 2 kFigure 6 .
Figure 6.Example of gradient orientation clustering (GOC).The blue arrow represents the gradient orientation of a point and the red squares represent the clustering path.

Figure 7 .
Figure 7.An example of single-tree selection based on density feature.The green pixels and the red pixels represent region 1 and region 2, respectively.

Figure 7 .
Figure 7.An example of single-tree selection based on density feature.The green pixels and the red pixels represent region 1 and region 2, respectively.

Figure 11 .
Figure 11.Matching results of 4-N GOC and LM + filtering methods.

Figure 9 .
Figure 9. Results for each plot of the dataset.

Figure 9 .
Figure 9. Results for each plot of the dataset.

Figure 10 .
Figure 10.The evaluation parameters for each plot of the study plots (the raster resolutions are 0.3, 0.5, and 0.7 m respectively, and the numbers indicate the plot IDs).

Figure 11 .
Figure 11.Matching results of 4-N GOC and LM + filtering methods.

Table 1 .
Region and forest type of experimental plots.
ML/C: multi-layered coniferous forests; SL/C: single-layered coniferous forests; ML/M: multi-layered mixed forests; SL/M: single-layered mixed forests.N Figure 1.Study areas distribution.Red triangles represent different study areas that this paper used, and yellow triangles represent the remaining areas of NEWFOR (NEW technologies for a better mountain FORest timber mobilization) project.

Table 1 .
Region and forest type of experimental plots.
R com = (N test − N match )/N test , where commission represents misrecognized of the test result.• R om (Omission rate): R om = (N ref − N match )/N ref , where omission represents leakage recognition of the reference data.
mean : Mean tree height difference of the matching trees • R extr (Extraction rate): R extr = N test /N ref • R mat (Matching rate): R mat = N match /N ref • R com (Commission rate): Root mean square of H mean .• RMS_V: Root mean square of V mean .• RMS_R extr : Root mean square of R extr .• RMS_R mat : Root mean square of R mat .• RMS_R com : Root mean square of R com .• RMS_R om : Root mean square of R om .
• RMS_M: Root mean square of matching score.

Table 6 .
Matching results of single-tree detection methods.Our single-tree detection method with 4-N GOC achieves a maximum RMS_M over the dataset.

Table 2 .
Matching Scores of all Plots and their RMS_M for different Gaussian filters.
RMS_M: root mean square matching score.

Table 3 .
Matching results of plot 10 with density threshold T d of 0, 1.2, 1.4, and 1.6 respectively.

Table 4 .
Effects of different parameters on the matching results in single-tree detection.

Table 5 .
Matching results of 4-N gradient orientation clustering (GOC) and 8-N GOC at different resolutions.4-N GOC achieves more satisfactory results at 0.5 m resolution.

Table 6 .
Matching results of single-tree detection methods.Our single-tree detection method with 4-N GOC achieves a maximum RMS_M over the dataset.

Table 6 .
Matching results of single-tree detection methods.Our single-tree detection method with 4-N GOC achieves a maximum RMS_M over the dataset.

Table 7 .
Matching results for different forest types.