Trunk Model Establishment and Parameter Estimation for a Single Tree Using Multistation Terrestrial Laser Scanning

Parameter estimation and effective tree trunk model establishment are significant to tree biomass measurements and three-dimensional (3D) reconstruction. In this study, we proposed a trunk model establishment and parameter estimation method for a single tree through multistation terrestrial laser scanning (TLS). This approach obtains a more accurate and dense point cloud and yields several improvements. First, to correct the tree position, the angle between the trunk centreline and the vertical line was extracted to normalize the tree point clouds. Subsequently, some tree model parameters were accurately calculated. The tree height was extracted based on individual differences rather than directly adopting the Digital Elevation Model (DEM) model; The B-Spline curve fitting was added to the diameters at breast height (DBH) measure, which better determined the shape of the DBH compared with circle fitting; The CBH was measured by multi-view voxel projection instead of Single-view fixed threshold segmentation, which confirmed the position of the CBH more accurately. To highlight the effectiveness and pertinence of the proposed method, we analysed three plots that the proposed method can achieve more accurate results, where the mean relative error (MRE) of the tree height, DBH and CBH were 1.89%, 1.74% and 2.91% satisfying the accuracy standard specified by the technical indicators (±5%). The experimental results validate the trunk model accuracy and the effectiveness of the parameter measurements (close to the field measurements). The proposed method is significantly different from the accuracy of the previous method. The proposed method has the best accuracy in Multistation Terrestrial Laser Scanning. Hence the proposed method has outstanding performance in the high-precision extraction of trunk parameters and the construction of trunk models which will perform forest inventory and construction of virtual forest more accurately.


I. INTRODUCTION
A. BACKGROUND Estimation of trunk parameters such as tree diameter or height is of great significance to tree growth state and biomass estimation analyses. The tree volume and biomass allometric equations have been established with extensive destructive sampling of trees [1]. Tree models require many dependable samples [2], but the number and spatial extent of samples collected in the field are always limited. These limitations The associate editor coordinating the review of this manuscript and approving it for publication was Lefei Zhang . primarily stem from the significant cost and labour expense involved in destructive sampling. With the rapid development of measurement equipment, the non-destructive measurement method has been widely used in tree surveying, which not only reduces the labour cost but also avoids considerable damage.
In recent years, many kinds of non-destructive tree measurement methods including hyperspectral [3], [4], infrared [5], two-dimensional image [6], ultrasonic [7], laser point cloud [8] and other methods have been widely used. UAV-sensed high-resolution imagery [9], 2D laser point clouds [10], Feller-buncher cutting sounds [11], VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ ARBOTOM R equipment (Measure root sound waves) [12] and Digital images [13], [14] have been adopted to measure the diameters at breast height (DBH). Jarocińska et al. [15] used AISA hyperspectral images to analyze the locations, species and levels of discolouration and defoliation of trees. Hu et al. [16] proposed a method for measuring the leaf geometric characteristics of poplar seedlings based on 3D visualization via the use of time-of-flight (ToF) and digital cameras. A new method was proposed for calculating the tree volume model with high precision by using the electronic theodolite [17]. Mugasha et al. [18] Constructed a variety of D-H relationship equations. Extents in which climate and forest stand variables explain the variation in H-D allometry were also assessed. However, hyperspectral and infrared are more suitable for stand biomass, chlorophyll index [19], [20], pest and disease monitoring [21], etc. UAV images or satellite high-precision images are suitable for biomass of forest stands and tend to check large areas of resources [22]. And, their measurement accuracy and data stability are difficult to guarantee compared with point cloud data. Measuring tree parameters through lidar data can realize forest monitoring (including tree species classification [23], forest structure [24], canopy height [25], parameter extraction, etc.) and mapping [26]. The point cloud data is particularly prominent in highprecision data extraction and morphological model construction. However, point cloud data is particularly prominent in high-precision data extraction and morphological model construction which can even clearly show the shape of a leaf [27]. Taking a single tree as an example, point cloud data can reflect the individual differences of each tree [28] and completely show the shape of a tree from the bottom to the top [29] which is an advantage over hyperspectral and high-score images. Through laser point cloud data, the tree parameters (accuracy in millimeters) can be collected [30], and the tree model can be reconstructed to build a virtual forest [31] for scene restoration and evolution.
The current carrying means of the laser system include airborne laser scanning (ALS), terrestrial laser scanning (TLS) and mobile laser scanning (MLS). ALS is always used to detect individual trees [39], [40] measure tree stem radii [41] and classify tree species [42]. Ramiya et al. [39] proposed an object-based labelling framework for delineating individual tree canopy clusters in urban green spaces from an airborne LiDAR point cloud. Fritz et al. [41] measured the tree stem radius using aerial SfM-MVS above an oak-dominated forest. The process assigns the points in the 3D point cloud image to each tree in the upward direction and separates not only tree trunks but also tree canopies. It is desirable to automatically detect each tree and analyse its structural parameters, which cannot provide accurate information about beam intensity. Torabzadeh et al. [42] fused imaging spectroscopy (IS) data with leaf-on and leaf-off small-footprint airborne laser scanning (ALS) data for tree species identification in a dense temperate forest in Switzerland. Although ALS has the advantage of a wide scanning range and fast scanning speed, some trunk information is easily lost in layered forests.
MLS involves fixing a laser sensor to a car or backpack. Xu et al. [43] proposed a method for individual stem detection in residential environments with MLS data. Li et al. [44] presented a dual growing method for the automatic extraction of individual trees from MLS data according to the general constitution of trees. A mobile laser scanner is flexible and convenient, but the point cloud is difficult to register.
TLS has also been widely used in forestry to automate tree diameter and height measurements [45], [46], biomass estimation [47], tree species [48] and stem detection [49]. Korň et al. [45] proposed a method of tree diameter estimation from TLS by fitting a circle. Tree trunks can be detected using the Hough transform or by conducting circle (cylinder) fitting to a cross-section at a certain height of 3D point cloud tree images. Cabo et al. [46] developed an automatic method based on the isolation and vertical continuity of stems to identify tree stems and estimate tree heights and diameters from TLS data. First, a height-normalized version of the point cloud was created, from which the stems were individualized. Subsequently, an iterative process was applied to the points at breast height for estimating diameters, and tree heights were calculated after denoising and clustering the points of each tree. Conto et al. [49] proposed stem denoising and modelling algorithms on single tree point clouds from terrestrial laser scanning, which were adapted from common TLS stem detection techniques. The researchers analysed the performances of three different stem denoising methods and three different stem modelling methods on TLS point clouds containing single trees and carried out detailed comparison experiments. The results indicated that the combination circle model of the Hough transformation stem denoising method and the iteratively reweighted total least squares modelling method had the best performance for stem modelling. However, TLS faces some challenges in measuring individual tree parameters. TLS is a discrete sampling technology, and its point density is range-dependent. These features mean that stems far from the scanner will be mapped with less details than those close to the scanner. A single laser scan can only observe one part of a stem, and the geometry is unfavourable from a mathematical perspective, which causes errors in the estimation of stem radius and location. Moreover, a single laser scan cannot solve the mutual occlusion. As a result, occluded stem sections or even full stems will fail to be mapped from a single scan. Now, new TLS devices are small, lightweight and fast scanning, so multiple scans from different standpoints can be easily realized in a forest environment. The multistation laser scanning can solve the problem of occlusion since this scanner can scan all objects.
In this study, the parameter measurements were presented based on multistation laser scanning. Trunk parameters, such as DBH and height, are the important parameters required for forest inventory, especially in the 3D reconstruction and biomass estimation of trees. The conventional instruments for the measurement of forest plots are callipers or girth tapes for trunk diameters and hypsometers for total tree heights. In recent years, laser data have offered new opportunities to obtain tree diameter and height [32], [33], [45], [50]; There are many studies on tree diameter estimation, which include circle fitting [45], [50], circle-ellipse fitting [51] and Hough transformation [52].
However, several problems remain in the aforementioned studies. (1) Some methods based on ALS may lose points on part of the trunk. (2) Some methods based on a single scan TLS may have unfixed point cloud densities according to the distance between the tree and laser. (3) Some methods are time consuming and destructive. (4) Some methods are inaccurate when fitting the DBH. (5) Previous TLS equipment is bulky and difficult to carry. (6) Some methods are influenced by the number of scans of the tree, the scanning angle, shadowing by trees, low vegetation and terrain.

