Abstract
Tree skeleton could describe the shape and topological structure of a tree, which are useful to forest researchers. Terrestrial laser scanner (TLS) can scan trees with high accuracy and speed to acquire the point cloud data, which could be used to extract tree skeletons. An adaptive extracting method of tree skeleton based on the point cloud data of TLS was proposed in this paper. The point cloud data were segmented by artificial filtration and -means clustering, and the point cloud data of trunk and branches remained to extract skeleton. Then the skeleton nodes were calculated by using breadth first search (BFS) method, quantifying method, and clustering method. Based on their connectivity, the skeleton nodes were connected to generate the tree skeleton, which would be smoothed by using Laplace smoothing method. In this paper, the point cloud data of a toona tree and peach tree were used to test the proposed method and for comparing the proposed method with the shortest path method to illustrate the robustness and superiority of the method. The experimental results showed that the shape of tree skeleton extracted was consistent with the real tree, which showed the method proposed in the paper is effective and feasible.
1. Introduction
As we know, the laser has advantages of good monochromaticity and strong direction. With the rapid development of laser measurement technology, there have been some products such as airborne laser scanning, mobile laser scanning, and terrestrial laser scanning.
In recent years, TLS has been more mature with the development of computer technology and digital image processing technology [1–5]. TLS could obtain the three-dimensional information of the object by measuring many sample points of the objects [6–8].
More and more researchers have applied TLS to forest inventory and tree structural parameter extraction [9–14]. These characteristics of trees, such as the diameter of breast height (DBH), tree height, crown width, and branch structure, are important to the forest inventory. TLS is characterized by high precision and large data volume, which is very suitable for data acquisition. It can effectively avoid waste of labor and time consumption and improve the accuracy of the collected information. Based on these point cloud data obtained by TLS, we can measure many parameters of trees, such as tree DBH [15], tree height, biomass, and extraction of tree skeleton [16–19].
The tree skeleton is a concise description of the three-dimensional structure of trees, which can describe the shape and topological structure of a tree [20–23], and tree skeleton contains information about the geometry of the trees, branching patterns, and branching angles. And tree skeleton can also reflect the growth process and characteristics of trees [24]. The tree skeleton is the basis for studying the morphological structure of the trees. How to effectively and accurately obtain the tree skeleton is the key factor to analysis of the morphology of trees.
Some algorithms were designed to extract the tree skeleton from point cloud data [25–31]. First of all, the skeleton nodes which can represent trunk and main branches would be extracted. And then skeleton nodes could be connected according to the tree structure. Finally, the tree skeleton would be generated.
In order to realize this idea, Xu et al. had presented a method about extracting tree skeleton from point cloud data based on the shortest path algorithm and clustering algorithm [32].
Yan et al. had used -means clustering method and cylinder detection method to get the tree skeleton [33]. In Delagrange and Rochon’s study [34], the skeleton is connected according to minimum spanning tree method. Wang et al. had presented a global optimization method based on perceptual structure [35, 36]. Bucksch and Lindenbergh had proposed a method of extracting tree skeleton from point cloud data based on octree [37].
In the paper, we proposed a method to extract the tree skeleton based on breadth first search (BFS) method from the point cloud data acquired by a single scanning of TLS. Firstly, leaves points are separated from the scanning point cloud artificially, and ground points are separated from the point cloud by using -means clustering based on normal vector of points [38]. Then an undirected graph is built according to the neighborhood relationship and every point is labeled by using BFS. After quantization operation, some points would be clustered to a block, followed by calculating the skeleton nodes in each block. However, some points that did not belong to the trunk or a significant branch (leaves, twigs, and noise) would have effect on the skeleton generation. Two ways were used to reduce these points. First one is called threshold filter, which ignores these points in a block whose size is less than the threshold. Another one is called proportion filter, which ignores these points in a block whose size is less than the proportion of its father block size. If two points belong to an edge which were in different blocks, these blocks have the connection. According to the connection of blocks, the skeleton nodes belonging to these blocks are connected. Finally, tree skeleton would be generated after connecting skeleton nodes. And the skeleton needs to be smoothed by using the Laplace method.
2. Method
Tree skeleton could describe many essential characteristics of trees, such as tree height, branch structure, and distribution. The aim of skeleton extraction is to get the morphological structure and branch structure of trees.
In our experiments, tree skeleton was built according to the following process also described in Figure 1. The proposed method is composed of three phases. Phase 1 is excluding the ground points from point cloud data by -means method. Firstly, the point cloud data of a tree was extracted from the point cloud data set of whole scene. Then the point cloud data belonging to ground and foliage was removed, and the point cloud data belonging to trunk and branches remained, which was called working data in the following context. The second phase is extracting the skeleton nodes from tree point cloud. An undirected graph was established by using KD-tree method, followed by clustering the blocks of points by quantizing. It could automatically judge whether filtering process needs to be taken by estimating the size of each block and deciding to choose the threshold filter or proportion filter. The tree nodes were found by using the clustering method based on the working data. Finally the third phase is connecting skeleton nodes to construct the tree skeleton according to graph theory and smoothing the tree skeleton. Then tree skeleton would be generated.
2.1. Data Preparation
In our experiments, TLS was used to scan the trees. Meanwhile the acquired point cloud data contains not only the trees but also the ground and the background. In order to simplify the analysis, the point cloud data of a single tree was manually extracted from the whole point cloud data in the study. And the point cloud data of leaves and ground should be removed because they have nothing to do with the extraction of tree skeleton.
The distances between the branches and the leaves are relatively small, which will affect the skeleton extraction. Therefore the point cloud data of leaves should be removed before extracting the skeleton. In the experiment, the removal was carried out by using the intensity value obtained while scanning. The ground and the trunk end were connected and it was not easy to directly separate them because of the complex ground topography and tree growth direction. In the experiment, the normal vector method and -means clustering method were used to separate the point cloud data of ground. According to the tree growth characteristics, the trunk and the ground generally form a larger angle. Therefore the normal vector of the ground and the normal vector of the trunk would be different in the direction. Normally, the direction of the trunk is close to the vertical direction, and the normal vector of trunk is close to the horizontal direction. And the normal vector of ground is close to the vertical direction. Therefore, the direction of normal vector could be used to distinguish the trunk and the ground.
In clustering, we had chosen the -means method, which is a nonsupervised clustering algorithm for multidimension data and does not require prior knowledge. We chose the 0 degrees and 90 degrees as cluster center. Then according to the distance, which is not the Euclidean distance, but an abstract link, here it refers to the normal vector value of point. If its value of the normal vector is close to the value of the selected ground point, the point is classified as the ground point cloud and the cluster center is recalculated. If its value of the normal vector is close to the value of tree points, the point is classified into the tree point cloud, and the cluster center was recalculated. After traversing all points, the point cloud data of ground would be separated from those of trunk end.
The objective function of -means in the proposed method can describe asIn the above equations, is the minimum sum of the distances between each point belonging to and , which means the result of clustering is the best choice. represents the normal direction angle of a point; is the size of point cloud. And is the number of clusters. In this paper, the value of is 2 which means ground and tree truck. is the center of one cluster. means that if a point belongs to , the value of is 1. Finally, the -means method would find the values of and , which can make the objective function converge to local minimum.
2.2. Generating Skeleton Nodes
Firstly, an undirected graph would be created based on the distance between scanning points. Each point was connected to all its neighbor points. The distances between points are different due to the different accuracy, different range, resolution, and other scanning parameters. It will cost a long searching time because of the large amount of point cloud data. We used the KD-tree which was used wildly in computing graphics to implement the nearest neighbor query and establish the relationship between points.
KD-tree is a binary tree search structure for -dimensional data; here takes 3, which refers to the three-dimensional data. The KD-tree method starts a binary searching from the root node in different dimensions. The nearest neighbors of a point were obtained by using the nearest neighbor query, and the two-way edges between the point and each neighbor were established.
The connected undirected graph could be defined as . Some vertices in were selected as root set according to (2). Each vertex in was a point of point cloud data. The root set included all points whose distance from the lowest point is less than in the -axis direction and marked by 1. is the coordinate of a point. is a point in root set. is the lowest point in the point cloud data. is a threshold and was set to 0.02 m in the experiments. As shown in Figure 2, all points that meet the requirements were classified as root points. And other points were traversed and numbered according to the relationship of undirected graphs.
In Figure 2, the value of vertex in root set was 1, and the values of other vertices in were calculated from root set by BFS and the following equation: is the father node of in BFS process. is the value of a vertex. is the value of father node of .
The calculating process of the values was demonstrated in Figure 3. At the beginning of the method, the values of blue vertices which belonged to root set were initialized to 1. Then the values of green vertices which were adjacent to vertices in root set were calculated. Next, the values of yellow vertices were calculated in the same way. And the points will be marked only once in the whole process.
(a) The points before drawing
(b) The points after drawing
When all values of vertices have been determined, the values should be quantized and mapped in a range. In our study, we set the range as 0 to 60. In Xu’s research [32], all the values would be linearly mapped in the range of 0 to 60 shown in Figure 4. However, the linear mapping mentioned above could cause the loss of details of tree branches, because the branches far away from the root are very thin and short. We proposed a nonlinear mapping method to reduce the loss of the small branches details that was shown in Figure 4. All the values would be quantized as is the quantized value of point , is the initial value of point , is the maximum value in the point cloud data, and brackets denote getting the largest integer no more than the value in the brackets.
Due to different types and structure of trees, the quantitative standard could be varied. In our experiments, all the points would be quantized into 60 intervals with different values. And each interval may contain parts of different branches, so we need to determine whether they belonged to the same branch.
Depth-first search (DFS) was used to determine which points belong to the same block. When searching the neighboring point for a point, the neighboring point () of the point and the point () have the same quantization value; then the two points are considered to belong to the same block. The points belonging to the same cluster should meet and are two points in the same interval. refers to determining whether and have connectivity. If , and are connected. If , and are not connected.
Figure 5 shows the result of clustering. As we can see, block 1 and block 2 belong to the same interval; however they were two different clusters without connecting relationship.
In the practical processing, some outliers would interfere with the generation of tree skeleton nodes. The outliers include the noisy points, the points of leaves, and small braches which often contain a few points. So a threshold value was proposed for the number of points to remove these outliers. And the setting of the threshold is related to the characteristics of the trees and the scanning parameters. We defined some point cloud clusters whose sizes were less than 30 or the ratio of size between the son block and the father block was too large which will be deleted to prevent the occurrence of abnormal skeleton.
The center point of each cluster was calculated according to (6). The center points which are centers of every part of tree were referred to as skeleton nodes. is a center point. is the size of a cluster; is the point in the same cluster.
The skeleton nodes always were calculated by using three-dimensional Euclidean distance. We can also use point cloud data intensity values to calculate the coordinates of the skeleton node; the intensity values of points which are in the middle of the trunk are usually greater than the intensity values of points on both sides in the point cloud acquired by the laser scanning, and the skeleton node usually exists in the trunk center. So intensity values can be used to calculate the coordinates of skeleton nodes.
2.3. Skeleton Production
The generation of tree skeleton is to connect the skeleton nodes. We use merge-find operation to connect skeleton nodes. The merge-find operations were two types of operation: find operations and union operations. Find operation can determine the block that a point belongs to. Union operation can merge two blocks into one subset. We visited each edge in the graph . Suppose two vertices are related to an edge which are labeled as and . and are the center points which belong to the clusters where and located. If and Find ()! = Find(), then we connect and , that is, connecting the skeleton nodes which belonged to two different clusters that have the connectivity which was shown in Figure 6. A skeleton of tree would be produced when all the skeleton nodes were connected.
However there were still some ghosts of skeleton produced, because some of the leaves and noise points may be recognized as branches. In order to filter out these skeleton nodes which were actually recognized incorrectly, and will not be connected when the distance of and is less than 0.2.
2.4. Skeleton Smoothing
It would be found that the tree skeleton established by the above method was not smooth. In the paper, we adopt the smoothing processing to increase the smoothness of the whole skeleton. On the premise of ensuring that the tree shape does not occur with distortion, we modified the three-dimensional coordinates of the skeleton nodes except bifurcation points. are the coordinates of the new point NEW_A. are the coordinates of the point . are the coordinates of the point . are the coordinates of the point .
NEW_A is the changed current skeleton node, is the current skeleton node, is the parent node of , and is the child node of . The values of NEW_ are the mean of their coordinates . The NEW_A’s three-dimensional coordinates would replace the current skeleton node A. The processing is shown in Figure 7.
The blue nodes are the original skeleton nodes, where the triangle node refers to the root node, the circular nodes are the branch nodes, and the red nodes are the new skeleton nodes that have been smoothed. It can be seen that some blue circular nodes’ linear relationship with other nodes is very poor, which will affect the skeleton reconstruction. Therefore it is necessary for these convex points to do smoothing processing.
3. Experimental Results
In the paper, the point cloud data of a toona tree and peach tree were used to test the proposed method. We distinguished the ground and tree trunk according to their normal vectors. And the point cloud data of ground were separated successfully by using -means clustering method. The experimental data and results of toona tree and peach tree are shown in Table 1.
The skeleton model generated by using the proposed method had good effect. Compared with the original point cloud, it had high reduction and was consistent; otherwise, we had got more details of branches in comparison with shortest path method. The experimental results were shown in Figures 8 and 9.
(a)
(b)
(c)
(d)
(e)
(f)
Figure 8 shows the skeleton extraction result of toona tree and peach tree. It could prove that the method has a certain degree of robustness. And it can be seen that the shape of skeletons is consistent with the shape of tree point cloud.
Next we mainly choose the toona tree as the example to illustrate the results of the proposed method. The parameters and results of the proposed method and the shortest path method for toona tree are shown in Table 2. Label a and Label b mark experimental information of the proposed method and the shortest path method, respectively.
Figure 9 shows the process results of extracting tree skeleton with the proposed method and the shortest path method. The point cloud data of the same toona tree was treated in two ways. It can be seen that the proposed method could extract the tree skeleton with many details by enlarging the branches of the skeleton extraction.
We should consider some factors when we use the above method. These factors decided the final tree skeleton model, which are the search radius (S-R), quantization scale (Q-S), threshold filter (T-F), and proportion filter (P-F).
Search radius (S-R) determines the size of undirected graph. It determines whether a point would be searched and used to build the undirected graph. Quantization scale (Q-S) determines the number of skeleton nodes to be generated. The larger Q-S will result in more skeleton nodes and therefore more details of branches.
Threshold filter (T-F) determines whether some clusters will be ignored. In the point cloud data of tree, there are often some nonskeleton points, such as leaves and noise points. These nonskeleton points maybe judged as skeleton points, which will interfere with the skeleton generation. In the process of clustering, these points will be put into one cluster. There is often a small amount of these points in the experiments. Therefore a threshold was used to filter these points. If the number of points in the cluster is less than the threshold, these points will not participate in the generation of tree skeleton.
Proportion filter (P-F) has the same function of threshold filter (T-F). Because of the scanning accuracy, the size of the tree itself, and other factors, the number of nonskeleton points may be larger than the filter threshold. Then T-F may not have the possibility of eliminating the nonskeleton points completely. In this case, P-F could be used to do further filtering. That is, when the number of points in the current cluster is less than the number of points in the father cluster block and the ratio between them is less than proportion filter threshold, then these parts of the point cloud will not be considered in the generation of skeleton.
The six configurations of the three-dimensional reconstruction of the tree skeleton are shown in Table 3, and the tree skeleton building results were shown in Figure 10. The sum of tree points (SP), the sum of tree skeleton nodes (SN), operation time (), and sum of tree skeleton edges (SE) were the four input parameters that were used to build tree skeleton.
(a)
(b)
(c)
(d)
(e)
(f)
SP refers to the amount of data obtained after the preprocessing, and SNs refer to the amount of skeleton nodes that we have calculated. Operation time () refers to the time of running the codes to build tree skeleton. SE refers to the amount of edges in the tree skeleton, that is, the standard of judging the branch details.
The results of building tree skeleton with different parameter sets (S-R, Q-S, T-F, and P-F) were shown in Figure 10. It can be seen that different parameter sets could result in different results based on the point cloud data of same tree.
The smoothing process proposed in the paper is based on Laplace smoothing. We choose those nodes that are not bifurcation nodes and adjusted the position of these nodes replaced with the average value of position coordinates of its father node and its child node and itself, so that some of the bumps can be smoothed out to prevent producing some irregular branches. And the method can increase the smoothness between ordinary nodes and the adjacent nodes to generate a more flexible and beautiful tree skeleton. The result of smoothing process was shown in Figure 11.
(a)
(b)
(c)
4. Discussions
The tree skeleton constructed by the method proposed in this paper basically restores the morphological characteristics of the trees which was shown in Figure 8.
We use different kinds of trees to verify the robustness of the method. This shows that the method has robustness by showing the skeleton extraction results of the toona tree and the peach tree.
And comparing the skeleton we had generated with the original point cloud, we can see that we better restore the morphological characteristics of the main branches of the tree. That is to say, the visual quality of our extraction effect is good.
And the method could retain more branch details of tree by comparison with the shortest path method seen in Figure 9.
In the process of building tree skeleton, the roles of the parameters are very important. It will get different experimental results with different parameter configurations that can be seen in Table 3 and Figure 10.
Larger S-R will get a larger SP. This means that more point cloud data could be searched and used to generate tree skeleton. Figure 10(b) has a larger S-R than Figure 10(a); we can see that it has more point cloud data. Larger quantization range will result in larger SN. We could get more intervals that would generate more tree skeleton nodes. Figure 10(c) has a smaller quantization range than Figure 10(a), so it has less skeleton nodes. And taking the T-F and the P-F will reduce the effect of nonskeleton points in tree skeleton reconstruction. Figure 10(d) does not take threshold filter and proportion filter compared with Figure 10(a), so it has more nonskeleton nodes data. The function of threshold filter is deleting some nonskeleton points to prevent generating nonskeleton nodes. Figure 10(e) takes threshold filter, so it has less nonskeleton nodes. Figure 10(f) takes proportion filter; it has also less difference with Figure 10(a). The function of proportion filter is to reduce the edges which would be used in generating tree skeleton. Figure 10(f) takes proportion filter, so it has less SE than Figure 10(d). We also can see that there is less difference between Figures 10(a) and 10(f); it explains that its role is relatively larger than threshold filter.
In the process of quantification, it is very important to select the appropriate quantization range according to the setting of scanning parameters and the number and intensities of the point cloud data of trees. Under the condition of efficiency and reconstruction effect, selecting the appropriate quantization value can keep more branch structure information. Different quantization range will produce different results. The more the quantization range, the more the skeleton nodes that we got. And the time needed is longer, while a longer building time will be needed to generate tree skeleton. If the point cloud data has a simple structure, increasing the value Q-S will have little effect on getting more details.
We have two ways to remove the nonskeleton point; the first one is to determine whether the number of points in cluster block exceeds the threshold; if it is less than the threshold, all points of the block will be ignored to calculate skeleton node. Another one is to judge the ratio between the number of points in the parent cluster and child cluster; if they are less than the proportion threshold, they will also be ignored. Taking the two filter processes will produce less nonskeleton details.
When we build tree skeleton by connecting the skeleton nodes, we should not only take into account the distance factors, but also pay attention to the connectivity of the clusters, because the different clusters that belong to different branches may have the same distance to the father cluster; that is, there is no connectivity between them. They need not be connected to prevent the emergence of unnecessary ring structure.
The necessity of smoothing the skeleton is that the skeleton only roughly reflects the trees in the morphological structure; it cannot completely restore the real tree structure. We can see that the skeleton is more smooth than the original skeleton. In addition, smoothing can be manipulated several times; we can also see the tree skeleton smoothed two times is smoother that the skeleton smoothed once.
5. Conclusions
The method proposed in the paper can extract tree skeleton based on the point cloud data. In the experiment, a toona tree and a peach tree were scanned by TLS and the point cloud data of the tree was used to extract the tree skeleton. The point cloud data of ground and leaves were removed. And the point cloud data of trunk and branches were used to generate the tree skeleton. Skeletons of the toona tree and the peach tree were generated successfully by using the proposed method. The method was proved to be effective and efficient by the results.
Parameters used in the method were discussed in the paper. Selection of the parameters is related to the tree types, structures, and scanning parameters. Optimization of the parameter selection will be a part of future work. Also skeleton-based tree rebuilding will be studied in the future work.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors gratefully acknowledge the financial support from the Fundamental Research Funds for the Central Universities (no. 2015ZCQ-LY-02) and the Beijing Natural Science Foundation (no. 6164037).