A Region-Based Hierarchical Cross-Section Analysis for Individual Tree Crown Delineation Using ALS Data

In recent years, airborne Light Detection and Ranging (LiDAR) that provided three-dimensional forest information has been widely applied in forest inventory and has shown great potential in automatic individual tree crown delineation (ITCD). Usually, ITCD algorithms include treetop detection and crown boundary delineation procedures. In this study, we proposed a novel method called region-based hierarchical cross-section analysis (RHCSA), which combined the two procedures together based on a canopy height model (CHM) derived from airborne LiDAR data for ITCD. This method considers the CHM as a three-dimensional topological surface, simulates stereoscopic scanning from top to bottom using an iterative process, and utilizes the individual crown and vertical structure of crowns to progressively detect individual treetops and delineate crown boundaries. The proposed method was tested in natural forest stands with high canopy densities in Liangshui National Nature Reserve and Maoershan Forest Farm, Heilongjiang Province, China. Its performance was evaluated by an accuracy procedure that considered both the relative position of treetops and overlapped area of crowns. The average overall accuracy achieved was 85.12% for coniferous plots, 83.86% for deciduous plots and 86.44% for coniferous and broad-leaved mixed forest plots. The results revealed that the RHCSA method can detect and delineate individual tree crowns with little influence from forest types and crown size. It could provide technical support for individual tree crown delineation in coniferous, deciduous and mixed forests with high canopy densities.


Introduction
Forests are an important natural resource which play a significant role in modulating stores and fluxes of water and carbon, maintaining ecological diversity, and regulating climate on the Earth's surface, and provides timber and other forest products constantly which are closely related to humans [1][2][3].Parameters of individual trees, such as tree location, tree species, tree height, diameter at breast height (DBH), and crown diameter, are vital for sustainable and precise forest management [4][5][6].These individual tree-related properties can also be utilized to estimate forest parameters at stand level, such as tree species composition, mean tree height, timber volume, canopy density, and mean basal area [7][8][9][10].
Traditionally, forest parameters at the individual tree level have been measured by means of field surveys [11], which are time consuming and cost-intensive when carried out over broad areas [12,13].Since the early 1960s, with the extensive application of remote sensing data in forest inventory [14][15][16] a variety of image processing techniques were developed for automatically detecting and delineating individual tree crowns, such as local maxima filtering with fixed or variable window size [17][18][19], valley-following [20][21][22], region-growing [23,24], and watershed segmentation [25,26].These methods make it possible to obtain individual tree-based attributes that can be directly used as input parameters for environmental modeling without the limitations of the sample sizes and inaccessible areas [27,28].In recent decades, airborne Light Detection and Ranging (LiDAR), also referred to airborne laser scanning (ALS), is becoming one of the most common remotely-sensed data sources for forest inventory analysis [29][30][31].LiDAR remote sensing can provide detailed and precise three-dimensional structural information of the forest area [32], thus improving the development of individual tree crowns delineation (ITCD) algorithms [4,33].Laser-based individual tree crown delineation can reflect accurate geometrical properties of trees, such as tree height and crown diameter which cannot be affected by illumination angle and shadows on multispectral imageries [34], and become attractive techniques for both forestry and remote sensing communities [13].
According to the inputs of ITCD algorithms, LiDAR data can be mainly divided into three categories: LiDAR point cloud, LiDAR-derived raster, e.g., canopy height model (CHM), and raster-point combined data [13].With the increased point density of LiDAR data, a growing number of ITCD studies have begun to use LiDAR point clouds directly as the data source.Point cloud-based methods take advantage of multi-return information to describe vertical characteristics under canopies [35], which cause low information loss and provide an opportunity to multi-layer forests structure modelling [36,37].For example, Ferraz et al. [38] firstly defined forest layers and then applied a clustering method based on mean-shift algorithm using LiDAR point cloud.Vega et al. [39] developed a multi-scale dynamic segmentation method called PTrees to extract individual trees from point clouds.However, the method that focused on full exploitation of 3-D points may be currently restricted by computational burden, which is difficult to be applied broadly [40,41].Since raster-based methods are easy to be exploited and improved by powerful knowledge-based approaches [6,13], it becomes the most commonly used method for individual tree crown delineation.CHM, a type of LiDAR-derived raster, is the dominant dataset for representing canopy surface in individual tree crown delineation [8].Besides CHM, other LiDAR-derived datasets, such as Height-Scaled Crown Openness Index (HSCOI) and Canopy Maximum Model, were also used for ITCD [42,43].Sometimes, researchers combined both raster and LiDAR point data for ITCD studies.For example, researchers applied CHM for the initial segmentation and subsequently refined the results using a LiDAR point cloud [4,44,45].
Recently, many ITCD methods have been developed and applied for various forests conditions [43,46,47].Thereinto, most studies focus on coniferous forest because (1) many study area were located in high latitude regions that are dominated by coniferous forests; and (2) the basic assumption of a conical crown shape which many algorithms were developed based on is more appropriate for conifers than other forest types [11].However, it is not easy to obtain satisfied accuracy under other forest conditions.Ke and Quackenbush [34] compared watershed segmentation, region-growing, and valley-following method for different forest types on different images and found that all methods works well on coniferous stand and have demonstrated lower accuracy for the deciduous stand on either image set.Other researchers also found that the results of individual tree crown detection are largely influenced by forest types [48,49].Over-segmentation and under-segmentation may be greatly increased in complex and closed forests.
Majority of previous methods extracted individual tree crowns by two steps: (1) treetop detection; and (2) crown boundary delineation.Most treetop detection methods assumed that a treetop is the local maximum brightness or height value and detected individual treetops using local maxima filter.However, during treetop detection, branches, or sub-crowns of a deciduous tree may be detected as small trees (commission error), and small trees with smaller crown size or tree height than average trees tended to be omitted (omission error).To avoid issues mentioned above, researchers have developed a series of local maxima filters with variable window sizes to detect individual treetops [16,19,43].Although local maximum filtering with variable window sizes has been widely developed and used, structure of tree crowns are not always regular and cannot be represented by a simple 2-D filtering window.In the crown delineation process, watershed segmentation and region-growing, as two of classic image segmentation methods, were often employed based on the markers produced by treetop detection [16,43].Although, researchers have improved delineation algorithms by priori knowledge and morphology of crowns [46,50], the two-steps methods could transfer errors from treetop detection to crown boundary delineation procedure.
Researchers also developed some ITCD algorithms without prior treetop detection.For example, Jing et al. [51] segmented multispectral images at different scales by watershed algorithm and then integrated the multi-layered segments.This approach effectively reduced the over-segmentation which are caused by branches of the deciduous tree, but required the dominant tree crown sizes beforehand.Liu et al. [52] developed a boundary refining program called Fishing Net Dragging (FiND) based on watershed segmentation and used a boundary classification model to merge the adjacent segments.These approaches do not rely on prior detected treetops, but require some complicated post-process procedures.
Beyond the ITCD algorithms, accuracy assessment for ITCD studies also varies greatly, which make it difficult to compare the performance of the algorithms [35].Some studies evaluated ITCD accuracy at both plot and individual tree level (e.g., [16,34]), while others only provided accuracy assessment at individual tree level (e.g., [52]).From two-step perspective, ITCD accuracy assessment can be divided into detection accuracy and delineation accuracy assessment.Pouliot et al. [53] proposed accuracy index (AI) to represent detection accuracy and used the root mean square error (RMSE) of the crown diameters between reference and delineated crowns for delineation accuracy.Ke et al. [11] used a confusion table to describe the detail of tree crown detection and evaluated accuracy from both producer's and user's perspectives.Some studies used overall accuracy that was evaluated by the number of matched reference and detected trees [4,51].While some researchers analyzed different cases of individual tree crowns matching [51,52].It is necessary to develop a comprehensive and stable accuracy assessment system which contains individual tree crown matching and other accuracy indices for ITCD studies.
In this study, we aimed to propose a novel ITCD algorithm for different forest types in northeastern China based on LiDAR data.The specific objectives were to (1) develop a one-step ITCD algorithm that combines the procedures of treetop detection and crown boundary delineation into an iterative process, which could utilize horizontal relationship among crowns in vertical direction to reduce commission and omission errors; (2) explore the proposed algorithm under different forest types with high canopy density, including coniferous forest, deciduous forest, and conifer-broadleaves forest, and compare it with the marker-controlled watershed algorithm; (3) combine beneficial and complementary aspects of existed accuracy assessment methods into a new ITCD accuracy assessment program to evaluate the accuracy of proposed algorithm.It could provide technical support for individual tree crown delineation in coniferous, deciduous, and mixed forests with high canopy density.