B. OBJECTIVES
In this study, we proposed a non-destructive measurement method to establish trunk model and estimate parameter for a single tree through multistation terrestrial laser scanning (TLS). Three main tree taxa in Nanjing Forestry University were selected as measurement objects, and Research Institute of Forest Resources Information Techniques, Chinese Academy of Forestry has assisted the integration and splicing of point clouds. The main objective is to develop an algorithm for quickly segmenting stems and estimating the tree height, DBH and CBH from a single tree point cloud. Based on the research [45], [53], a single tree point cloud was extracted by combining the Faster R-CNN with the nearest neighbour clustering algorithm [53]. A modified circle fitting method was introduced to estimate DBH, CBH and trunk height based on the point cloud obtained by multistation laser scanning. The objectives mainly included (1) calculating tree height from the difference between the highest point and lowest point; (2) constructing the DBH model by combining circle fitting based on adaptive differential evolution algorithm (ADE) and cubic B-spline fitting; (3) segmenting the stem from a single tree point cloud using the voxel projection method and measuring the CBH; and (4) establishing the trunk tape model using the polynomial-fitting method.
Through the above goals, the main contributions of this article include: (1) According to Multistation Terrestrial Laser Scanning, we get fast, accurate and complete point cloud data. This method can help us obtain more accurate forestry resource data.
(2) The individual differences, B-spline fitting, and multiview voxel projection have been added to trunk parameter extraction, this allows us to obtain more accurate tree morphology parameters. DEM does not fully consider the relationship between trees and terrain. Therefore, we propose a centerline intersection method to extract the lowest point in the tree height calculation; The B-spline curve combined with ADE replaces the traditional circular fitting, which further improves the accuracy of DBH parameters and better builds the trunk shape;By increasing the judgment angle and adaptive threshold, the influence of small branches on CBH parameter extraction is reduced, and the accuracy of CBH judgment is improved. With the extraction of the crown parameters to be completed, the virtual forest will be better constructed and the trees will be more realistic in the future.
(3) In addition, this research can also be used for smallscale forest inventory and reconstruction and inversion of tree morphological parameters.
The rest of this article is organized as follows. Section II explains the proposed parameter extraction method. Experimental results are displayed in Section III to demonstrate the performance of the proposed method. Finally, Section IV gives someconclusions and suggestions.

