A Novel Approach to Characterizing Crown Vertical Proﬁle Shapes Using Terrestrial Laser Scanning (TLS)

: Crown vertical proﬁles (CVP) play an essential role in stand biomass and forest ﬁre prediction. Traditionally, due to measurement difﬁculties, CVP models developed based on a small number of individual trees are not convincing. Terrestrial laser scanning (TLS) provides new insights for researching trees’ CVPs. However, there is a limited understanding of the ability to accurately describe CVPs with TLS. In this study, we propose a new approach to automatically extract the crown radius (CR) at different heights and conﬁrm the correctness and effectiveness of the proposed approach with ﬁeld measurement data from 30 destructively harvested sample trees. We then applied the approach to extract the CR from 283 trees in 6 sample plots to develop a two-level nonlinear mixed-effects (NLME) model for the CVP. The results of the study showed that the average extraction accuracy of the CR when the proposed approach was applied was 90.12%, with differences in the extraction accuracies at different relative depths into the crown (RDINC) ranges. The TLS-based extracted CR strongly correlated with the ﬁeld-measured CR, with an R 2 of 0.93. Compared with the base model, the two-level NLME model has signiﬁcantly improved the prediction accuracy, with R a2 increasing by 13.8% and RMSE decreasing by 23.46%. All our research has demonstrated that TLS has great potential for accurately extracting CRs, which would provide a novel way to nondestructively measure the crown structure. Moreover, our research lays the foundation for the future development of CVP models using TLS at a regional scale.


Introduction
The crown vertical profile (CVP) of a tree is the minimum boundary surrounding the crown that connects the tree top and the branch vertices [1]. It reflects the shape and size of the crown and directly describes peoples' overall impression of a tree. The CVP has become one of the most critical forest structure parameters, both at the sample plot and tree levels, and it lays the foundation for exploring the spatial attributes of tree crowns, such as the crown volume and surface area [2,3]. It can reflect the horizontal distribution of the leaf quantity and biomass in the tree crown [4] and can also function as an indicator for forest fire prediction models [5][6][7]. Therefore, the accurate simulation of the trees' CVP holds great significance for forest management and fire prediction.
Numerous models have been developed to describe trees' CVPs [8]. The crown radius (CR) is an essential variable for developing CVP models; however, the challenge of accurately obtaining CR information remains unresolved. In the past, it was frequently necessary to destructively harvest trees first, and then measure the branch attributes, and finally, use trigonometric relationships to compute the crown radius [9,10]. This The study area is in Mengjiagang Forestry Farm, Huanan County, Jiamusi City, Heilongjiang Province, which was established in February 1956 and is one of the first stateowned forestry farms to be established in Heilongjiang Province after the founding of the People's Republic of China, subordinate to the Jiamusi Forestry Bureau. Moreover, the administrative area belongs to the territory of Huanan County. The geographical coordinates are 130 • 32 ~130 • 52 E, 46 • 20 ~46 • 30 N, as shown in Figure 1. The landforms are mainly low hills with gentle slopes and an average elevation of 250 m. The forestry operation covers a total area of approximately 16,000 hectares, with a total forest accumulation of 1.41 million cubic meters. The forestry operation is mainly based on planted forests, of which 2/3 of the entire operation area is planted forests, and 1/3 is natural secondary forests, with a forest coverage rate of 81.7%. The tree species are primarily conifers, including Pinus koraiensis, Picea asperata, Pinus sylvestris var. mongholica, Larix olgensis, etc.

Field Measurement Data
In September 2020, we investigated 6 planted Korean pine sample plots of 0.06 ha in size (20 m × 30 m) at the Mengjiagang forest farm. These samples had different attributes (density, age, and elevation). All individual trees with a DBH over 5 cm in the sample plots were measured. Specifically, the DBH of the tree was measured using diameter tape, the tree height and crown base height (HCB) were measured with an ultrasonic altimeter (Vertex IV, Haglöf Sweden, Sweden), and the crown width (CW) of the tree in different directions was measured with steel tape. In addition, we also measured the geographic coordinates of all trees in the sample plots by combining real-time kinetic (RTK) technology and a global navigation satellite system (GNSS). More detailed statistical information regarding the sample plots is shown in Table 1.