Study Areas
In this study, nine plots in two study areas (study area I and II) were selected according to forest types.Study area I is primary forests dominated by conifers and study area II is secondary forests dominated by deciduous forests and coniferous and broad-leaved mixed forests.Both of study areas are the typical forests in Northeastern China.All plots were set to 100 × 100 m.

LiDAR Data
The LiDAR data were acquired on September 2009 and 2015 using a LiteMapper 5600 airborne LiDAR system (Integrated Geospatial Innovations company, Kreuztal, Germany) and CAF's LiCHy airborne observation system (The Chinese Academy of Forestry, Beijing, China) for study area I and study area II, respectively [55,56].These two sets of LiDAR data were both acquired by the Chinese Academy of Forestry (Beijing, China) and Northeast Forestry University (Harbin, Heilongjiang, China).The characteristics of LiDAR data for these two study areas are summarized in Table 1.

LiDAR Data
The LiDAR data were acquired on September 2009 and 2015 using a LiteMapper 5600 airborne LiDAR system (Integrated Geospatial Innovations company, Kreuztal, Germany) and CAF's LiCHy airborne observation system (The Chinese Academy of Forestry, Beijing, China) for study area I and study area II, respectively [55,56].These two sets of LiDAR data were both acquired by the Chinese Academy of Forestry (Beijing, China) and Northeast Forestry University (Harbin, Heilongjiang, China).The characteristics of LiDAR data for these two study areas are summarized in Table 1.

Reference Data
Since it is difficult to obtain accurate field data for validation, the reference data used in this study were visually interpreted from high spatial resolution orthoimages for accuracy assessment, which was a common method used in ITCD studies (e.g., [12,16,34]).The spatial resolution of orthoimageries were 0.5 m and 0.2 m for study area I and study area II, respectively.Both of the orthoimageries had three broad spectral bands: blue (with a center wavelength of 450 nm), green (550 nm), and red (625 nm).The reference data, including treetops and tree crowns, were determined using ArcGIS 10.2.