A. MATERIALS
The TLS can obtain massive amounts of point clouds efficiently and quickly without any intervention. A TLS (Trimble TX8) was used to collect plant point cloud data. The weight, measurement accuracy, and distance range are 11.0 kg, 1 mm, and 120 m, respectively, and the point scanning speed is 1000000 points per second. The main scanning principle of TX8 is to scan the tape with a rotating prism. The TX8 has a high scanning speed and only takes three minutes to complete a scanning task. The selected experimental trees were located on the campus of Nanjing Forestry University (32 • 08'N, 118 • 81'E). Three plots were selected as study sites. Rows of poplar trees, liriodendron Chinensis and ginkgo were all planted at intervals of 5 m-10 m, and their diameters were approximately 25 to 65 cm, 10 to 40 cm and 10 to 35 cm, respectively. The scene is shown in Fig. 1. To obtain accurate information on each target tree, a multistation scanning method was adopted to scan each tree face, where we chose different station layouts according to the size of the trees and the degree of occlusion. Taking the poplar plot as an example, we adopted 42 stations to scan the poplar forest, and the station arrangement is shown in Fig. 2, where ' ' represents the station and ' ' represents the poplar tree. After scanning, the data scanned from different stations were integrated into a single coordinate system by registering the point clouds with Trimble Realworks. This program registered the point cloud of each station with other stations and determined whether the point cloud registration was reliable based on the surrounding positioning points, errors, credibility and coincidence. Taking station 1 in a poplar plantation as the example, we matched station 1 with the other 41 stations. Then, we fused 4 matched stations whose errors were less than 10 mm and credibility was 100% with station 1 (other stations fused in a similar way). Finally, the coincident point clouds were synthesized as full coverage scanning data from the object trees, which included three plots. The study was  executed on a computer equipped with a 12 core CPU running at 3.70 Hz with 32 GB of RAM, a Nvidia GTX1080Ti GPU processor and a Windows 10 operating system. The programming environments were Visual Studio 2019, PCL-1.10.1 and Unity3d.

B. METHOD 1) FRAMEWORK OF THE METHOD
Trunk parameters such as height, DBH, and CBH are important for calculating the biomass and realizing the 3D reconstruction of a tree. Wang et al. [53] extracted the individual tree segmentation by combining the Faster R-CNN of the deep learning method, which recognized the locations of the tree trunks, and the nearest neighbour clustering algorithm, which extracted the individual tree crown segmentation. Koreň et al. [45] proposed a method for tree diameter estimation from TLS by the circle fitting method. Tree trunks can be detected using the Hough transform or by conducting a circle (cylinder) fitting to a cross-section at a certain height of the 3D point cloud tree images. Bu and Wang [51] proposed an adaptive circle-ellipse fitting method based on the point cloud transect. Because of the complexity of the forest environment and the trunk shape, it is difficult to describe the shape of a trunk with a circle or ellipse. In this study, we proposed a novel method to extract the parameters of the trunk. The segmentation method in [53] was used to extract a single tree, and the statistical filtering method in [16] was employed to filter noises and outliers. Because most of the trees are not completely vertical to the ground, the centreline of the preliminary trunk part is fitted and tree point clouds are corrected by the angle between the line projecting the centreline to the XY plane and the X axis and another angle between the line projecting the centreline to the XZ plane and the Z axis. Then, the tree height is the Z coordinate difference between the optimal highest and lowest points of the tree, where the highest point of the tree is found by the local extreme method, and the lowest point of the tree is the trunk centreline and the fitted ground plane intersection. Next, a cross-section, which is approximately 1.3 m above the terrain, is extracted from the trunk point clouds. Subsequently, these point clouds are projected onto the XY plane, and the circle fitting based on ADE and cubic B-spline fitting methods are used to obtain the DBH. After extracting the DBH, the multiangle voxel projection method with adaptive threshold is introduced to segment the trunk and measure the crown bottom height (CBH). To assist in constructing the trunk model, the taper equation is fitted by extracting the taper parameter of the stem. The method framework is shown in Fig. 3.

2) HEIGHT-NORMALIZED POINT CLOUD
In the natural environment, woodlands are often rugged, and tree growth is heliotropic and random. Hence, it is difficult for trees to grow vertically, even in planted forests. To reduce the errors caused by the irregular growth of trees and rugged woodlands, a height-normalized point cloud is required. Cabo et al. [46] built a direct digital elevation model and generated flat ground by setting a threshold for the expected altitude. The triangulated irregular network (TIN) method was used to establish a digital terrain model (DTM), where the z-value of each point was normalized to normalize the point cloud [54]. In these studies, the error was reduced by flat ground, while the method lacked focus for a single tree. This study introduced a centreline correction method to normalize the tree point clouds by using trees and ground simultaneously. First, the centreline of the trunk needed to be obtained. The cylinder fitting method [45] was employed to extract partial stem points with z coordinates that ranged from 0.395 m to 0.405 m and 1.295 m to 1.305 m. Hence, two centres (x 1c , y 1c , z 1c ) and (x 2c , y 2c , z 2c ) at different heights and the centreline between these two centres were obtained. Second, the centreline is projected to the XY plane. All the points are rotated by Equation 1 according to the angle difference θ 1 between the projecting centreline and the X axis in the plane, as shown in Fig. 4a. In succession, the centreline is projected to the XZ plane, and the new points are rotated θ 2 in Equation 2, as shown in Fig. 4b. Finally, the point cloud after correction is obtained.
The processing of the tree height normalization is as follows: (1)Setting all cloud points as P, (x i , y i , z i ) ∈ P, i = 1, 2, . . . , N , N is the number of tree points. Then, the two centres (x 1c , y 1c , z 1c ) whose z coordinate range from 0.395 m to 0.405 m and (x 2c , y 2c , z 2c ) and whose z coordinate range from 1.295 m to 1.305 m are extracted by cylinder fitting.
(2)The centreline connecting (x 1c ,y 1c ,z 1c ) and (x 2c ,y 2c , z 2c ) is projected to the XY and XZ planes, respectively. The angle θ 1 is the tilt angle between the centreline projection and the X axis, and the angle θ 2 is the tilt angle between the centreline projection and the Z axis, as shown in Fig. 4a and b. (3)Calculating the new point clouds as P , where x i , y i , z i ∈ P by rotating the θ 1 angle around the Z axis by using the following matrix: (4)The normalized point clouds are obtained as P , where x i , y i , z i ∈ P by rotating the θ 2 angle around the Y axis by using the following matrix: Fig. 5a illustrates that there is an inclination angle between the trunk and the horizontal plane. After correction using the mentioned processes, the tree trunk is almost vertical to the horizontal plane, as shown in Figure 5b.