Field Measurement Data
In September 2020, we investigated 6 planted Korean pine sample plots of 0.06 h size (20 m × 30 m) at the Mengjiagang forest farm. These samples had different attrib (density, age, and elevation). All individual trees with a DBH over 5 cm in the sam plots were measured. Specifically, the DBH of the tree was measured using diameter t the tree height and crown base height (HCB) were measured with an ultrasonic altim (Vertex IV, Haglöf Sweden, Sweden), and the crown width (CW) of the tree in diffe directions was measured with steel tape. In addition, we also measured the geogra coordinates of all trees in the sample plots by combining real-time kinetic (RTK) tech ogy and a global navigation satellite system (GNSS). More detailed statistical informa regarding the sample plots is shown in Table 1. The sample trees were selected with the following method. Each plot's trees w arranged according to their DBH in ascending order. The cumulative percentage of basal area was then divided into five even intervals, each representing 20% of the basal area, after the cumulative basal area had been calculated and divided by the basal area. Finally, the average diameter of the trees in each interval was calculated, the tree with a value closest to the average diameter was selected as the sampling tre total of 30 sample trees were selected to destructively measure the branch attributes  The sample trees were selected with the following method. Each plot's trees were arranged according to their DBH in ascending order. The cumulative percentage of the basal area was then divided into five even intervals, each representing 20% of the plot basal area, after the cumulative basal area had been calculated and divided by the plot basal area. Finally, the average diameter of the trees in each interval was calculated, and the tree with a value closest to the average diameter was selected as the sampling tree. A total of 30 sample trees were selected to destructively measure the branch attributes. We took care to cut down sample trees and tried to keep the crown undamaged. After felling, the tree height and HCB were measured using a tape measure, and the distance from the HCB to the tree top was the crown length (CL). Then, the whole tree was cut into 1 m sections at the beginning of the stump, and the sections were numbered. When measuring the branching factor, stem sections were placed vertically on the ground to maintain the branches in their most natural growth state and group all the branches with a vertical spacing of less than 20 cm into the same whorl. The following branching factors ( Figure 2) were measured: the branch length (BL), which is the length measured tightly against the branch; the branch chord length (BCL), which is the straight line distance between the start and end of the branch; the branch angle (BA), which is the angle between the chord length and the trunk of the tree; the branch depth into the crown (DINC). Additional crown parameters were also calculated: the crown radius (CR), CR = BCL × sin (BA) and the relative depth into the crown (RDINC), RDINC = (DINC -BCL × cos (BA))/CL. A total of 2581 branches were measured from 30 sample trees. We selected the largest branch (this branch has the largest CR) in each whorl of all sample trees, 616 in total, and used them as reference data to validate the extraction of the CR using point cloud data. Descriptive statistics of the sample trees and branches are shown in Table 2.
statistics of the sample trees and branches are shown in Table 2.  Figure 2. Illustration of tree factors and branch factors. CW is the crown width, CL is the crown length, DINC is the depth into the crown of the branch, HCB is the crown base height, HT is the tree height, BL is the branch length, BCL is the branch chord length, BA is the branch angle, and CR is the crown radius.

TLS Data
On a clear and windless day, we collected point cloud data from 6 sample plots and 30 sample trees via Trimble TX8 (Trimble, CA, USA) terrestrial laser scanning (Table A1). Regarding the sample plots, we scanned sample plots from five positions (Figure 3, right) to ensure the integrity of the point clouds of the sample trees within the plots and to avoid missing point clouds due to mutual occlusion between the sample trees [63]. Specifically, Figure 2. Illustration of tree factors and branch factors. CW is the crown width, CL is the crown length, DINC is the depth into the crown of the branch, HCB is the crown base height, HT is the tree height, BL is the branch length, BCL is the branch chord length, BA is the branch angle, and CR is the crown radius.

TLS Data
On a clear and windless day, we collected point cloud data from 6 sample plots and 30 sample trees via Trimble TX8 (Trimble, CA, USA) terrestrial laser scanning (Table A1). Regarding the sample plots, we scanned sample plots from five positions (Figure 3 right) to ensure the integrity of the point clouds of the sample trees within the plots and to avoid missing point clouds due to mutual occlusion between the sample trees [63]. Specifically, a scanning station was placed in the center of the sample plot, and the scanning time was set to 10 min. Four corners of the sample plot were each placed in a scanning station, and the scanning time was set to 3 min. Notably, when the scanning stations is placed at the sample plot corner points, the TLS should be placed at a distance of 5 m outside the sample plot corner points to ensure the integrity of the sample trees near the sample plot corner points. Regarding the purpose of the co-registration of multi-station cloud data, 10 reflection target balls (RETB) were placed evenly around the sample plot, with the height of the RETB was staggered as much as possible to avoid mutual blockage between the RETB and to ensure that at least three or more RETB at the same time can be observed at each scanning site. The average time spent collecting TLS data for each sample plot was 45 min. the sample tree at an angle of 120° as far away as possible, and the distance between each scanning station and the sample tree was approximately the size of the tree's height. Five RETBs were evenly placed within 2 m of the sample trees, and at least two RETBs could be observed at each scanning station to allow later co-registration of the point cloud data. Weeds and low bushes that may cause shading around the sample trees were removed before scanning began, and the scanning time was set to 10 min for each scanning station.

Pre-Processing the Point Cloud Data
The co-registration of multi-station point cloud data was implemented using Trimble RealWorks 11.2 software (Trimble, CA, USA), which accompanied Trimble TX8. A singlepoint matching method was used for co-registration using RETB as the reference target. When the error was less than 5 mm after refinement during the matching process, it was considered to be successful. To solve the problem of data redundancy, we used 'Octree' to downsample the point cloud data [64,65]. An improved progressive triangulated irregular network densification (IPTD) filtering algorithm [66,67] was used to separate ground points from non-ground points. The digital elevation model (DEM) was then generated via irregular triangular mesh interpolation to facilitate the normalization of the point cloud. TLS works via a bottom-up scanning method that clearly identifies the stem; therefore, we apply a bottom-up approach for individual tree segmentation. The trunk was first For sample trees, we scanned sample trees from three positions (Figure 3 left) to ensure that the point cloud could cover the whole tree completely, but not cause redundancy in the data [53]. The scanning stations were evenly distributed around the perimeter of the sample tree at an angle of 120 • as far away as possible, and the distance between each scanning station and the sample tree was approximately the size of the tree's height. Five RETBs were evenly placed within 2 m of the sample trees, and at least two RETBs could be observed at each scanning station to allow later co-registration of the point cloud data. Weeds and low bushes that may cause shading around the sample trees were removed before scanning began, and the scanning time was set to 10 min for each scanning station.