Data Preprocessing
Data preprocessing in our study consisted of three steps: (1) classify the raw LiDAR point data into ground and above-ground points using TerraSolid's TerraScan software (Terrasolid Ltd., Helsinki, Finland); (2) derive CHM using the classified point clouds, that is the difference between Digital Surface Model (DSM) and Digital Terrain Model (DTM) interpolated from surface points and ground points, respectively; (3) correct and smooth the CHM by using a pit-filling algorithm [57] and Gaussian filtering with a window size of 3 × 3 pixels.
In original CHM, some unnatural black or grey holes caused by abnormal or sudden changes of the height values exists regularly [58].These pits may cause potential errors in the treetop detection and crown delineation [3,59].In this study, the original CHM was filled by a semi-automated pit-filling algorithm outlined by Ben-Arie et al. [57].Then, a Gaussian filtering with window size of 3 × 3 pixels was used for image smoothing and noise removal according to trial-and-error.

Algorithm Description
In this study, we proposed a novel method called region-based hierarchical cross-section analysis (RHCSA) based on a canopy height model (CHM) derived from airborne LiDAR data for ITCD.Like other tree crown delineation algorithms (such as watershed segment, level-set, and valley-following), this algorithm considered tree canopies as a three-dimensional topographic surface that a tree crown can be manifested as a mountain-like uplift on a CHM.A tree apex is equal to a mountaintop and the height values decrease continuously from the treetop to crown the boundary.From a three-dimensional perspective, the RHCSA algorithm vertically scanned the CHM from top to bottom by a horizontal plane, see Figure 2a-c.In the process, the higher individual tree (tree B) produced a cross-section earlier (Figure 2d) than the shorter tree (tree A). Figure 2e,f indicated that the two cross-sections dilated with the decrease of scanning height.The cross-sections of tree A and B appeared at previous step must be contained in the subsequent cross-section.Cross-sections of individual tree crown regions that do not contact with others typically appear circular, as shown in Figure 2d,e, and the cross-sections produced by multiple contacted trees often appear irregular shape, like Figure 2f.Thus, RHCSA is a discretization method that slices the CHM by a series of equidistant horizontal planes to simulate stereoscopic scanning from top to bottom.After level cutting, a CHM was resolved into horizontal crown regions at different heights in vertical space.The nth cross-section region (CSR) in the ith layer of level cutting was defined as CSR n i , and the set containing all the CSR n i in ith level was defined as a CSR i .The containment relationship of CSR i among successive cutting levels can be defined as Equation ( 1) according to the vertical structure of raster data: Based on the containment relationship, RHCSA automatically determined whether a CSR n i was emerged for the first time, or produced by a single crown or an intersection of several crowns, then segmented the CSR n i that contained more than one tree from top to bottom.In RHCSA algorithm, each level cutting represented one iteration.Individual tree crowns and treetops will be extracted completely until all iteration ends (until level cutting reach the last layer).
Remote Sens. 2017, 9, 1084 6 of 21 segmented the CSR i n that contained more than one tree from top to bottom.In RHCSA algorithm, each level cutting represented one iteration.Individual tree crowns and treetops will be extracted completely until all iteration ends (until level cutting reach the last layer).

Iterative Process
In the iterative process of RHCSA algorithm, the iteration time is defined by users which can be calculated as follow: where I is the iteration time, H starting and H ending represent starting height and ending height in scanning process, respectively.H starting is equal to the maxima pixel-value of the input CHM, and H ending is a user-defined parameter which could be decided by the crown base height of study area.
H step is the step of level cutting.The detailed of vertical crown scanning may be improved with the step subdivision, but the efficiency of this algorithm will decrease.In this study, we set up a small step of 0.1 m to make sure that no obvious uplift exists during scanning.H ending was set to 2 m in this study in order to remove soil and grass.The iteration process of this method is illustrated in Figure 3.The main steps of one iteration are described as follows, including level cutting, identification of first emerged individual tree regions, classification of cross-section regions, segmentation of fusion regions, and adjustment of segments and markers.

Iterative Process
In the iterative process of RHCSA algorithm, the iteration time is defined by users which can be calculated as follow: where I is the iteration time, H starting and H ending represent starting height and ending height in scanning process, respectively.H starting is equal to the maxima pixel-value of the input CHM, and H ending is a user-defined parameter which could be decided by the crown base height of study area.
H step is the step of level cutting.The detailed of vertical crown scanning may be improved with the step subdivision, but the efficiency of this algorithm will decrease.In this study, we set up a small step of 0.1 m to make sure that no obvious uplift exists during scanning.H ending was set to 2 m in this study in order to remove soil and grass.The iteration process of this method is illustrated in Figure 3.
The main steps of one iteration are described as follows, including level cutting, identification of first emerged individual tree regions, classification of cross-section regions, segmentation of fusion regions, and adjustment of segments and markers.

 Level cutting
The iterative process started with level cutting which can be seen as an image binarization procedure at different levels.Cross-section image (CSI) produced by image binarization can be defined as Equation ( 3): where p means the pixel values on CHM and represents the threshold of image binarization at the i th level which can be calculated by Equation (4).On the CSI, crown regions have values of ones and non-crown regions have values of zeros.Then, cross-section regions were separated from CSI by connected-component labeling which is an algorithmic application of graph theory used in computer vision to detect connected regions in binary digital images [60].