3) TREE HEIGHT MEASUREMENT
Tree height is an important parameter for tree visualization and biomass estimation. Previous studies [55], [56] usually used simple and straightforward algorithms to derive the tree height by calculating the distance between the highest laser scanning point and the lowest point. However, the growth environment of trees is complex, and trees are easily covered by bushes or uneven ground surfaces, so it is difficult for the laser altimeter to find the highest and lowest points directly.
To solve these problems, this paper determined the optimized lowest point H min to be the intersection of the centreline of the trunk and the fitted ground plane. Extracting the optimized highest point as H max by the local extreme method [57], the distance between H max and H min is the tree height H as shown in Fig. 6.
The tree height calculation process is as follows: (1) The ground plane is fitted by Equation 3: where x and y are the X-coordinates and Y-coordinates of ground points, respectively, and ω = (a, b) represents the fitting coefficients.
(2) New centres x 1c , y 1c , z 1c and x 2c , y 2c , z 2c , which range from 0.395 m to 0.405 m and 1.295 m to 1.305 m, respectively, are extracted to search the trunk centreline.
(3) The vertical planes ω 1 (perpendicular to the XY plane) and ω 2 (perpendicular to ω 1 ) are calculated through these two centre points. The intersection line of the two vertical planes, which is presented in Equation 4, is the trunk centreline.
where (x , y ) are the X-coordinates and Y-coordinates of the trunk points and ω 1 = (a 1 , b 1 , c 1 ) and ω 2 = (a 2 , b 2 , c 2 ) are the parameters to be fitted. VOLUME 8, 2020 (4) Calculate the intersection (x L , y L , z L ) of the fit ground plane using steps (1), (2) and (3). The centreline of the trunk was denoted as the lowest point H min , and the highest point H max was obtained by the local extreme method. Finally, the tree height H was calculated using the following equation: To obtain the optimal mathematical model above, many optimal algorithms were used for parameter fitting, such as the least square, particle swarm optimization (PSO), gradient descent and evolutionary difference methods. Many experimental results show that the evolutionary difference method is fast. In this study, all parameters ω, ω 1 and ω 2 were fitted using the ADE [58].

4) DBH MEASUREMENT BASED ON THE CIRCLE FITTING BASED ON ADE AND CUBIC B-SPLINE
DBH is an important parameter for building a trunk model. Many studies [30], [54] employed a circle model to establish the DBH model. The processed point clouds in these studies did not usually form a regular circular shape, which led to large errors in fitting the breast diameter. The point clouds collected through the multistation scanning have enough points, where the circle fitting based on ADE method combining the B-spline method was proposed in this paper.
Because there are some noises in the process of point cloud collection, it will affect the accuracy of DBH. The noise filtering method based on the distance density function [16] was first employed to filter the noise points. Fig. 7a illustrates the partial trunk cloud points ranged from 1m to 1.5m. As shown in Fig. 7b, there are some noise points projected from Fig. 7a to XY plane. The noise points are obviously filtered in Fig.7d projected from Fig. 7c filtered using the method [16] to XY plane. Subsequently, the points P DBH P DBH ⊂ P were extracted for DBH fitting by projecting the filtered points  to the XY plane whose Z coordinates ranged from 1.295 m to 1.305 m as shown in Fig. 8a. Then, the circle model is developed as follows: where ν = (x c , y c , R, A), x c , y c R and A are the centre X-coordinate, Y-coordinate, radius and elliptic coefficient, respectively. Then, the ADE [58] is adopted to fit the parameters of ν. The fitting circular curve is illustrated in Fig. 8a. Fig. 8b shows that the ellipse coefficient is added to adaptively fit the circular curve at the same time, which can be beneficial to deal with the case where the trunk DBH is non-circular. However, there are still more points outside the curve. To improve the accuracy of the diameter measurement, cubic B-splines are used to further fit these points. The threshold is set as T D = 0.02R, and the distance d c from each point to the centre of the circle is calculated by the Euclidean distance. If d c > T D , this point is considered to not be on the circle. As shown in Fig. 9a, the ' ' points are not on the circle. If there are more than five ' ' points continuously, the B-spline is further used to fit the extension points on both sides of the ' ' and ' ' points, which is shown in Fig. 9a. Fig. 9b and c show that outliner points generate corresponding control points for curve fitting. Regardless of the adjustment of the circle fitting parameters, some outlier points cannot fit the fitting curve. Using the B-spline fitting method and curve fitting to redraw the curve passing through these points makes the DBH parameter extraction more accurate. Fig. 9d shows that two new curves, instead of the original fitting curve, are obtained. The total fit curve length L and the final DBH as D are calculated by using the following equations: where L 1 , L 2 , · · · , L n are the lengths of the new curves, n is the number of the new curves and L c is the remaining length of the circle fitting.