Pre-Processing the Point Cloud Data
The co-registration of multi-station point cloud data was implemented using Trimble RealWorks 11.2 software (Trimble, CA, USA), which accompanied Trimble TX8. A singlepoint matching method was used for co-registration using RETB as the reference target. When the error was less than 5 mm after refinement during the matching process, it was considered to be successful. To solve the problem of data redundancy, we used 'Octree' to downsample the point cloud data [64,65]. An improved progressive triangulated irregular network densification (IPTD) filtering algorithm [66,67] was used to separate ground points from non-ground points. The digital elevation model (DEM) was then generated via irregular triangular mesh interpolation to facilitate the normalization of the point cloud. TLS works via a bottom-up scanning method that clearly identifies the stem; therefore, we apply a bottom-up approach for individual tree segmentation. The trunk was first detected using a clustering algorithm, and then the shortest path algorithm based on metabolic ecology theory was used to segment the crown [68] for the purpose of segmenting the complete individual tree. The individual tree segmentation step was performed using Lidar360 software (Beijing Digital Green Soil Technology Co., Ltd., Beijing, China), and the segmentation results contain DBH, CW and tree height information, as well as the coordinate positions of the individual trees. Figure 4 summarizes the workflow of this study. detected using a clustering algorithm, and then the shortest path algorithm based on metabolic ecology theory was used to segment the crown [68] for the purpose of segmenting the complete individual tree. The individual tree segmentation step was performed using Lidar360 software (Beijing Digital Green Soil Technology Co., Ltd., Beijing, China), and the segmentation results contain DBH, CW and tree height information, as well as the coordinate positions of the individual trees. Figure 4 summarizes the workflow of this study.

Automatic Detection of the HCB
The crown base height (HCB) is generally defined as the vertical distance from the ground at the base of a continuous living crown [69][70][71][72], and the precision of the HCB determines the integrity of the crown. In this study, the principle of automatic HCB detection is based on variation in the point cloud frequency distribution over the vertical profile of individual trees [1,11,73]. When the point cloud frequency increases dramatically along the tree height direction for the first time and is above a certain threshold, we determine this position as the HCB ( Figure 5). To obtain details, first, the individual tree  The crown base height (HCB) is generally defined as the vertical distance from the ground at the base of a continuous living crown [69][70][71][72], and the precision of the HCB determines the integrity of the crown. In this study, the principle of automatic HCB detection is based on variation in the point cloud frequency distribution over the vertical profile of individual trees [1,11,73]. When the point cloud frequency increases dramatically along the tree height direction for the first time and is above a certain threshold, we determine this position as the HCB ( Figure 5). To obtain details, first, the individual tree point clouds were sliced from bottom to top at a thickness of 0.1 m, and then the number of point clouds within each slice was counted to calculate the percentage of the number of point clouds within each slice to the total number of point clouds in the individual tree. Afterward, the vertical distribution of the point cloud frequency was smoothed using Gaussian regression, followed by the calculation of the local maximum and minimum values. The HCB is defined as the local minimum before the first local maximum above the mean point cloud frequency. To verify the effectiveness of automatic HCB detection, we calculated the detection accuracy of the HCB with the field-measured HCB as the reference value. point clouds were sliced from bottom to top at a thickness of 0.1 m, and then the num of point clouds within each slice was counted to calculate the percentage of the numbe point clouds within each slice to the total number of point clouds in the individual t Afterward, the vertical distribution of the point cloud frequency was smoothed us Gaussian regression, followed by the calculation of the local maximum and minimum ues. The HCB is defined as the local minimum before the first local maximum above mean point cloud frequency. To verify the effectiveness of automatic HCB detection, calculated the detection accuracy of the HCB with the field-measured HCB as the re ence value. is an il tration depicting the process of automatically detecting HCB. The gray points are the point cl frequency within each slice. The green line is the smoothed point cloud frequency curve. The and blue points represent the local minima and local maxima detected on the curve, respectiv The vertical dotted line represents the mean point cloud frequency and the horizontal black represents the detected HCB (i.e., the local minimum before the first local maximum above the m point cloud frequency).

Automatic Detection of the Tree Trunks
Regarding field measurements, the horizontal distance from the branch tip to trunk surface is commonly regarded as the CR at a certain crown level. It is essentia detect trunk points from the tree crown point cloud data to more accurately compare analyze the differences between the CR extracted using point cloud data and the ac CR measured in the field. In this study, we applied the random sample consensus (RA SAC) algorithm [74,75] for the detection of trunk point clouds in a 2D projection pla RANSAC is an iterative algorithm that correctly estimates the mathematical mo is an illustration depicting the process of automatically detecting HCB. The gray points are the point cloud frequency within each slice. The green line is the smoothed point cloud frequency curve. The red and blue points represent the local minima and local maxima detected on the curve, respectively. The vertical dotted line represents the mean point cloud frequency and the horizontal black line represents the detected HCB (i.e., the local minimum before the first local maximum above the mean point cloud frequency).