Identification of the first emerged individual tree regions
During top-to-bottom level cutting, an individual tree crown follows a "growing" process among levels, that is, CSR of an individual tree emerges firstly and then grows with the iterations progressing.In one iteration, we defined the CSR i n that do not contain any CSR at the previous iteration ( CSR i-1 ) as the first emerged individual tree region (FEITR).Then, the centroid of the CSR i n belonging to FEITR was defined as a marker.An overall dataset called M FEITR was created to save all markers of FEITR produced at each level.M FEITR is defined as Equation ( 5): where centroid() is the function to determine the center of the current region, CSR i n is the nth crosssection region in the ith layer of level cutting.∩ represents intersection, and ∅ represents an empty set.

 Classification of cross-section regions
In this step, all cross-section regions were classified into two categories: individual tree region (ITR) and fusion region (FR).ITR is the cross-section region that belongs to an individual tree, and FR represents a cross-section of multiple contacted trees.According to the relationship of between CSR i n and M FEITR , CSR i n can be classified by Equation ( 6): • Level cutting The iterative process started with level cutting which can be seen as an image binarization procedure at different levels.Cross-section image (CSI) produced by image binarization can be defined as Equation ( 3): where p means the pixel values on CHM and T i represents the threshold of image binarization at the ith level which can be calculated by Equation (4).On the CSI, crown regions have values of ones and non-crown regions have values of zeros.Then, cross-section regions were separated from CSI by connected-component labeling which is an algorithmic application of graph theory used in computer vision to detect connected regions in binary digital images [60].

• Identification of the first emerged individual tree regions
During top-to-bottom level cutting, an individual tree crown follows a "growing" process among levels, that is, CSR of an individual tree emerges firstly and then grows with the iterations progressing.In one iteration, we defined the CSR n i that do not contain any CSR at the previous iteration (CSR i−1 ) as the first emerged individual tree region (FEITR).Then, the centroid of the CSR n i belonging to FEITR was defined as a marker.An overall dataset called M FEITR was created to save all markers of FEITR produced at each level.M FEITR is defined as Equation ( 5): where centroid() is the function to determine the center of the current region, CSR n i is the nth cross-section region in the ith layer of level cutting.∩ represents intersection, and ∅ represents an empty set.

• Classification of cross-section regions
In this step, all cross-section regions were classified into two categories: individual tree region (ITR) and fusion region (FR).ITR is the cross-section region that belongs to an individual tree, and FR represents a cross-section of multiple contacted trees.According to the relationship of between CSR n i and M FEITR , CSR n i can be classified by Equation ( 6): where count() is the function that determines count.ITR n i is the nth ITR in the ith level and FR n i is the nth FR in the ith level.When CSR n i has the only one intersection with M FEITR , it can be classified as ITR n i , otherwise it is classified as FR n i .

• Segmentation of fusion regions
After classification of cross-section regions, we employed marker-controlled watershed segmentation to segment FRs into separated ITRs.Watershed segmentation is one of the most popular techniques for segmentation and treats each concave tree crown in the inverted CHM as a catchment basin [61].We segmented FRs using M FEITR as markers.As a result, all FRs were divided into ITRs and there is no CSR with more than one marker inside.
However, branches and sub-crowns of deciduous trees may be misclassified as FEITR markers that are included in M FEITR and will cause commission errors in this step.Thus, we considered the containment relationship between M FEITR and CSR and the horizontal shapes of crown regions produced by level cutting to reduce commission errors.Similar to other algorithms [15,16,51], circularity was applied as an index to reduce over-segmentation, see Equation ( 7): where A is the area of the CSR n i and r is the largest distance between the centroid and the border of the CSR n i .At a certain height, an individual tree crown does not contract with other crowns, its cross-section usually appears circular.While a cross-section of a tree cluster is likely to be composed of overlapped circles, and the edge of the cluster also becomes somewhat lobate [51].Thus, simple rules have been added in the segmentation process to reduce commission errors.
where A and C are the area and circularity of the FR n i , respectively; A T and C T are thresholds of area and circularity, respectively.A T was designed for controlling crown area and avoiding the special condition that a CSR of a large closed tree cluster appears to be circular.In this study, A T was empirically set to 500 pixels.C T was designed to make sure that the shape of every crown is not far from circular.The parameter C T was set to 0.85 in this study as typical values mentioned in Wolf and Heipke [15].The sensitivity of parameters C T under different forest conditions will be discussed in Section 5.2.

• Adjustment of segments and markers
After classification and segmentation of the CSR, a pixel-based binary morphology opening operation was applied to refine the crown boundaries at each level and remove the objects that are smaller than the structuring element (SE) [51].We employed the opening operation with a disk SE width of three pixels after watershed segmentation to refine the segment boundaries and remove the irregular segment objects produced by small noise in the CHM.
During the segmentation and refinement procedure, the containment relationship between M FEITR and ITR may be changed.The ITR which was directly reclassify from FR may have more than one markers.The marker in M FEITR need to be adjusted before next iteration.Each individual tree crown should contain only one marker and only the highest one was preserved.Through this procedure, invalid markers (i.e., false treetops) caused by large branches and the noise in the CHM was removed efficiently.After the adjustment procedure, the RHCSA algorithm skips to the next iteration and continues level cutting the CHM.