5) CBH MEASUREMENT BASED ON VOXEL PROJECTION
In recent years, studies [59]- [62] have been devoted to the development of methods and techniques for extracting CBH by using 3D point clouds. Despite the aforementioned progress in extracting CBH, few studies have resolved the interference caused by side branches, and previous studies used a fixed threshold as the segmentation boundary. This paper designed a multiview voxel method, as shown in Fig. 10, and employed an adaptive algorithm [63] to determine the CBH threshold, which could solve the interference caused by side branches, as shown in Fig. 10d. Then Fig. 11  shows the comparison of the effects of fixed threshold voxel segmentation and multi-view adaptive threshold voxel segmentation. Before voxel projection, the voxel grid system is constructed as follows. Every voxel is a building block in the shape of a cuboid, and its geometry is defined by the length (l), width (w), and height (h). The location of a voxel in the 3D voxel grid system is indexed by the column (i), row (j), and layer (k). The valid voxel is defined as the case in which the number of points contained in the cell block of the mesh is greater than or equal to 3. In the multiview voxel method, the voxel projection of each layer (Z k , k = 1, 2, . . . , M , where M is the number of layers) is obtained by constructing 3 cm * 3 cm * 5 cm (l * w * h) valid voxels as K i,j . For each layer, the number of voxel projections along the height direction was calculated, as shown in Fig. 10.
The first step of the multiview voxel method is voxelization based on point clouds divided into 4 quadrants that have been centred on the trunk centre point of each layer. As shown in Fig. 10, statistically, the valid voxels in each layer and each quadrant are projected onto the Z axis and displayed as a voxel projection line chart. Then, the adaptive threshold method is used to obtain the threshold of the voxel projection histogram in each quadrant, and the first point whose projection coordinate exceeds the threshold is H c (c = 1, 2, 3, 4). The smallest one in the four quadrants is the final CBH H c .
where D is the DBH, H is the tree height and [] represents the integer function.
(2) Detect the relationships among all points in each voxel in each layer. The j-th valid voxel of the i-th floor K i,j is VOLUME 8, 2020 determined as follows: where n i=1 P i,j is the number of points in the j-th voxel in the i-th layer.
(3) Sum the valid voxels of each layer as Z i using the following equation: where n is the number of i-th layer voxels, and then, project them to the XZ plane in quadrants, as shown in Fig. 10.
(4) Calculate the threshold (T c ) using an adaptive algorithm and find H 1 , H 2 , H 3 , H 4 , which are the first points of each quadrant larger than the threshold (T c ), and calculate the CBH with Equation 11:

6) TRUNK TAPER MODEL
The crown of a broad-leaved tree usually has dense branches and leaves, which make it difficult for laser scanners to directly obtain information about the crown. Hence, the taper of the stem is an important auxiliary parameter for fitting the trunk model. The covered trunk of the crown can be fit by the taper equation. Dozens of fractal equation models can be divided into the simple taper model, the segmented taper model, and the variable parameter taper model. Type d 2 = f (D, H , h) and type d = f (D, H , h) are the two main types. Considering that the stems of the tree species studied in this paper were straightforward and had relatively simple growth characteristics, a multivariate nonlinear regression analysis was performed to obtain the optimal taper equation of the trees with the ADE algorithm.
Taper point cloud segmentation was performed before fitting the taper curve. The points P t x i , y i , z i ∈ P between H min + 0.4m and H c were extracted. Then, a curved rectangular parallelepiped from the DBH centre point to 0.1 * D on both sides of the X axis and from the DBH centre point to 1.5 * D on both sides of the Y axis was obtained. After projecting these points onto the YZ plane, the pitch equation was employed to obtain the taper fitting curve.
The process of the taper fitting is as follows: Extract points P t x i , y i , z i ∈ P as shown in Fig. 12, which meet the following conditions: Fit the eligible point clouds through the following taper equation and obtain the fitting curve, as shown in Fig. 11.
where h is the relative height from the ground, d is the diameter at the height h of the trunk, D is the diameter of the breast, H is the height of the tree, and a 0 is the parameter to be determined.

III. RESULTS
In this study, numerous experimental analyses were performed to evaluate the effectiveness and practicability of the proposed method. We utilized the Trimble X8 to collect point clouds of the poplar plantation, liriodendron Chinensis plantation and ginkgo plantation at Nanjing Forestry University by multistation scanning samples with different heights and DBH, as shown in Table 1 and Table 2   points. The artificial Liriodendron chinense forest collected 40 stations, and the total point cloud was 129467345 points (The point cloud synthetic credibility is 100% and the initial point cloud error is less than 2 mm). The data collection and fusion used Trimble X8 and Trimble RealWorks respectively. Finally, 93 of 107 trees on the research plot were used in the evaluation of the proposed method. In Multistation Terrestrial Laser Scanning, the average integrity of the point cloud of each tree was 97%. Simultaneously, we used callipers, altimeter (CQG-1 Direct reading altimeter) and a tower ruler to measure the actual DBH, tree height and CBH of the sample trees.
In section III-A, the comparison experiments of optimization algorithms for trunk parameter model fitting were discussed. In section III-B, the performance of the proposed method was analysed by comparing the trunk's measured parameters with the actual parameters.