Automatic Detection of the Tree Trunks
Regarding field measurements, the horizontal distance from the branch tip to the trunk surface is commonly regarded as the CR at a certain crown level. It is essential to detect trunk points from the tree crown point cloud data to more accurately compare and analyze the differences between the CR extracted using point cloud data and the actual CR measured in the field. In this study, we applied the random sample consensus (RANSAC) algorithm [74,75] for the detection of trunk point clouds in a 2D projection plane. RANSAC is an iterative algorithm that correctly estimates the mathematical model parameters from a set that contains "outliers" [76]. A basic assumption for the RANSAC algorithm is that the data are composed of "inliers" and "outliers". The "inliers" are data that contribute to the model parameters, and the "outliers" are data that do not fit the model.
A tree trunk's projection plane is generally considered to be circular or elliptical. These shapes can be perfectly adapted to the detection of trunk point clouds by iteratively fitting circles using the RANSAC algorithm. In this study, the trunk point cloud was sliced with a thickness of 50 cm to avoid projection plane distortion caused by excessive trunk tilt, and then RANSAC was applied to detect the trunk point cloud in the two-dimensional Remote Sens. 2023, 15, 3272 9 of 25 projection of the sliced point cloud. RANSAC is, by nature, an iterative algorithm [77,78]. Thus, detecting the trunk point cloud from the sliced crown point cloud requires a large amount of calculations and may cause the fit to fail since the crown point density ('outlier') is much greater than the trunk point density is ('inliers'). Therefore, to improve the operation efficiency and success rate of the algorithm, we continued to refine the sliced point clouds to ensure that the trunk point clouds were dominant in the refined sliced point clouds. Specifically, we use a 50 cm thick point cloud slice below the HCB height as the initial slice from which the initial circle center and diameter were detected. When we were detecting a new slice, we constructed a buffer ( Figure 6) with the fitted circle center and diameter from the previous slice as the new circle center and search radius, respectively, to refine the point cloud in the new slice.
that contribute to the model parameters, and the "outliers" are data that do not fi model.
A tree trunk's projection plane is generally considered to be circular or ellip These shapes can be perfectly adapted to the detection of trunk point clouds by itera fitting circles using the RANSAC algorithm. In this study, the trunk point cloud was with a thickness of 50 cm to avoid projection plane distortion caused by excessive tilt, and then RANSAC was applied to detect the trunk point cloud in the two-dimens projection of the sliced point cloud. RANSAC is, by nature, an iterative algorithm [7 Thus, detecting the trunk point cloud from the sliced crown point cloud requires a amount of calculations and may cause the fit to fail since the crown point density ('ou is much greater than the trunk point density is ('inliers'). Therefore, to improve the ation efficiency and success rate of the algorithm, we continued to refine the sliced clouds to ensure that the trunk point clouds were dominant in the refined sliced clouds. Specifically, we use a 50 cm thick point cloud slice below the HCB height a initial slice from which the initial circle center and diameter were detected. When we detecting a new slice, we constructed a buffer ( Figure 6) with the fitted circle cente diameter from the previous slice as the new circle center and search radius, respect to refine the point cloud in the new slice.

Automatic Detection of the Crown Outer Contour Points
The CVP model is a function of the CR at any height as the dependent variabl the distance from the CR to the treetop as the independent variable, and it is a curvi equation used to model the pattern of vertical variation in the shape of the tree crow Therefore, it is essential to obtain the CR at different heights to simulate the CVP. I study, we separated the crown point cloud (Figure 7a) from the individual tree point according to the HCB in Section 2.5.1, and then used the point cloud sliced proje method for the extraction of the crown radius using MATLAB R2020b (MathWorks tick, MA, USA) software. With reference to Gao [61] and Quan [1] et al., we chose 0.5 the slice thickness (Figure 7b). The convex hull algorithm [79] is a common algorithm for finding the outermost points from a scattered set of points, and many studies been conducted and used this algorithm for the extraction of tree crown attributes 84]. We use the 2D convex hull algorithm to find the outer contour point from the c

Automatic Detection of the Crown Outer Contour Points
The CVP model is a function of the CR at any height as the dependent variable and the distance from the CR to the treetop as the independent variable, and it is a curvilinear equation used to model the pattern of vertical variation in the shape of the tree crown [1]. Therefore, it is essential to obtain the CR at different heights to simulate the CVP. In this study, we separated the crown point cloud (Figure 7a) from the individual tree point cloud according to the HCB in Section 2.5.1, and then used the point cloud sliced projection method for the extraction of the crown radius using MATLAB R2020b (MathWorks, Natick, MA, USA) software. With reference to Gao [61] and Quan [1] et al., we chose 0.5 m as the slice thickness (Figure 7b). The convex hull algorithm [79] is a common algorithm used for finding the outermost points from a scattered set of points, and many studies have been conducted and used this algorithm for the extraction of tree crown attributes [80][81][82][83][84]. We use the 2D convex hull algorithm to find the outer contour point from the crown point cloud after projection (Figure 7c) and subsequently define the difference between the distance from the outer contour point to the trunk centroid and the trunk radius of that layer as the CR. The trunk detection method is described in Section 2.5.2.
TLS takes a bottom-up scanning approach; so, some 'noise' will inevitably appear in the upper part of the tree crown, resulting in an abnormally extracted CR. In this study, we applied the 'boxplot method' [85] to remove outliers from the extracted CR. After removing the outliers from the CR, the largest radius was selected as the maximum CR for that layer. The principle of the 'boxplot method' to detect outliers is to label the points above the upper limit and below the lower limit as outliers ( Figure 8). The formulas for calculating the upper and lower limits are as follows: where Q3 is the third quartile, which is also known as the "upper quartile" and is equal to the 75th percentile of all values in the sample ordered in ascending order; Q1 is the first quartile, which is also known as the "lower quartile" and is equal to the 25th percentile of all values in the sample ordered in ascending order; IQR is the interquartile range.
Remote Sens. 2023, 15, x FOR PEER REVIEW 10 of 26 point cloud after projection ( Figure 7c) and subsequently define the difference between the distance from the outer contour point to the trunk centroid and the trunk radius of that layer as the CR. The trunk detection method is described in Section 2.5.2. TLS takes a bottom-up scanning approach; so, some 'noise' will inevitably appear in the upper part of the tree crown, resulting in an abnormally extracted CR. In this study, we applied the 'boxplot method' [85] to remove outliers from the extracted CR. After removing the outliers from the CR, the largest radius was selected as the maximum CR for that layer. The principle of the 'boxplot method' to detect outliers is to label the points above the upper limit and below the lower limit as outliers (Figure 8). The formulas for calculating the upper and lower limits are as follows: where 3 is the third quartile, which is also known as the "upper quartile" and is equal to the 75th percentile of all values in the sample ordered in ascending order; 1 is the first quartile, which is also known as the "lower quartile" and is equal to the 25th percentile of all values in the sample ordered in ascending order; is the interquartile range. To reduce the influence of the trees' height on the extraction accuracy, the maxim CR extraction accuracy was analyzed in this study in the RDINC range [86]. RDI ranges were grouped at 0.05 intervals, and then the radius values calculated from the la est branches of each whorl in 30 sampled trees served as reference values to compare a