Output of Individual Treetops and Crowns
According to above procedures, false treetops were removed from M FEITR , which represented the set of all emerged treetops at each level.After the iteration, M FEITR was the output of the detected treetops.ITR images represented the horizontal boundary of crowns at different heights.According to the containment relationship in the RHSCA algorithm, the ITR image at the last level contained all the previous ITRs and should be the output of delineated individual tree crowns.

Accuracy Assessment
Based on some existing accuracy assessment methods (e.g., [13,52]), an evaluation system was proposed to comprehensively assess the accuracy of ITCD studies.The accuracy assessment scheme started with individual tree crowns matching between delineated and reference crowns, then quantitative assessment metrics were calculated on the basis of matched crowns.

Individual Tree Crowns Matching
The accuracy assessment of individual tree crown delineation considers correspondence between delineated crowns and reference crowns on a location-by-location basis [11,13].Since individual tree crowns matching may be affected by the matching order, this study defined correspondence of "individuals" from the perspectives of reference crowns and detected crowns.The reference crown perspective represents how well each reference crown was delineated, while the detected crown perspective reflects whether a detected crown is represented by reference crowns.From each perspective, different correspondences were determined by considering both the overlapped crown area and the position of treetops.The accuracy metrics from the two perspectives were defined as Table 2.

Quantitative Assessment
In this study, producer's accuracy (PA) and user's accuracy (UA) were introduced from traditional accuracy assessment of classification [26,34], corresponding to the reference crown perspective and detected crown perspective, respectively.Producer's accuracy indicates the probability that a reference tree is correctly delineated, as defined in Equation ( 8).User's accuracy indicates the probability that a delineated tree will be correctly represented by reference trees, as defined in Equation ( 9): where N PM , N PN M are the number of 1:1 matches and near-matches, respectively, from reference crown perspective.N Re f is the total number of reference crowns.N UM , N UN M are the number of 1:1 matches and near-matches, respectively, from detected crown perspective.N Det is the total number of detected crowns.
Overall accuracy (OA) can describe the relationship between reference crowns and detected crowns on both perspectives.In this study, overall accuracy was calculated by the harmonic mean of PA and UA, as defined in Equation (10), which was similar to the definition of F-score in previous ITCD studies [35,39].
Since 1:1 match and near-match crowns from reference and detected crown perspectives are different, we introduced an "overall match" that consider the two cases from both perspectives.In this study, if a crown belongs to a 1:1 match or near-match crown from both detected and reference crown perspectives, it is defined as an "overall match".Detection accuracy (RMSE(P)) was defined as root mean square error (RSME) of the Euclidean distance between the reference and detected treetops of the overall match crowns [62], see (Equation ( 11)).Delineation accuracy (RMSE(D)) was defined as the difference of crown diameter between delineated and reference crowns [34], see (Equation ( 12)), where, dist() is the function of calculating Euclidean distance, N o_match is the total number of the overall match.Re f P i and DetP i represent the position of the reference and detected treetops for the ith overall crowns, respectively.Re f D i and DetD i represent the crown diameter of the reference and the delineated crowns for ith overall crown, respectively.

Individual Tree Crowns Matching
This study conducted the proposed RHCSA algorithm in nine plots of 100 × 100 m, which belonged to coniferous forest, coniferous-broadleaved mixed forest, and deciduous forest.A total of 2768 reference crowns were manually delineated according to high spatial resolution orthoimages.A total of 2898 crowns were delineated by the proposed algorithm and slightly over-segmented reference crowns.Since both the overlapped area and treetop locations were considered in the proposed accuracy assessment approach, the number of each cases from the reference crown perspective is different from that from detected crown perspective [52].Table 3 summarizes the result of individual tree crowns matching in the nine plots from both reference and detected crown perspectives.In the coniferous (Plots 1-3) and mixed forest plots (Plots 4-6), the numbers of reference crowns were always greater than that of detected crowns, while the results of deciduous stands are in direct contradiction.This is because deciduous trees with irregular branches tended to be over-segmented and conifers with small crowns tended to be omitted.In general, 1:1 matches are dominant in all metrics from both perspectives and the number of near-matches accounted for about 10% of the number of 1:1 matches.This illustrated that most reference crowns can be delineated well and a majority of delineated crowns can represent reference crowns.The number of "merge" errors (about 40) was much larger than that of "split" errors (about 10) in coniferous stands.In deciduous and mixed forest stands, "split" errors increased, especially in deciduous forests.However, "merge" errors of deciduous forests were much lower than that of the other two from both reference and detected crown perspectives.This is because a conifer tree that usually has a small crown tends to be merged with adjacent trees.However, deciduous trees that usually have larger crowns are not easy to be omitted.It is observed that no omission errors occurred except for Plot 3 (only two).This indicated that the RHCSA method has great potential to control omission and commission error.Due to the small number (less than five), "Multi-intersected" and "Mis-located match" were not the main errors in the nine plots, either.
Figures 4 and 5 illustrated the results of the detected individual treetops and tree crowns by the proposed methods for the nine plots.In these figures red points and blue points present detected treetops and reference treetops, respectively; red lines and blue lines represent the boundaries of detected individual tree crowns and reference individual tree crowns, respectively.The blue points that were not corresponding to red points in Figure 4 may belong to omission errors, split or multi-intersected errors.Overall, the proposed algorithm detected most individual treetops and successfully delineated most individual tree crowns in three types of forests.However, there are still some omitted trees and some over-segmented crowns in each plot.
Remote Sens. 2017, 9, 1084 11 of 21 usually has a small crown tends to be merged with adjacent trees.However, deciduous trees that usually have larger crowns are not easy to be omitted.It is observed that no omission errors occurred except for Plot 3 (only two).This indicated that the RHCSA method has great potential to control omission and commission error.Due to the small number (less than five), "Multi-intersected" and "Mis-located match" were not the main errors in the nine plots, either.Figures 4 and 5 illustrated the results of the detected individual treetops and tree crowns by the proposed methods for the nine plots.In these figures red points and blue points present detected treetops and reference treetops, respectively; red lines and blue lines represent the boundaries of detected individual tree crowns and reference individual tree crowns, respectively.The blue points that were not corresponding to red points in Figure 4 may belong to omission errors, split or multiintersected errors.Overall, the proposed algorithm detected most individual treetops and successfully delineated most individual tree crowns in three types of forests.However, there are still some omitted trees and some over-segmented crowns in each plot.