A. COMPARISON ANALYSIS OF OPTIMIZATION ALGORITHMS
In this study, the tree tilt correction model, the tree height measurement model, the DBH model, the CBH model, and the taper model are established. However, the parameters of every model need the optimization algorithm to approach the optimal value. At present, there are least square (LS) algorithms, gradient descent (GD) algorithms, genetic algorithms (GA), PSO algorithms and ADE algorithms. The algorithms in this article are all done in C++ and visual studio 2019. We used these algorithms to compare and analyse the parameters of the above models. A large number of experimental results show that the ADE algorithm is not only fast but also closer to the actual value. Taking the DBH fitting as an example, we use the comparative analysis method to analyse the DBH fitting optimization algorithm and the effectiveness of the algorithm.
We use the method in section II-B-4 to extract the DBH point cloud data (1.295 m-1.305 m) from 93 sample point clouds and project these points to the XY plane. The number of point clouds with different sizes of DBH is different, but for the DBH of the same trunk, the ADE fitting time is the fastest under the condition of the same fitting accuracy. To compare the optimization algorithm, we utilize the circle fitting model for the DBH of these trunks to compare the effectiveness of LS, GD, GA, PSO and ADE. Table 3 indicates the fitting results of different optimization algorithms for the same tree with the actual DBH of 16.426 cm. The last four algorithms obtain the same DBH value (15.760 cm), although ADE has the shortest running time (0.797 s). To improve the accuracy of the DBH measurement, we combine the circle based on ADE model and cubic B-spline model proposed in section II-B-4 to fit the DBH of the tree trunk. When using a regular circle fit, the fitted curve cannot pass through all points. In this case, the point cloud of knots or non-circular part of the tree will usually deviate from the fitted curve. The cubic B-spline fits a new curve to cope with this situation, and this curve will pass through all the knots or non-circular point clouds, thus ensuring that the DBH parameter is closer to the real situation. In Table 3, the fitted DBH value increased from 15.693 cm to 15.929 cm, which is closer to the true value of 16.426 cm. Simultaneously, the method is also efficient, with a fitting time of 0.876 s. Hence, compared with other algorithms combined with the cubic B-spline model, the fitting time is still the shortest, and the fitting DBH is closer to the actual value. Similarly, ADE is also applicable to the optimization of tree height and CBH parameters to approximate the best value.

B. PERFORMANCE OF THE TRUNK PARAMETER MEASUREMENTS
To better verify the effectiveness of the model, we carried out a large number of comparative experimental analyses on the performance of the measurement parameters. This section analyses the performance of the proposed method by comparing the measured parameters of the trunks from the proposed method with the actual parameters. The main parameters include DBH, CBH and tree height, which were verified using root mean square error (RMSE), R 2 , VOLUME 8, 2020 Chi-Square, mean relative error (MRE%), F-test etc by SPSS (Statistical Product and Service Solutions).

1) DBH MEASUREMENT PERFORMANCE
In the process of DBH estimation, the comparison experimental analysis was implemented between the proposed method of DBH measurement and the circle fitting method [45]. The trunk is not always round, especially in the measurement of DBH, and knots, local depressions, etc. occasionally occur. Fig. 13 shows the fitting effect of three different shape point clouds of DBH by using circle fitting and the proposed algorithm. Fig. 13a and b show the point cloud fitting results of local deformations at the DBH of the poplar tree trunk of approximately 1.3 m. As shown in Fig. 13a, the circle model is used to fit the local deformation of the trunk, which will cause a measurement error. Hence, we use the DBH model proposed in section II-B-4 to fit the point clouds. Fig. 13b indicates that our algorithm is closer to the real point cloud. Fig. 13c and d show the ginkgo tree fitting result. Fig. 13c shows the result by using the circle model, but there is obviously a part of the convex point cloud that is not on the circle. For these convex point clouds, we further use cubic B-spline curves to fit them as shown in Fig. 13d. Fig. 13e and f show the fitting result of the liriodendron tree.
As shown in Fig. 13e, there are some concave point clouds that do not belong to the circle. For these concave point clouds, we further use a cubic B-spline curve to fit them, as shown in Fig. 13f. The results indicate that the proposed model is more relevant to the actual points and smaller error value than the circle model.
To better analyse the performance of the proposed DBH model, the accuracy of the two fitting methods are analysed and discussed for the 93 sample trunks, including 34 ginkgo trees, 30 liriodendron trees and 29 poplar trees. Then, two comparison experiments were carried out: one experiment compared the fitting DBH value of trees from the proposed method with the actual manual measurement using callipers, and the other experiment compared the fitting DBH value of trees from the circle model fitting method with the actual manual measurement, as illustrated in Fig. 14.  where the red circles ( ) represent the measured DBH using this proposed method. The corresponding brown line is a fitting function between the measured DBH and the actual DBH; the slope is 0.976, the offset is 0.412, the RMSE is 0.532 cm, the R-square value is 0.998 and the fitting function is y = 0.976 * x + 0.412. The blue points ( ) and the corresponding purple line represent the measured DBH and actual DBH using the circle model fitting function; the slope is 0.958, the offset is 0.498, the RMSE is 0.754 cm, the R-square is 0.996 and the fitting function is y = 0.958 * x + 0.498. The grey line function equation is Y = X as a reference curve. According to the experimental data analysis, the combined circle model based on ADE and cubic B-spline method is better than the circle model fitting and the RMSE value decreased by 29.4%. The correlation between the fitted value and the actual value is checked by the Chi-Square test and MRE%. The Chi-Square test value of the circle model and cubic B-spline methods is 0.452, which is lower than the Chi-Square test value of 0.926 from the circle model. The MRE% of the circle model and B-spline method is 1.871% lower than that of the circle model, whose MRE% value is 2.706%. Finally, the F-test result p=0.99 shows that there is no significant difference (p>0.05) between the actual and measured values based on the proposed method under a confidence interval of α = 0.05. These results indicate that MRE of the proposed model is within ±5%, satisfying the accuracy standard specified by the technical indicators. And the proposed can be used to effectively extract DBH more than the circle model method.