Accuracy Analysis of the Extracted CR
To reduce the influence of the trees' height on the extraction accuracy, the maximum CR extraction accuracy was analyzed in this study in the RDINC range [86]. RDINC ranges were grouped at 0.05 intervals, and then the radius values calculated from the largest branches of each whorl in 30 sampled trees served as reference values to compare and analyze the maximum CR actually measured in each RDINC range and the maximum CR extracted from the point cloud data. If the largest branches of two whorls are too close to each other in the actually measured data of the sampled trees, the maximum CR values calculated from these two branches will fall into the same RDINC range, and only the largest CR of the two will be compared in this case. If the maximum branches of two whorls are too far apart, some RDINC ranges will have no actually measured radius values. Therefore, not all extracted CRs have actually measured values for comparison. In this study, only those parts of the actually measured and extracted values that correspond to each other were compared. The evaluation indexes of the extraction accuracy are as follows: where P is the extraction accuracy, n is the number of samples, X i is the extracted CR, and x i is the field-measured CR.

Development of the Crown Vertical Profile Model
The Weibull function is a popular function used in forest modeling [87,88]. Ferrarese [8] and Gao [89] et al. have modified the Weibull function in order to make it more applicable to the trend of CVP variation in different tree species, respectively. In this study, a modified 3-parameter Weibull function was applied to model the Korean pine CVP. The formula for the 3-parameter Weibull function is as follows: where CR is the crown radius, a, b, and c are the model parameters, and RDI NC is the relative depth into the crown. The shape of the CVP varies for different sizes of trees and has been reparametrized to improve the generalizability of the model [89]. We fitted the parameters of Model (6) with the measured data of the 616 largest branches of 30 sampled trees, and then analyzed the correlation between each parameter and the individual tree factors (DBH, H, CL, HD, and CW). Finally, the factors with a high level of correlation were introduced into Model (6). The formulation of the reparametrized model is as follows: where CR is the crown radius; a 1 , a 2 , b 1 , b 2 , c 1 , c 2 , and c 3 are the model parameters; RDI NC is the relative depth into the crown; the DBH is the diameter at the breast height; the CW is the crown width; the CL is the crown length. The data used for the development of forest growth models often have a hierarchically nested structure [90], such as individual tree data nested within sample plot data. Therefore, the ordinary least squares regression model developed can only reflect the overall average and cannot explain the differences between individuals due to some differences in the study subjects over time or space [91]. In this context, we developed a two-level (plot-level and nested tree level) mixed-effects model [59] to account for the random disturbance to the CVP of individual trees within different plots. In addition, we also developed one-level (plot level and tree level) mixed-effects models to illustrate the differences in the effects of sample plot conditions and tree size on the CVP via comparative analysis, respectively. The mixed-effects model was formulated as follows: where i, j, and k are the plot level, the tree level within the plot level, and the observation of the crown radius, respectively; M is the number of sample plots; M i is the number of sample trees within the ith sample plot; n ij is the number of crown radii on the jth tree within the ith sample plot; Y ijk is the observed value of the kth crown radius on the jth tree within the ith sample plot; A ijk , B i,jk , and C ijk are known as the design matrices; χ is the fixed effect vector; β i is the sample plot level random effect vector; β ij is the sample tree level random effect vector; ε ijk is the error term.

Model Evaluation and Validation
In this study, we performed the random sampling of all data, with 75% of the data being used to fit the model, and the remaining 25% being used to validate the prediction of the model. R 2 a and RMSE were applied to evaluate the fitting effect of the model. The larger R 2 a is and the smaller RMSE is, the better the fitting effect of the model is. MAE and RMAE were applied to validate the prediction effect of the model; the smaller MAE and RMAE are, the better the prediction effect of the model is. The formula for each evaluation index is as follows: where y i is the observed value,ŷ i is the predicted value, n is the number of samples, p is the number of model parameters, and y i is the average of all the observed values.
Regarding the mixed-effects model, the fixed-effects component was validated in the same way as for the OLS model, with an emphasis on the prediction of random parameters. We applied a quadratic sampling method to calculate random parameter values and used the empirical best linear unbiased predictor (EBLUP) [92] to calculate random effects in the model validation according to the following formula.
whereβ k is the random effect parameter;D is the variance-covariance matrix of the random effect parameter;Ẑ k is the design matrix;R k is the within-group variance-covariance matrix; andε k is the difference between the observed and predicted values calculated based on the fixed-effect parameter. Then, the results were evaluated using corresponding statistical indicators.

HCB Extraction Results
We extracted the HCBs from a total of 283 trees from six sample plots, and Figure 9 shows the results of the extracted HCBs compared with the field-measured HCBs. In general, there was excellent correlation between the TLS-based extracted HCBs and the field measurements, with an R 2 of 0.905 and an RMSE of 0.52 m. All the scatter points are mainly concentrated in the lower part of the 1:1 reference line, which proves that the extracted HCBs are generally lower than the HCBs measured in the field are. Table 3 shows the statistical information from the HCB extraction results from different plots, and the results demonstrate that there is minimal difference in the extraction accuracies between different plots, with P being higher than 94.4%, and the biases all being negative. These results prove that the HCBs are underestimated in all plots.