Quantitative Assessment
Quantitative assessment based on PA, UA, OA, detection accuracy, and delineation accuracy was conducted in this study.Table 4 shows that the RHCSA method had good performances and obtained PA, UA, and OA of around 85% for all nine plots.The accuracies of mixed forest plots were highest (i.e., PA: 87.02%, UA: 85.88%, OA: 86.44%), and followed by coniferous plots (i.e., PA: 85.35%, UA: 84.91%, OA: 85.12%).The lowest accuracies were obtained in deciduous plots (i.e., PA: 85.24%, UA: 82.55%, OA: 83.87%), but still above 80%.For the detection and delineation accuracy assessment, the detection accuracy (i.e., RMSE(P)) and delineation accuracy (i.e., RMSE(D)) varied from about 0.6-1 m and about 0.5-0.7 m for nine plots, respectively.The RHCSA method could accurately describe individual tree positions and crown diameters for overall matched trees.To further explore the performance of the RHCSA in three different forest types, ANOVA (Analysis of Variance) was conducted for the five accuracy metrics (i.e., PA, UA, OA, detection and delineation accuracy).Table 5 illustrated the means of the accuracy metrics in three forest types and the p-values of ANOVA.The null hypothesis of the ANOVA was that there was no statistically significant difference of every accuracy index among the three forest types.From Table 5, it showed that the p-values of all accuracy metrics were much greater than 0.05, which indicated that no statistically significant differences of accuracies were found among the three forest types.Thus, the proposed RHCSA algorithm had very stable performances for different forest types.

Comparison with Previous Studies
According to the above results analysis, it is evident that a region-based hierarchical cross-section analysis has obtained high and consistent accuracy for all three forest types.Under the natural forest condition with high canopy density, the RHCSA method was less affected by the variation of crown size and required less prior information.From Table 5, the highest overall accuracy (86.44%) was derived from the closed coniferous-broadleaved mixed forest, which is different from the other ITCD studies where the algorithms performed well in coniferous stands [46,52].That is because the RHCSA algorithm could detect "dumpy" protrusions on CHM and not be limited by the typical conical structures.Thus, RHCSA performed well in coniferous-broadleaved mixed forest and had stable performances for the three forest types.
For comparison, we implemented a local maxima filter (LMF) [63] and the marker-control watershed segmentation (MCWS) [43] to extract individual tree crowns in Plots 3, 6 and 9 (coniferous, mixed forest, and deciduous stands).Firstly, we adopted the same preprocessing steps as that we used in RHCSA.Secondly, 5 × 5 and 7 × 7 local maxima filters were used to detect treetops in Plots 3, 6 and 9, respectively, according to crown diameter.Then we employed a marker-control watershed segmentation using detected treetops as markers.The results are summarized in Tables 6 and 7.  From Table 6, according to the total number of reference crowns and detected crowns, the MCWS method generally overestimated the number of individual trees.For example, the MCWS method detected 461 individual crowns, compared to 361 reference crowns in Plot 6 The counts of "split" error of MCWS were quite high from both two perspectives, especially in deciduous and mixed forest plots.Compared to the results of RHCSA method in the same plot, the number of "split" error decreased by approximately 70%.It indicated that RHCSA method can effectively reduce the over-segmentation, which was a common error for the MCWS method.This is because the RHCSA method combines the procedures of treetop detection and crown delineation into an iterative process which considers the horizontal relationship between the crowns in vertical space, instead of considering an appropriate window size for a given tree crown.
From Table 7, the overall accuracies of the RHCSA method were consistently higher than those of the MCWS method, and obtained 85.36%, 84.46%, and 83% in mixed forest, coniferous and deciduous stands, respectively.Different from the RHCSA method, the highest overall accuracy (74.35%) of the MCWS method was acquired in coniferous plot (Plot 3), and followed by mixed forest (64.29%,Plot 6) and deciduous forest (57.24%,Plot 9).The MCWS method is ideal for detecting and delineating conifers with uniform crown sizes [11] and difficult to perform well in deciduous forests because of various sizes and shapes.For detection and delineation accuracy, there was no significant difference between the two methods.Both methods could provide accurate information on tree location and crown diameter for overall matched trees (less than 1 m for detection accuracy; around 0.5 m for delineation accuracy).
To better delineate individual tree crowns, many study utilize complicated user-defined parameters in ITCD algorithms according to field or empirical data [51,64].For example, Jing et al. [65] proposed multi-scale crown delineation method and correctly delineated 69%, 65%, and 73% of individual tree crowns in coniferous, deciduous, and mixed forests, respectively.This method had great advantages in eliminating false treetops caused by branches (commission error).However, the ITCD accuracy greatly depended on the selection of the scale.A few representative scales may not be adequate for complex forests with various trees species.In RHCSA algorithm, level cutting was applied to avoid over-segmentation caused by various crown sizes instead of scales, and no complicated parameters were required.The main parameter that directly affect the accuracy and efficiency in the RHCSA algorithm is level cutting step.The default step was set to 0.1 m which is adequate for most situations, because the horizontal shape of crowns varies little in every 0.1 m on the CHM.The efficiency of the ITCD algorithm depends on iteration times (level cutting step) and image size without the influence of the number of trees.It took about 60-70 s to complete ITCD for the each 1 ha plot, using Matlab R 2014a on a computer with a 2.5 GHz Intel Core i7-4710 HQ CPU.When we shorted the step into 0.01 m, the time consuming increased to about 600 s while the overall accuracy has an increase of only 2%.Readers could ask for the application code of the RHCSA algorithm for academic use.The sample data (Plots 3, 6 and 9) is available online for academic use (mat format).