2) CBH MEASUREMENT PERFORMANCE
In the process of CBH estimation, a comparison experimental analysis was implemented between the proposed method of CBH measurement and the voxel projection method [61]. For the same 93 samples in section III-B-1, the accuracies of the two methods are analysed and compared with the actual CBH measured by manual measurement with a tower ruler. Fig. 15 presents the results of the CBH measurements, where the red circles ( ) represent the measured CBH using this proposed method. The corresponding brown line is a fitting function between the measured CBH and the actual CBH; the slope is 0.989, the offset is 0.126, the RMSE is 0.092 m, the R-square value is 0.993, and the fitting function is y = 0.989 * x + 0.126. The blue points ( ) and the corresponding purple line represent the measured CBH and the fitting function using the voxel projection method [38]; the slope is 1.011, the offset is -0.045, the RMSE is 0.109 m, the R-square is 0.989 and the fitting function is y = 1.011 * x−0.045. The grey line function equation is Y = X as a reference curve. According to the experimental data analysis, the multiview voxel projection method is better than the voxel projection method [61] in extracting CBH and the RMSE value decreased by 15.6%. The correlation between the fitted value and the actual value is checked by the Chi-Square test and MRE%. The Chi-Square test value and MRE% of the multiview voxel projection method are 0.120 and 2.720%, which are both 0.150 and 5.013% lower, respectively, than the values from the voxel projection method [38]. Finally, the F-test result p=0.97 shows that there is no significant difference (p>0.05) between the actual and measured values by the proposed method under a confidence interval of α = 0.05. These results indicate that MRE of the proposed method is within ±5%, satisfying the accuracy standard specified by the technical indicators. And the proposed method can be used to extract CBH more effectively than the traditional method and build more accurate tree models.

3) TREE HEIGHT MEASUREMENT PERFORMANCE
In the process of tree height estimation, a comparison experimental analysis was implemented between the intersection method of tree height measurement and laser altimetry method [13]. For the same 93 samples in section III-B-1, the accuracies of the two methods are analysed and compared with the actual tree height measured by manual measurement with hand held altimeter. Fig. 16 presents the results of the tree height measurements, where the red circles ( ) represents the measured tree height using the proposed method. The corresponding brown line is a fitting function between the measured tree height and the actual tree height measured by manual measurement using a hand held altimeter, the slope is 0.991, the offset is 0.012, the RMSE is 0.327 m, the R-square value is 0.996, and the fitting function is y = 0.991 * x + 0.012. The blue points ( ) and the corresponding purple line represent the measured tree height and fitting function using the laser measuring instrument method; the slope is 1.039, the offset is -0.417, the RMSE is 0.479 m, the R-square is 0.989 and the fitting function is y = 1.039 * x − 0.417. The grey line function equation is Y = X as a reference curve. According to the experimental data analysis, the proposed method is better than the laser measuring instrument method [13] at extracting tree height and the RMSE value decreased by 31.7%. The correlation between the fitted value and the actual value is checked by the Chi-Square test and MRE%. The Chi-Square test value and MRE% of the intersection method are 0.142 and 1.075%, which are both 0.564 and 3.011% lower, respectively, than those of the laser measurement method [13]. Finally, the F-test result p=0.99 shows that there is no significant difference (p>0.05) between actual and measured values by the proposed method under a confidence interval of α = 0.05. These results indicate that MRE of the proposed method is within ±5%, satisfying the accuracy standard specified by the technical indicators. And the proposed can be used to extract tree height more effectively than the traditional method and the laser measurement method.

IV. DISCUSSION AND CONCLUSIONS
A. DISCUSSION