HCB Extraction Results
We extracted the HCBs from a total of 283 trees from six sample plots, and shows the results of the extracted HCBs compared with the field-measured HCB eral, there was excellent correlation between the TLS-based extracted HCBs and measurements, with an R 2 of 0.905 and an RMSE of 0.52 m. All the scatter points a concentrated in the lower part of the 1:1 reference line, which proves that the HCBs are generally lower than the HCBs measured in the field are. Table 3 s statistical information from the HCB extraction results from different plots, and t demonstrate that there is minimal difference in the extraction accuracies between plots, with P being higher than 94.4%, and the biases all being negative. The prove that the HCBs are underestimated in all plots.

CR Extraction Results
In this study, we extracted a total of 392 CRs from the point cloud data of 30 sampled trees and 3757 CRs from the point cloud data of 283 individual trees from six sample plots. The analysis of the extraction accuracy was then carried out using the branch factor data measured in the field from 30 sampled trees as reference values. The results of the extracted CRs based on point cloud data are shown in Table 4, where we applied a random sampling method to group the data into modeling and validation datasets at a 3:1 ratio.

Removal of Outliers of CRs
We evaluated the differences between the directly extracted CRs (Method 1) and the field-measured CRs and the differences between the extracted CRs after removing outliers (Method 2) and the field-measured CRs, respectively. The detailed results are reflected in Figure 10. Overall, the trend of the CRs extracted via the two methods is consistent with the trend of the CVP shape along the RDINC. The differences were mainly in an RDINC range from 0 to 0.3 (Figure 10c). More specifically, in the RDINC range of 0~0.03, the directly extracted CRs are greater than those of the field measurement (Figure 10a), while the extracted CRs after the outliers were removed are already very similar to the field measurements (Figure 10b). In the end, we took the extracted CRs after removing the outliers as the result.
trees and 3757 CRs from the point cloud data of 283 individual trees from six sample plots. The analysis of the extraction accuracy was then carried out using the branch factor data measured in the field from 30 sampled trees as reference values. The results of the extracted CRs based on point cloud data are shown in Table 4, where we applied a random sampling method to group the data into modeling and validation datasets at a 3:1 ratio. We evaluated the differences between the directly extracted CRs (Method 1) and the field-measured CRs and the differences between the extracted CRs after removing outliers (Method 2) and the field-measured CRs, respectively. The detailed results are reflected in Figure 10. Overall, the trend of the CRs extracted via the two methods is consistent with the trend of the CVP shape along the RDINC. The differences were mainly in an RDINC range from 0 to 0.3 (Figure 10c). More specifically, in the RDINC range of 0~0.03, the directly extracted CRs are greater than those of the field measurement (Figure 10a), while the extracted CRs after the outliers were removed are already very similar to the field measurements ( Figure 10b). In the end, we took the extracted CRs after removing the outliers as the result.  accuracy of the CRs is 90.12%. There are variations in the extraction accuracy of CRs in different RDINC ranges. The best RDINC range is 0.15-1, and the extraction accuracies of CRs in this range are all stable at approximately 90%; the worst range is 0-0.15, with the extraction accuracies of CRs in this range all being lower than 80%. The lowest extraction accuracy of CRs is only 73.24% in the range of 0-0.05.  Figure 11 illustrates the relationship between the extracted CRs and the field-measured CRs. As shown in the figure, the scatter points are more concentrated near the 1:1 reference line, which proves that there is a strong correlation between the CRs extracted from the point cloud data and the field-measured CRs. An additional fact is that, overall, all the scatter points are mainly distributed in the upper part of the 1:1 reference line, which proves that the extracted CRs from the point cloud data led to an overestimation of the CRs.

Mixed-Effects Model
Mixed-effects model parameters were solved using the 'nlme' package in R software. We used all combinations of random effect parameters (128 in total), and the model failed to converge when the number of random effect parameters was greater than three. Table   Figure 11. Comparison of TLS-based extracted CR with field-measured CR.

Mixed-Effects Model
Mixed-effects model parameters were solved using the 'nlme' package in R software. We used all combinations of random effect parameters (128 in total), and the model failed to converge when the number of random effect parameters was greater than three. Table 6 reports the evaluation of the optimal model when different combinations of the number of random effect parameters were used. Based on the AIC minimum principle [90], we selected the optimal combination of parameters for each mixed-effects model. At the plot level, the model fit performed best when the random effects were applied to parameters, c1, c2, and c3. At the tree level, the model fit performed best when random effects were applied to parameters, b1 and c2. The model fit performed best when random effects were applied to parameters, a2, b1, and c2 considering both the plot and tree levels. In addition, the LRT test results demonstrated a significant difference between the base model and the mixed-effects model (p < 0.001).  Table 7 shows the parameter estimators and fit statistics of the optimal model. The parameter estimators of all models passed the significance test. It was observed from the results of the fitted statistics that the R a 2 of the model was significantly improved and the RMSE was significantly reduced after the introduction of random effects, which proves that the mixed-effects model significantly improves the model fitting effect. As shown by the comparison results of considering only the plot effect, tree effect, and both the plot and tree effects, the fitting effect of both the tree effect and two-level mixed effect is significantly better than that of the plot effect. In addition, the fitting effect of the two-level mixed effect is also improved compared with that of the tree effect.

Model Validation
In this study, we visualize the residual results of the base model and the mixed-effects model, and the results are shown in Figure 12. It is clear from the results that the residual distribution range of the base model is significantly larger than that of the mixed-effects models, proving that the model's predictive power is significantly improved by introducing random effects. In addition, it is worth noting that the two-level mixed-effects model performs better than the one-level (plot-level or tree-level) mixed-effects model does, with a more concentrated distribution of residuals. Considering plot-level or tree-level random effects alone can improve the predictive power of the model, but only a combination of plot-and tree-level random effects can maximize the performance of the model.