Sensitivity Analysis of C T
Another major parameter that affected the overall accuracy in the RHCSA algorithm was the circularity, which was applied for segmenting fusion regions in the iterative process.To further investigate the impact of circularity under different forest conditions, we adjusted the circularity threshold (C T ) from 0.5 to 1 with a step of 0.5 and analyzed the variation of overall accuracy in all plots.Figure 6 showed the impact of the circularity threshold on the overall accuracy of ITCD for three forest types.For coniferous forests, the overall accuracy increased rapidly with the increase of C T until reaching a sill, see Figure 6a.A small C T value tended to prevent fusion regions from further segmentation and resulted in more "merge" errors.Figure 6a illustrated that the C T with the highest overall accuracy for coniferous forests was about 0.9, which indicated that the crown of a conifer tree appears close to a circle [25].Although the overall accuracy increased with the increase of C T for mixed and deciduous forests as well, the increase rates were not as high as that of coniferous forest.The C T that led to the peak values of the overall accuracy for mixed forest and deciduous forest were slightly lower than that for coniferous forest, that is, about 0.8-0.85 and 0.7-0.75,respectively (see Figure 6b,c).This is because that coniferous crowns are more circular than deciduous and mixed forests, while deciduous trees have larger and more irregular crowns.Although a higher C T could lead to higher accuracy in coniferous stands, the average circularity will decrease and a lower C T could obtain higher accuracy with the number of deciduous trees increased in a plot.
slightly lower than that for coniferous forest, that is, about 0.8-0.85 and 0.7-0.75,respectively (see Figure 6b,c).This is because that coniferous crowns are more circular than deciduous and mixed forests, while deciduous trees have larger and more irregular crowns.Although a higher C T could lead to higher accuracy in coniferous stands, the average circularity will decrease and a lower C T could obtain higher accuracy with the number of deciduous trees increased in a plot.

Limitations
Although RHCSA is a promising ITCD algorithm, several limitations still exist.First of all, RHCSA considers CHM as mountain-like topographic surfaces, some flat crowns and suppressed trees without a dominant protrusion on CHM are difficult to detect and delineate.It becomes serious when a small tree is under adjacent tree crowns and causes "merge" error.Second, since RHCSA was developed for the forests in Northeastern China, it performed well for boreal temperate forest.The performances of other forests, like bamboo forests and rainforests, are unknown due to lack of data.In the future, the RHCSA method could be validated using other types of forest data from other countries.
Additionally, the topographic effect is another factor that influenced the accuracy of the ITCD method [6,66].Since the RHCSA algorithm detection and delineation of individual tree crowns relied on the "peak" and "valley" structures shown on CHM, it could result in errors if the structure changed due to extreme topographic factors (e.g., steep slope).However, the terrain is relatively flat in our study area, and no significant correlation was found between the slope and the overall accuracy based on correlation analysis.From a topography perspective, the RHCSA algorithm performed well for a relatively flat area.The performance influenced by other topographic factors (e.g., slope) should be investigated in the future.

Conclusions
Individual tree crown delineation using airborne LiDAR techniques plays a significant role in precise forestry and forest inventory and analysis.In this study, we developed a region-based hierarchical cross-section analysis (RHCSA) for treetop detection and crown delineation using LiDAR-derived CHM.This novel method used level cutting method to detect individual trees and delineate its horizontal boundaries at each level through an iterative process.The proposed algorithm is a one-step individual tree crown delineation (ITCD) algorithm with a few user-defined parameters and not influenced by treetop detection accuracy.Additionally, we also developed a comprehensive accuracy assessment scheme which could provide a guidance for ITCD studies and future applications using single tree approach.