1) COMPARISON WITH EXISTING METHODS
In this paper, a new method for extracting the 3D point cloud parameters of trunks was proposed. Different from the method in which the tree morphological parameters are collected in a single way [61], our trunk parameter extraction method is combined with multiple methods, which is a suitable method for each parameter. The method mainly includes improvements in tree height, DBH and CBH. First, we proposed an intersection method to extract tree height instead of laser altimeters [55]. The laser altimeter cannot handle the errors caused by the tilt when it encounters tree growth. By finding the intersection of the trunk centreline and the fitted ground plane, we can more accurately extract the lowest point of the tree without being disturbed by tree occlusion. The experimental results in Section III-B-3 show that the intersection method reduced the RMSE by 31.7% compared to the laser altimeter measurement. The second is that we use a curve compensation instead of a single circle fit. Monte Carlo, optimal circle [45] and the Random Hough Transform [54] have been used to fit the circle. However, the circle fitting encounters a noncircular diameter, and it is impossible to further iteratively improve the fitting. We used the curve to fit the abnormal points to further improve the fitting of the DBH curve by the B-spline method. This kind of mixed curve composed of the B-spline curve and circle fitting curve replaces the original single circle fitting curve, which improves the accuracy and robustness of the fitting. The experimental results in Section III-B-1 show that the intersection method reduced the RMSE by 29.4% compared to the circle model fitting measurement. Third, as discussed before, the voxel projection method is not sufficiently precise, especially for clusters of interference point clouds on the trunk. Using the fixed threshold [62] to extract CBH and Lorey's mean as a CBH estimator [59], their applicability is poor for a single tree. The proposed method performs a multiangle voxel projection, and the CBH point is found by adaptive threshold segmentation. Multi-angle voxel projection improves the ability to interfere with interference points compared to voxel projection. Adaptive threshold segmentation improves the accuracy of CBH recognition. This method effectively reduces the influence of interference points on CBH parameter extraction and improves the accuracy and robustness when finding CBH. The experimental results in Section III-B-2 show that the intersection method reduced the RMSE by 15.6% compared to the voxel projection method [61].
Moreover, to verify the validity and accuracy of the parameters of the trunk model, we used Unity3d software to build visual models of the ginkgo, liriodendron, and poplar trees. Compared with the step-by-step fitting of the skeleton method extracted by Du et al. [64], we generated tree models by simulating the extracted parameters. When the branches are blocked, the skeleton method is easily affected by the blocking of the leaves, which makes model fitting impossible. Through the parameter model, we can fit the invisible part of the tree model by known parameters and trends, thereby improving the efficiency and applicability of model building. The results showed that the 3D trunk is near real trunk shape.

2) RECOMMENDATIONS
The results show the applicability of the proposed method. The proposed method can be applied to the visual system for trees to increase the accuracy of the visualization model and biomass estimation. However, there are two problems. One problem is that this method is not suitable for young trees. Therefore, for other point cloud data, analyses and comparisons should be performed to ensure the extraction effect. The other problem is that when the point cloud is too dense or there are too many interference points, the B-spline fitting will be disordered and has a certain impact on the parameter extraction accuracy for fitting DBH.

3) FUTURE WORK
Because the current scheme fails to obtain complete tree parameters, one further extension should take tree crown information into consideration. Additionally, in the process of searching for samples, many woodland trees were seriously occluded by overlap and it was difficult to locate and extract the trees. In future work, we will research the segmentation method for the point cloud of overlapping tree trunks and analyse the parameters of tree crowns to perform a complete tree morphological model reconstruction to further improve biomass estimations. Simultaneously we will continue to carry out research on visual forests using Unity3D for forest visualization.

B. CONCLUSIONS
This paper presented an effective extraction method for trunk parameters and model establishment based on multistation laser scanning point cloud data. Compared with other nondestructive forestry measurement methods, this study focuses on improving the accuracy of parameter measurement and determining the shape of individual trunks. First, We transformed the point cloud in 3D space to correct the tree point cloud. Then, We consider the individual differences of each tree to measure the tree height. The intersection of centerline and ground plane was proposed to determine the lowest point, and the local extreme value method to determine the highest point. Subsequently, a new fitting model by combining the circle based on ADE model and cubic B-spline model was proposed to more accurately extracted the DBH. The multidirection voxel projection method with adaptive threshold segmentation was constructed for CBH extraction, which improved the accuracy of parameter extraction by segmenting the viewing angle. Finally, we built a tree trunk taper model to assist in building the tree model. Based on the three sample cases, we carried out many comparison experiments and analyses of trees of different sizes and different kinds. Especially for the measurement of DBH, we performed a comparative analysis between the circle fitting method and the proposed method. The proposed method considered not only circular fitting but also irregular curve fitting. Moreover, we employed the adaptive threshold segmentation method and Multi-angle projection method rather than the fixed threshold to find voxel projection threshold points. The accuracy of parameter extraction was improved by adding judgment conditions and overall-local relationships. We compared the parameters of the trunk measured using the proposed method and some existing methods with the actual parameters obtained by manual measurement. The results showed that the RMSE value of the DBH, CBH and tree height decreased by 29.4%, 15.6% and 31.7%, respectively. According to the trunk models and measurement parameters, we tried to reconstruct the 3D trunk model. The results showed that the 3D trunk is near the real trunk shape. The results of the comparison experiments verified that multistation TLS extracted accurate point cloud data; The trunk model established in this paper can accurately extract the geometric parameters of the trunk, which can be applied to forest visualization and biomass estimation.
HUAIQING ZHANG received the B.S. degree from Central South Forestry College, in 1996, and the Ph.D. degree in forest management from the Chinese Academy of Forestry, in 2001. From 2001 to 2020, he worked with the Institute of Resource Information, Chinese Academy of Forestry, where he has been a Researcher, since 2010. His main research interests include forestry information technology, virtual reality, and 3-D visualization.
PINGPING LI received the B.S. and M.S. degrees in agronomy from Zhejiang Agricultural University, China, in 1982 and 1985, respectively, and the Ph.D. degree in agronomy from Nanjing Agricultural University, China, in 1995. She was a Teacher with the Department of Agriculture, Nanjing Agricultural University, from 1985 to 1997. She worked as a Research Professor with the School of Mechanical Engineering, Jiangsu University, from 1997 to 2011. Since 2011, she has been working with the College of Biology and Environment, Nanjing Forestry University, China, as a Professor. Her main research interests include computer vision, agricultural information engineering, and facility horticulture. VOLUME 8, 2020