Discussion
In the traditional approach, destructive harvesting is often required to obtain branch factors to simulate trends in the shape of the CVP [9,57,58,61,89]. TLS, as an active remote sensing technique, can truly and effectively restore the 3D structure of individual trees using massive point cloud data [33,48,73,93], which provides a new perspective for describing the trees' CVP. In this study, we propose a method to automatically extract crown radii at different heights based on TLS point cloud data. The validity of the method was then verified using field measurement data from destructive harvesting. Finally, we developed a two-level nonlinear mixed-effects model for simulating the changing trend of planted Korean pine CVPs using the CRs extracted from point cloud data.
The focus of the study in this paper is on the crown portion of the tree; so, it is essential to accurately determine the HCBs. Numerous studies have demonstrated the effectiveness of automatic HCB detection using ALS point cloud data at both stand and individual tree scales [11,71,72]. In contrast, the automatic detection of HCBs using TLS point cloud data has not been reported very often [94]. In this study, HCBs were automatically detected in the plantation based on TLS point cloud data, and the results of the study showed a good correlation between the detected HCBs and the field measurements (Figure 9), which is consistent with the reported results in the literature [11]; however, we produced better R 2 and RMSE values. This is expected because TLS is more suitable for the accurate extraction of HCB information from forest stands because of its ability to pro-  Table 8 shows the results of verifying the predictive power of the base model and the mixed-effects model using independent validation datasets. It can be observed from the table that the two-level mixed-effects model has the best predictive power. The Each model was ranked as follows according to their predictive power: two-level mixed-effects model > tree-level mixed-effects model > plot-level mixed-effects model > base model. The two-level mixed-effects model increased the R a 2 by 3.94%, 10.64%, and 13.8%, reduced the RMSE by 4.61%, 18.56%, and 23.46%, reduced the MAE by 6.48%, 19.52%, and 24.48%, and reduced the RMAE by 9.41%, 17.83%, and 19.37%, respectively, compared to those of the other three models.