Limitations
Although RHCSA is a promising ITCD algorithm, several limitations still exist.First of all, RHCSA considers CHM as mountain-like topographic surfaces, some flat crowns and suppressed trees without a dominant protrusion on CHM are difficult to detect and delineate.It becomes serious when a small tree is under adjacent tree crowns and causes "merge" error.Second, since RHCSA was developed for the forests in Northeastern China, it performed well for boreal temperate forest.The performances of other forests, like bamboo forests and rainforests, are unknown due to lack of data.In the future, the RHCSA method could be validated using other types of forest data from other countries.
Additionally, the topographic effect is another factor that influenced the accuracy of the ITCD method [6,66].Since the RHCSA algorithm detection and delineation of individual tree crowns relied on the "peak" and "valley" structures shown on CHM, it could result in errors if the structure changed due to extreme topographic factors (e.g., steep slope).However, the terrain is relatively flat in our study area, and no significant correlation was found between the slope and the overall accuracy based on correlation analysis.From a topography perspective, the RHCSA algorithm performed well for a relatively flat area.The performance influenced by other topographic factors (e.g., slope) should be investigated in the future.

Conclusions
Individual tree crown delineation using airborne LiDAR techniques plays a significant role in precise forestry and forest inventory and analysis.In this study, we developed a region-based hierarchical cross-section analysis (RHCSA) for treetop detection and crown delineation using LiDAR-derived CHM.This novel method used level cutting method to detect individual trees and delineate its horizontal boundaries at each level through an iterative process.The proposed algorithm is a one-step individual tree crown delineation (ITCD) algorithm with a few user-defined parameters and not influenced by treetop detection accuracy.Additionally, we also developed a comprehensive accuracy assessment scheme which could provide a guidance for ITCD studies and future applications using single tree approach.
Results showed that RHCSA method obtained stable and high accuracy for different forest types, including coniferous forest, coniferous-broadleaves forest and deciduous forest.The accuracies of mixed forest plots were highest (i.e., producer's accuracy (PA): 87.02%, user's accuracy (UA): 85.88%, overall accuracy (OA): 86.44%), followed by coniferous plots (i.e., PA: 85.35%, UA: 84.91%, OA: 85.12%) and deciduous plots (i.e., PA: 85.24%, UA: 82.55%, OA: 83.87%).The detection and delineation accuracy (i.e., RMSE(P) and RMSE(D)) were relatively small (less than 1 m).RHCSA method effectively improved OA by more than 10% comparing with the traditional marker-control watershed segmentation (MCWS) algorithm, especially for deciduous forest (more than 20%) and mixed forest (about 20%).Both RHCSA and MCWS algorithms could provide accurate information on tree location and crown diameter for overall matched trees (less than 1 m for detection accuracy; around 0.5 m for delineation accuracy).However, there are still some limitations in RHCSA method.Since the RHCSA process only considers the segmentation of canopy surface, some intersected crowns and small crowns under the surface are difficult to detect and delineate.Next, RHCSA has been designed and tested for boreal temperate forest.It is still necessary to investigate the performances of other forests (e.g., bamboo forests and rainforests).In addition, although it is suitable to use the RHCSA algorithm for a relatively flat area, the effect of topographic factors (e.g., slope) on accuracy should be explored in the future.In future, ITCD studies could combine basic image segmentation algorithms and tree morphological characteristics to improve the accuracy of the crown delineation based on LiDAR point cloud data.

Figure 2 .
Figure 2. A 3D schematic diagram of RHCSA algorithm based on a CHM containing tree A and tree B. (a-c) represents vertical scanning process on CHM; (d-f) are the corresponding cross-sections produced from (a) to (c), respectively.

Figure 2 .
Figure 2. A 3D schematic diagram of RHCSA algorithm based on a CHM containing tree A and tree B. (a-c) represents vertical scanning process on CHM; (d-f) are the corresponding cross-sections produced from (a) to (c), respectively.

Figure 4 .
Figure 4. Detected treetops using the RHCSA algorithm for all nine plots, which C T was set to 0.85.The blue points represent reference treetops and the red points represent detected treetops.

Figure 4 .
Figure 4. Detected treetops using the RHCSA algorithm for all nine plots, which C T was set to 0.85.The blue points represent reference treetops and the red points represent detected treetops.

Figure 5 .
Figure 5. Detected tree crowns using the RHCSA algorithm for all nine plots, which C T was set to 0.85.The blue boundaries represent reference crowns and the red boundaries represent detected crowns.

Figure 5 .
Figure 5. Detected tree crowns using the RHCSA algorithm for all nine plots, which C T was set to 0.85.The blue boundaries represent reference crowns and the red boundaries represent detected crowns.

Figure 6 .
Figure 6.The impact of C T on the overall accuracy of ITCD in three forest types: (a) coniferous forest, (b) coniferous-broadleaved mixed forest, and (c) deciduous forest.

Figure 6 .
Figure 6.The impact of C T on the overall accuracy of ITCD in three forest types: (a) coniferous forest, (b) coniferous-broadleaved mixed forest, and (c) deciduous forest.

Table 1 .
Characteristics of LiDAR data used for two study areas.

Table 1 .
Characteristics of LiDAR data used for two study areas.

Table 2 .
Summary of accuracy metrics for individual tree crowns matching from reference crown and detected crown perspectives.

Table 3 .
Accuracy assessment of individual tree crowns matching on both reference and detected crown perspectives using RHCSA method.

Table 4 .
Quantitative assessment results of the RHCSA algorithm.

Table 5 .
Mean and ANOVA results of the five accuracy metrics for three forest types.

Table 6 .
Comparison of individual tree crowns matching between MCWS and RHCSA.RCP: reference crown perspective; DCP: detected crown perspective; 2 total: the total number of reference crowns and detected crowns on reference crown and detected crown perspective, respectively. 1

Table 7 .
Comparison of accuracy metrics between MCWS and RHCSA.