Discussion
In the traditional approach, destructive harvesting is often required to obtain branch factors to simulate trends in the shape of the CVP [9,57,58,61,89]. TLS, as an active remote sensing technique, can truly and effectively restore the 3D structure of individual trees using massive point cloud data [33,48,73,93], which provides a new perspective for describing the trees' CVP. In this study, we propose a method to automatically extract crown radii at different heights based on TLS point cloud data. The validity of the method was then verified using field measurement data from destructive harvesting. Finally, we developed a two-level nonlinear mixed-effects model for simulating the changing trend of planted Korean pine CVPs using the CRs extracted from point cloud data.
The focus of the study in this paper is on the crown portion of the tree; so, it is essential to accurately determine the HCBs. Numerous studies have demonstrated the effectiveness of automatic HCB detection using ALS point cloud data at both stand and individual tree scales [11,71,72]. In contrast, the automatic detection of HCBs using TLS point cloud data has not been reported very often [94]. In this study, HCBs were automatically detected in the plantation based on TLS point cloud data, and the results of the study showed a good correlation between the detected HCBs and the field measurements (Figure 9), which is consistent with the reported results in the literature [11]; however, we produced better R 2 and RMSE values. This is expected because TLS is more suitable for the accurate extraction of HCB information from forest stands because of its ability to produce a higher point cloud density in the lower forest canopy compared to that of ALS [95]. The results of HCB detection in sample plots with different stand densities (Table 3) indicated that the extraction method was not affected by stand density and was highly applicable and robust. The detection of HCBs via TLS underestimated the HCBs overall, with a bias of −0.19 m. One reasonable explanation is interference from the presence of dead branches on individual trees [8]. Wang [96] and Béland [97] have conducted research and reported that geometric or radiometric intensity features can be used to distinguish between leaves and branches on coniferous and broadleaf trees. Therefore, in future studies, considering the automatic detection of HCBs using geometric [96,98] and intensity [97,99] information will greatly improve the accuracy of HCB extraction.
In previous literature that used point cloud data to study CVP models, CRs were calculated as the horizontal distance from the crown's outermost point cloud to the vertical line where the apex of the tree is located [1]. This ignores the effect of the degree of tree tilt and results in a large error between the calculated CR values and the actual CR values. In this study, we simulated the measurement of destructive harvesting in the field [9,58,61], in which the CRs were calculated as the distance from the crown's outermost point cloud to the trunk surface. We detected circles from the crown point cloud projection using the RANSAC algorithm and treated them as tree trunks; the RANSAC algorithm has been demonstrated to be very suitable for extracting trunk information [74,75,100]. This reduced the error caused by the large degree of tree tilt to a certain extent and, at the same time, it can be better compared with the field measurement data to analyze the accuracy of CR extraction.
The results of this paper indicate a strong correlation between the CRs extracted using TLS data and actual field measurement data ( Figure 11). Although TLS overestimates CRs in general, such a bias is within acceptable limits considering the time and effort required to obtain CRs in the field. The extraction results of the maximum CR in different RDINC ranges were generally better, but there were still large differences in the extraction accuracies between the crown top and the crown bottom, where the extraction accuracy of the crown bottom was better than that of the crown top ( Figure 10). Quan [101] and Xu [102] et al. used UAV-LS point cloud data to extract CRs at different heights of Larix olgensis and Cunninghamia lanceolata (Lamb.) Hook., respectively, but they did not provide specific extraction accuracies because they did not have actual field measurement data for reference. Li [86] et al. extracted the trunk diameters of Larix olgensis at different heights in the vertical direction based on TLS data, and they found that the accuracy of the trunk diameter extraction decreased sharply when the relative height exceeded 0.9. In this study, the extraction accuracy of CRs at different heights was analyzed explicitly with reference to the branch factor data measured in a field containing 30 planted Korean pine trees harvested via destructive harvesting ( Table 2). The extraction accuracy was low in the RDINC range of 0-0.15, stabilized from an RDINC of 0.15, and it remained at approximately 91%, which is consistent with the findings of Li [86]. Although trunk diameter extraction was not performed in this study, it was studied in the vertical direction using point cloud data; so, it can be compared with the results of Li's study for mutual verification. TLS applies a bottom-up scanning method [95,103], which gradually becomes less effective as the trees become taller. Due to the mutual occlusion of tree crowns, it further increases the noise of the top crown point cloud, which eventually leads to a decrease in the accuracy of the CRs extracted at the top of the tree crown. It has been demonstrated that the fusion of UAV-LS and TLS data can restore the crown structure more clearly and completely [104]. Therefore, it will be our future goal to improve the extraction accuracy of tree crown structure information by combining LiDAR data from multiple platforms in future research.
Most of the data for forest stand growth modeling have a hierarchical structure, and nonlinear mixed-effects models are an effective approach to address this problem [90,91] and have been widely used in the development of tree CVP models [61]. The subject of this paper was 283 individual trees from six sample plots with different attributes, which makes it very suitable for developing a mixed-effects model. With the aim of adequately explaining the random disturbances caused by sample plots and sample trees on the trees' CVP, we developed tree-level, plot-level, and two-level mixed-effects models. In the process of fitting a mixed-effects model, the model can have difficulty converging as the random effects increase [105]. We have adequately considered all combinations of random effect parameters, and the results show that the model does not converge when there are more than three random effect parameters. The mixed-effects model has a significantly better fit than the base model does (Table 8), demonstrating that the introduction of random effects can greatly improve the prediction accuracy of the model, which is consistent with the research results of Gao [61] and Liu [105] et al. Ferrarese [8] et al. reported that the CVP was not strongly influenced by the tree size and sample site conditions. This is because their study samples were biased toward more open stand conditions and bare canopies. The plot-level mixed-effects model we developed reflects differences between the sample plots, such as the stand density, indicating that all sample plot factors have some effect on the CVP. The tree-level mixed-effects model outperformed the sample plot-level mixed-effects model, which indicates that the differences in the CVP models were mainly due to random disturbances caused by different size sample trees. The two-level mixed-effects model had the best performance compared to those of the plot-level mixed-effects model and the tree-level mixed-effects model, which demonstrates that considering both the plot-level and tree-level effects was more effective at adequately explaining the random influence on the tree crown outer profile model.
Roeh [106] and Cluzeau [58] chose to use an indirect method to describe the CVP, as they considered that there were many difficulties in directly measuring the crown radius at different heights. TLS can provide detailed information regarding the crown structure, which is very suitable for directly describing the CVP [8], and thus, successfully overcomes the abovementioned difficulties. The previously developed CVP models based on field measurements usually require the destructive harvesting of trees, and as a result, the study sample is very limited, e.g., (Marshall [9], 36 western hemlock; Doruska [107], 34 loblolly pine; Linnell Nemec [49], 45 amabilis fir, 60 lodgepole pine, and 60 white spruce). In this study, we automatically extracted the crown radii of different heights from 283 individual trees in six sample plots using TLS without destructive sampling. It produced a large sample size and a representative dataset. A two-level mixed-effects model for the CVP was then developed based on this dataset, which greatly improved the predictive capability of the model. Overall, the results of this paper demonstrate that TLS can replace the traditional method of measuring crown structures and improve the ability to describe the CVP.
It is worth emphasizing that our proposed method is more applicable to complete crown point clouds. The subject of this study is a planted Korean pine forest, where there is less crossover between crowns; so, the crown point cloud data we obtained are complete. Using the proposed method, the crown radius of each individual tree can be extracted more successfully, which is conducive to the development of high-precision crown vertical profile models. In more complex mixed forests, crowns are difficult to be completely segmented due to the serious crossover phenomenon between the crowns. Since the crown point cloud itself is incomplete, the extracted crown radius also does not match the reality, which will seriously affect the prediction accuracy of the developed model. Therefore, a more robust crown segmentation algorithm suitable for high-density mixed forests would improve the applicability of the proposed approach. It is also an interesting topic for our future research.

Conclusions
TLS has great potential for accurately measuring tree crown structures and can also be used as a measurement tool in the field. In this study, we demonstrate the feasibility of the application of TLS for the development of plot-level CVP models. We propose an approach that enables the nondestructive measurement of CRs at different heights. It enables us to acquire larger amounts of data in less time and, in addition, to repeat the acquisition of data, which introduces a way to monitor dynamic changes in forest crown structures. The experimental results regarding 30 felled trees showed that the average extraction accuracy of the method was 90.12%. The extraction accuracy of the top part of the tree is poor due to the limitation of the device itself, and future research will focus on a combined UAV and TLS to improve the extraction accuracy of the CRs of the top part of the tree.
Finally, we developed mixed-effects models based on automatically extracted crown radii from 283 individual trees in six sample plots, and the results showed that the introduction of random effects significantly improved the prediction accuracy of the models. Each model was ranked as follows according to their predictive capability: two-level mixedeffects model > tree-level mixed-effects model > plot-level mixed-effects model > base model. It was demonstrated that the tree size was the main factor influencing the shape of the CVP. However, considering both the sample plot conditions and tree size could more adequately explain the differences in the CVP. Therefore, it is essential to develop CVP models on a regional scale.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.