Automatic and Self-Adaptive Stem Reconstruction in Landslide-Affected Forests

Abstract

Terrestrial laser scanning (TLS) is a promising technique for plot-wise acquisition of geometric attributes of forests. However, there still exists a need for TLS applications in mountain forests where tree stems’ growing directions are not vertical. This paper presents a novel method to model tree stems precisely in an alpine landslide-affected forest using TLS. Tree stems are automatically detected by a two-layer projection method. Stems are modeled by fitting a series of cylinders based on a 2D-3D random sample consensus (RANSAC)-based approach. Diameter at breast height (DBH) was manually measured in the field, and stem curves were measured from the point cloud as reference data. The results showed that all trees in the test area can be detected. The root mean square error (RMSE) of estimated DBH was 1.80 cm (5.5%). Stem curves were automatically generated and compared with reference data, as well as stem volumes. The results imply that the proposed method is able to map and model the stem curve precisely in complex forest conditions. The resulting stem parameters can be employed in single tree biomass estimation, tree growth quantification and other forest-related studies.

1. Introduction

Forest attributes, such as stem location, diameter at breast height (DBH), height, basal area, stem curve and volume, are critical to timber industry management (e.g., [1]), assessing the potential of wild fire hazard (e.g., [2]) and natural biodiversity conservation (e.g., [3,4,5]). Among them, the stem curve, which describes the diameter at any height along the stem, is of great importance [6]. Studying the stem curve is crucial not only for forest management and biometrics, but also for research on the geomorphological environment. In particular, understanding the role of trees and forests is increasingly important in high mountain areas, as sloping regions are often characterized by shallow and very slow moving landslides [7]. The so-called “drunken trees”, which means that tree stems are displaced from their vertical alignments, are caused by such soil movements [8]. In recent decades, there has been an increased occurrence of soil erosion in many alpine regions [9]. Climate change adaption and natural disasters in these regions are receiving considerable attention [10]. For example, recently-developed dynamic slope stability models (e.g., [11,12]) studied the interdependencies of different processes, including deforestation or afforestation. Detailed stem volume and biomass information at single tree level can greatly benefit and facilitate such studies [13]. Moreover, the stem curve can be used to quantify the tree growth anomaly in landslide forests (e.g., [14]). Some vital attributes can be derived, such as taper [15] and inclination angle [14].

The interest of studying the stem curve, or stem form, or taper dates back long before the 1950s [16]. There are a number of ways to measure the stem curve. Conventional forest inventory approaches for stem curve measurement are labor intensive and sometimes harmful to trees. Trees are often felled and then measured from the stump to the top level. In addition, mathematical models also are developed to predict stem curves [17]. These methods use only fixed measurements, usually DBH, total height and several measurements on the stem. The stem curve is then interpolated by mathematical functions (e.g., [18]). The relation of different components of trees is often estimated by allometric functions (e.g., [19,20,21,22]) based on DBH, tree height or other parameters of interest. However, these functions are not always valid, because they were developed from specific local morphological or climate conditions and do not provide the position at any height along the tree stem. Precise stem models would help to improve the accuracy of such estimations [23].

Terrestrial laser scanning (TLS) has been widely used in forest-related studies (e.g., [15,24,25,26,27,28]). A laser scanner acquires a high density point cloud, which enables the stem form reconstructions in a nondestructive way. Tree attributes, such as stem location, DBH, basal area and volume, can be effectively retrieved from the TLS point cloud using some automatic algorithms (e.g., [25,26,27,28,29,30,31]). Previous studies mostly focused on retrieving DBH or tree height, whereas a limited number of publications explored the applicability of modeling the stem curve (e.g., [32,33]). Circle fitting (e.g., [34]) and cylinder fitting (e.g., [35,36]) are the primary strategies often mentioned in the literature. A robust stem curve retrieval method based on cylinder fitting was described by the authors in [37], but their method focused on a forest in a flat environment. Nonetheless, most of the previous works focused on forests in flat environments or planned forests. To our best knowledge, similar approaches have not yet been applied to high mountain natural forests or shallow landslide-affected forests. Such forests are often characterized by steep terrain with a multi-layered canopy structure, including dense understory, mixed tree species and deformed stem shapes. The stem formation is often influenced by factors, such as the site fertility, spacing and light conditions, wind and landslide events [38]. Therefore, the stems are often growing in a manner deviating from the vertical direction and have irregular forms (Figure 1). Special methods have to be developed to adapt to these challenges.

The current study focuses on using TLS for automatic stem curve modeling in landslide-affected forests. The objective of this study is to expand the stem curve investigation to landslide-affected forests and to discuss the concomitant challenges. Finally, we present a novel algorithm for stem modeling in mountain forest environments.

2. Study Data

The study site is located in the federal state of Vorarlberg, Austria (47.224°N, 9.725°E). The Vorarlberg Alps features widespread translational landslides [39]. Our specific study site covers a plot of approximately 31 m × 16 m, inside a small forest located near the rupture surface of a shallow landslide. The site is characterized by steep terrain (>30°) with a multi-layered canopy structure, including dense understory, mixed forests and dead tree branches (Figure 1). Tree stems are overall anfractuous due to the effects of soil movement. The dominant tree species is Norway spruce (Picea abies (L.) H. Karst.). The study site consists of 27 trees with a DBH larger than 5 cm, whereas other small trees are neglected in this study. The density of trees is approximately 544 stems/ha. The DBH was manually measured in the field using a measuring tape. The DBH value in this study site ranges from ∼9.4 cm–∼61.2 cm, with a mean value of 32.8 cm and a standard deviation of 14.3 cm.

The TLS measurement was carried out in October 2015, using a Riegl VZ-2000 scanner (RIEGL Laser Measurement Systems, Horn, Austria). The scanner has a vertical view angle of 100° (+60°/−40°) and a full 360° horizontal view angle, with an effective measurement rate up to 400,000 points per second.

Table 1. Specifications of the Riegl VZ-2000 scanner.

Specifications	Riegl VZ-2000
Max. vertical field of view (°)	100
Max. horizontal field of view (°)	360
Accuracy (mm) at 150 m range	8
Points per second (max)	396,000
Beam divergence (mrad)	0.3
Max. resolution (°)	0.0015

summarizes the detailed specifications of the scanner. Seven scans were done according to the terrain accessibility in order to achieve a good laser scanning coverage (i.e., the percentage of a stem cross-section that is covered by the TLS point cloud) of all tree stems from different directions (Figure 2). Each tree was visible from the perspective of several scans depending on occlusions. The average TLS coverage rate for a single tree was 77.7% (Figure 2). The full coverage is defined as the complete stem (i.e., full cross-section) being scanned. Figure 2 displays the coverage rate for each tree. Reflectors were placed on the tree stems for the purpose of registration. The placements of reflectors ensured that one scanning position can be integrated with other scans. Afterwards, the seven scans were registered using Riegl’s RiSCAN PRO software (RIEGL Laser Measurement Systems, Horn, Austria). The overall registration accuracy is ±7.5 mm. Consequently, seven scans were merged, and no processing was applied to the data (e.g., removal of undergrowth, manual identification of trees), meaning that all branch and outlier points near a stem were retained.

In addition to TLS and DBH field measurements, stem curves were manually measured from the acquired point cloud using the CloudCompare software (CloudCompare 2.6.2, 2016). The diameter and the location of the stem center at the height of 0.65 m above ground were measured. The following measurement height was at 1.3 m (i.e., the height of DBH), then every one meter above until reaching the top of the stem or where no points can be identified as part of a stem. The diameter was determined by averaging the values from two directions (W-E, N-S), and the location of the stem center was the intersection.

3. Methods

The proposed approach automatically models tree stems from the point cloud. The anticipated workflow contains four procedures. First, the original point cloud is down sampled by a leveled histogram sampling method [40]. Second, the terrain points are identified and removed by a hierarchical approach [41], so that trees are assigned with normalized heights above the ground. Third, stem locations are recognized from the spatial distribution of stem points and labeled as individual groups. Finally, the stems are modeled by fitting a series of cylinders, based on a 2D-3D robust cylinder fitting strategy. In this step, the stem parameters, such as the DBH, diameter and location at various heights (i.e., stem shape) and stem volume, can be retrieved automatically. The overall workflow is given in Figure 3.

3.1. Down Sampling

Previous studies (e.g., [42]) showed that advanced sampling techniques are able to reduce point cloud processing time and retain the outcome quality simultaneously. Our original point cloud contains over 205 million points. The considerable redundancy calls for down sampling approaches. In this study, the point cloud sampling was performed with a sampling routine named leveled histogram sampling [40]. The routine was selected based on its speed and effectiveness. The leveled histogram sampling aims to collect scanned points evenly from a distribution describing the point cloud with a selected metric. Here, the 3D distance of the points from the scanner was used. The algorithm requires two control parameters, the sample size and histogram bin width. For bins with a high point number, the point selection is performed in a uniform fashion, as individual points are not assumed to have critical information in them. Bins with a low number of points are included as such. The total number of selected points in all bins will correspond with the required sample size in the end. In this study, 100 histogram bins are selected. These corresponded approximately to 0.40 m bin widths on average for the tested scans. A 7.5% sampling rate was chosen based on practical tests [43].

3.2. Terrain Model Derivation

The determination of the digital terrain model (DTM) was done by a hierarchical approach [44]. Starting from a thinned point cloud, the lowest points within 4 m × 4 m raster cells were used for robust moving planes’ interpolation. For filling the gaps in the derived model, a triangulated model was used, which was derived from the lowest points within the 4 m × 4 m raster cells. The determined elevation model is used for normalizing the elevations of the point cloud.

The robust filtering approach was originally developed for DTM extraction from airborne laser scanning data (e.g., [41]). One main assumption is that terrain points are located at the lowest. Thus, starting with a rough surface model, it can be assumed that it is more likely that a point below the “rough surface model” belongs to a terrain point than points above the “rough surface model”. Based on this assumption, a weight function can be defined (details can be found in [41]). This procedure is done in an iterative way starting with a coarse raster and ending with a fine one. Therefore, we started first with a cell size of 4 m × 4 m and did the refinement with 2 m × 2 m, 1 m × 1 m and, finally, 0.2 m. The entire workflow is implemented as a batch script using the software Orientation and Processing of Airborne Laser Scanning data (OPALS) developed at the Vienna University of Technology [45] and does not need any user interactions.

3.3. Stem Location Detection

There are a number of methods to delineate single tree stems in the TLS point cloud from forest plots. Hough transform and circle searching were explored in [26]; points between 1 and 2 m above ground were selected to find tree positions. The authors in [28] used the spatial clustering method for points lower than 4 m above the ground. Points were projected to the horizontal plane and grouped by their spatial distances. The work in [46] searched trees in the TLS range images. Salient features and the normal vector method were also used as an operative approach to identify points from different tree components based on their spatial distribution patterns [30,47].

In this study, we propose a two-layer projection approach to determine the location of each tree, instead of directly delineating stem points, due to the difficulty or specific challenges in identifying stem points in landslide-affected forests. The severe occlusion causes trouble for the circle searching and range image methods. Local foliage aggregations may resemble the stems, and the tree stem shapes are irregular (e.g., Figure 1 and Figure 4); thus, spatial pattern recognition with normal vectors becomes often impractical. Hence, we define the first layer as a full layer, in which all points from the laser scans after terrain removal are retained. The second layer is a sub-layer, which only includes points between 2 m and 4 m above the ground. Such a height range is useful to diminish the disturbance from dense undergrowth and tree crowns. Points in the sub-layer were down sampled by grid boxes with a 2-cm side length in order to make the point density uniform for all trees, which is crucial for calculating the density map and normal vectors. Points within each grid box were merged by averaging their locations.

One can define a point cloud as $P=\{x_i,y_i,z_i,\parallel i=1..n\}$; $P\subset \mathbb{R}$, $x_i,y_i,z_i$ are the coordinates of a point $p_i$ in P. The normal vector n at a point $p_i\in P$ can be estimated as the eigenvector to the smallest eigenvalue ${\lambda }_i$ of the covariance matrix given by:

$$Co{v}_p=\frac{{\sum }_{i=1}^K(p_i-\overline{p}){(p_i-\overline{p})}^T}{K},$$

(1)

where $\overline{p}$ is the barycenter of the K points. A low value of the z-component of n ($n_z$) (z-normal) means that a point is approximately on a vertical structure [30]. In this study, the normal vectors are calculated for the sub-layer. The z-normal is estimated from the covariance matrix. A point then can be represented by $p_i(x_i,y_i,z_i,n_{z_i})$. The sub-layer is then projected onto a horizontal plane of 2 cm × 2 cm cells. The size of the cell is determined to separate stems. For each cell, the number of corresponding points N (i.e., the density map) and the median value of z-normals (i.e., the z-normal map) for all points in this cell are determined. We use the median value so that the impacts from some outliers and poorly calculated normal vectors can be minimized. For a cell belonging to a stem, a high value of N together with a low value of the average z-normal mean that this represents a vertical planar surface. Nonetheless, in mountain forests, locally aggregative dense foliage could also result in garbled characters with stems (Figure 4). Therefore, it is difficult to separate stem points with foliage points merely from projection density or normal vectors. Moreover, the determination of the the thresholds for density and normal vectors (or fatness) is unclear, because they depend on the overall laser point density and local field conditions. The authors in [28,30] did not mention how they selected their thresholds. In light of that, we further develop a quantitative method to combine the density map and the z-normal map and determine their thresholds for identifying tree locations.

Figure 5 shows the relation between the density and z-normal for points in Figure 4. In general, such a relation curve will have an “L” shape, because the cells represent foliage that will normally have lower densities with high z-normal values. Therefore, the curve drops quickly in the beginning. When the density increases, the z-normal value will eventually become small enough until reaching a stable state, where the density is high enough so that the cells are dominated by stems. Such a curve can be well approximated by the Gaussian curve (bell-shape). A second order Gaussian curve:

$$f\left(x\right)=a_1{e}^{-}\frac{{(x-b_1)}^2}{2c_1^2}+a_2{e}^{-}\frac{{(x-b_2)}^2}{2c_2^2}$$

(2)

is accurate enough. The thresholds are the values where the curve becomes flat (i.e, the maximum curvature). Consequently, the density and z-normal thresholds can be determined, and all cells having a larger density and a smaller z-normal value will be identified as stems.

The corresponding points are further grouped by the density-based spatial clustering of applications with noise (DBSCAN) algorithm [48] on x-y coordinates. This algorithm finds all connected points from arbitrary shapes within 2 cm (radius of the neighborhood) and clusters them as groups. It does not require the number of clusters as an input, different from other methods (e.g., the k-means algorithm). Further, the natural occlusions may separate the points belonging to the same stem to different groups. Thus, the groups that are distributed closer than 40 cm are merged as one.

Notably, in theory, the aforementioned DBSCAN detection method only works well for vertical stems or trees with similar inclination angles; while in complex forests where some trees are vertical and some are not, such an assumption sometimes holds wrong. Our method assumes that at least part of the stem points are recognized and grouped (i.e., the parts with similar inclination angles). The x-y extent of every group is enlarged by a certain factor (e.g., a factor of 1.5 in this study) for every group, so that the surrounding points from the full layer can be retained if only parts of the stem are detected. Moreover, the group enlargement is important as a significant portion of trees in landslide-affected forests is typically inclined. This is also the reason we only estimate the location of every stem, instead of delineating stem points directly from the TLS point cloud. The range extension will probably also incorporate points from branches, foliage and even the understory. Nonetheless, the proposed stem modeling technique is able to mitigate such effects in a subsequent processing step.

3.4. Stem Reconstruction

A tree stem can be modeled by fitting a series of cylinders (e.g., [37]) or circles at various heights (e.g., [26]). The latter seems problematic in mountain forests because the severe occlusions could block the laser from reaching the stem at many heights. For the cylinder fitting approach, the challenges arise from the dense understory and foliage and irregular stem shapes in our study. Here, we propose a 2D-3D robust cylinder fitting and self-adaptive growing strategy to reconstruct the whole stems. A circular cylinder can be parameterized as:

$$\left|\right|(P-Q)\times a\left|\right|-r=0,$$

(3)

where $P=(x_i,y_i,z_i)$ is a point fulfilling this equation; thus, it is on the cylinder surface. $Q={(x_a,y_a,z_a)}^T$ is a point on the cylinder axis, and r is the radius. a denotes the direction of the cylinder axis and has a unit length. The infinite cylinder can be defined by 5 parameters, radius r and 4 more to define the axis.

The solution of the nonlinear problem is to minimize the residual ${\sum }_{i=1}^n{v}_i^2$ where $v_i=\left|\right|(P-Q)\times a\left|\right|-r$. It can be solved by the Gauss–Newton or Levenberg–Marquardt schemes. However, the performance of iterative algorithms heavily depends on the choice of the initial guess. Besides, there is always a chance that they would be trapped in a local minimum [49]. Further, standard least square, however, is not able to eliminate the outliers originating from branches, leaves and wind. The authors in [37] used a robust method with Tukey’s estimator. We solve this problem by introducing a 2D-3D approach. Outliers are preliminarily first removed by a 2D region growing method on the density map. Robust estimation of the initial guess values of a cylinder is also determined on 2D by a random sample consensus (RANSAC)-based [50] circle inscribing approach. Finally, the cylinder is fitted and extended in 3D.

3.4.1. Starting Cylinder

For each tree, we start from a 20-cm section between 2 m and 4 m above the ground, which is consistent with the stem location detection procedure. A section is selected if it contains more points than other sections. The section with the most points is more plausible to achieve a good fitting, as possible occlusion may present between 2 m and 4 m above the ground. The critical issue is then how to eliminate the effects of outliers. Moreover, the determination of the initial guess for the cylinder axis position and radius is of great importance for solving the nonlinear least square fitting. Robust estimation should be developed at this stage.

A region growing algorithm is firstly applied to remove some conspicuous outliers for each tree (Figure 6) on the density map. We only consider the cells that have valid densities. The cell that has the largest density is defined as the starting seed. The starting seed is assigned with a group. A neighbor cell is added to the same group if its density is larger than 20% of the seed. After all qualified neighbor cells are added to the seed group, a new seed is selected (i.e., the cell that has the largest density of the rest of the cells). The procedure is continued until all cells are visited. The resultant group that contains the most points is regarded as the dominant structure, thus kept as the stem or at least part of the stem. Figure 6 shows an example of two traversed groups for a stem section.

Some branches may also comprise a large amount of points and have similar density with stems; thus, they cannot be eliminated from the aforementioned region growing method. We further introduce an RANSAC-based approach to achieve a robust estimation of an initial guess for the cylinder axis position and its radius. The approach is done in 2D by inscribing a circle to the stem cross-section and then transform the circle to 3D.

The RANSAC algorithm iteratively estimates the parameters of a mathematical model by randomly sampling the data. Its effective iteration time M is defined by:

$$M=\frac{log(1-p)}{log(1-{(1-\lambda )}^m)}$$

(4)

where p is the probability to achieve a good model; m is the minimum amount of data required for determining the model. For circle determination, m equals 3; λ denotes the estimated percentage of error in the data, which is determined to be 40% in this study because of the complex field conditions. For each iteration S, a circle is determined by three randomly-selected points. Next, the distances d from all other points to the circle are computed by:

$$d_i=\sqrt{{(x_i-x_0)}^2+{(y_i-y_0)}^2}-R$$

(5)

where $x_0$ and $y_0$ are the center location and R is the circle radius. We define asymmetric distance thresholds to identify the inliers and outliers (Figure 7). The asymmetric thresholds work equally for an circle inscribing technique. A point is considered as an inlier if −1 cm $⩽d_i⩽+2$ cm. The number of inliers is marked with $N_i$. Further, we constrain that $\frac{N_i}{N_w}<0.2$, where $N_w$ denotes the amount of points with $d_i<-1$ cm, as laser pulses are not able to penetrate into stems; thus, no points should be spread inside the cylinder. However, this assumption holds wrong under certain circumstances, such as due to co-registration error, texture of the tree bark, long ranges or wind effect. Therefore, we allow 20% tolerance. This restriction is important for stems that are only partially covered by TLS, as the example given in Figure 7. We would like to stress that the circle inscribing is only used to estimate the starting values for solving the nonlinear cylinder fitting problem. After M iterations, the consensus set $S_i$, which contains the largest $N_i$, is selected, and the circle associated with $S_i$ is considered as the solution. One can think that a circle is exactly the projection of a vertical circular cylinder. Therefore, the determined circle can be easily transformed to a cylinder by introducing a height along the z-axis. Consequently, a robust initial guess for the following cylinder fitting is generated.

Cylinder fitting features the same strategy as the circle fitting, except that it is performed in 3D. Thereby, we simply apply the RANSAC again on top of the standard least square solution. The work in [51] proved that seven is the minimum amount of points to achieve a unique solution of cylinder fitting. Therefore, a consensus set of seven points is randomly selected and fitted with a cylinder.

3.4.2. Cylinder Growing

The fitted starting cylinder is the basis for tracking the whole stem. The position of the next cylinder can be found by adjusting the starting cylinder upwards or downwards. We here take upwards growing as an example (Figure 8). The cylinder is shifted by a length of L. Here, L is 20 cm to give enough flexibility to accommodate the curvature of the stem. The axis of the shifted cylinder is then rotated to the vertical direction in order to simplify the subsequent processing by:

$$v_t=R\times v_a,$$

(6)

where $v_t={[0,0,1]}^T$ is the unit vertical vector. $v_a$ denotes the unit vector representing the axis of shifted cylinder. R is the rotation matrix. All points are rotated by R, as well. Consequently, a new coordinate system is introduced. The new axis is tilted vertically by a gradually increased angle from 0–α, with 20% of α for every step. α is determined proportionally to the height of the cylinder. The largest α is assigned to the lowest cylinder with 20°, because the stem tends to grow in a more straight direction on the upper part. The vertically tilted axis is then rotated on the horizontal plane clockwise (or anticlockwise) by a successive step angle of β (30°) until a full turn. The coordinate of the upper vertex of the axis is then:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill x_u& =x+Lsin\alpha cos\beta \hfill \\ \hfill y_u& =y+Lsin\alpha sin\beta \hfill \\ \hfill z_u& =z+Lcos\alpha \hfill \end{array}\end{array}$$

(7)

where x, y and z are the coordinates of the upper vertex of the directly shifted cylinder, which is determined from the previous cylinder. Furthermore, the radius r of the shifted cylinder is also adjusted by downscaling the value by 20%, gradually by steps of 4%, by assuming that the stem will become thinner at the upper parts, and vice versa. The shifted cylinder is tilted, rotated and scaled by testing all combinations of values of α, β and r. The goodness of fitting is determined by identifying the amount of qualified points with the same asymmetric criteria of the circle inscribing step. Consequently, the position of the next cylinder can be resolved by evaluating the goodness of each fitting. Finally, the determined cylinder is transformed back to the original coordinate system.

This elongation is continued upwards until there are not enough points to form a cylinder. For example, in the canopy, there are no sufficient points that are conforming the tested cylinder. The downwards elongation follows the same strategy until reaching the ground.

3.5. Post Processing

The cylinder fitting and growing allow the reconstruction of each tree stem with cylinders in every 20-cm section. We further apply a moving window approach to smooth the connection parts of successive cylinders for each stem. As shown in Figure 9, a moving window with a 20-cm height is placed on the connection part of two cylinders, which contains part of the data from Section 1 and part from Section 2. All contained data are tested against Cylinder 1 and Cylinder 2, respectively. The radius of one of the two cylinders that better fits all contained data will be selected. Consequently, each stem section at a specific height can be described by the radius and the center location. A truncated cone is then generated to represent each section, instead of a circular cylinder.

3.6. Evaluation

The automatically-reconstructed stems are compared with the reference measurements. The accuracy of DBH, location and diameter at various heights and the stem volume is evaluated. The stem volume is determined by:

$$V={\displaystyle \sum _{i=1}^m}{v}_i={\displaystyle \sum _{i=1}^m}\frac{\pi h}{3}(r_1^2+r_2^2+r_1{r}_2)$$

(8)

where h denotes the section height, $r_1$ and $r_2$ are the top and bottom radii and m is the number of sections.

The estimation accuracy will be evaluated by the bias and root mean square error (RMSE):

$$\mathrm{Bias}=\frac{1}{k}{\displaystyle \sum _{i=1}^k}(y_i-\widehat{y_i}),$$

(9)

$$\mathrm{RMSE}=\sqrt{\frac{1}{k}{\displaystyle \sum _{i=1}^k}{(y_i-\widehat{y_i})}^2},$$

(10)

$$\mathrm{RMSE}(\%)=100\times \frac{y}{\overline{y}},$$

(11)

where k is the number of observation data, $\widehat{y_i}$ denotes the reference value and $\overline{y}$ is the mean value of the variable.

4. Results

The results of the automatic stem location estimation are given in

Table 2. Results of the stem detection rate and reconstruction rate.

	Stem Detection	Stem Reconstruction
True positive	27	25
False positive	0	0
False negative	0	2
Completeness	100%	92.6%

. “Completeness” refers to the percentage of trees that are detected and reconstructed out of 27 trees, respectively. The results of stem reconstruction are summarized in

Table 3. Bias and RMSE of automatically estimated DBH, diameter, location and volume. * Results are evaluated from all heights along all stems.

	Bias	RMSE	RMSE (%)
DBH (cm)	0.03	1.80	5.50
DBH (TLS coverage >80%, cm)	0.05	1.50	4.58
Diameter (cm) *	0.13	2.45	8.94
Location (cm) *	1.60	2.09	-
Volume (cm3)	84.88	451.38	7.07

. All trees were successfully recognized. We constrain a bias of DBH worse than ±5 cm as reconstruction failure. Thereby, 25 out of 27 stems are successfully reconstructed. The estimated DBH shows a mean bias of 0.03 cm and an RMSE of 1.80 cm (5.50%) compared to the reference values. Figure 10 shows a high correlation between field-measured DBH and automatically modeled DBH values from TLS. The $R^2$ value is 0.99.

The stem curve is evaluated from two variables, the stem diameters and the locations of stem centers. TLS often cannot cover the stem on the upper part of a tree, especially inside the canopy, because of the occlusions in the forest. Therefore, the accuracy of the stem curve is compared based on the fact that both manual measurement and automatic modeling are performed as high up the stem as possible. What is worth mentioning, the uppermost model from TLS may even reach higher heights than manual measurements, because our self-adaptive method is able to identify a group of points that in fact belong to a stem, but are not necessarily clustered in a cylinder or at least a portion in space. Nonetheless, the average relative heights of stem top between TLS modeling and field measurement is 98% in this study, showing a good agreement. The mean bias of stem diameters from all trees is 0.13 cm and the RMSE is 2.45 cm (8.94%). The accuracy of the center of the stem has a mean bias of 1.60 cm and an RMSE of 2.09 cm. The box plots of these two variables are given in Figure 11 and Figure 12, respectively.

Stem volume is calculated from the stem curve using Equation (8). The result of the volume comparison also shows a good accuracy, with a relative RMSE of 7.07%.

5. Discussion

In this study, tree stems are automatically detected and modeled by a robust cylinder fitting scheme. Figure 13 shows an overall representation of our results. In the following subsections, the parameters used in the algorithms and possible error sources that could affect the results of the presented algorithms are discussed. At the same time, we address several challenges about the application of TLS in a landslide region, which potentially also apply to other high mountain heterogeneous forests.

5.1. Parameters

In the stem detection and modeling steps, various parameter were determined. Parameters used for stem detection include the search radius of the DBSCAN algorithm and the stem group extent enlargement factor. Stem modeling-related parameters include cell size, circle inscribing parameters, cylinder growing step length L and cylinder shifting and rotating parameters α and β. However, none of these parameters is sensitive to the methodology. Furthermore, they can be determined according to the datasets and field conditions.

The search radius of the DBSCAN algorithm depends on the point density [52]. The criterion is that it should be able to connect adjacent stem points. The stem group extent enlargement factor is used to retain points from the upper part of the stem in case the stem is inclined. We chose a factor of 1.5 because the trees in our study region are not heavily inclined. It can be increased for heavily curved trees or those with more irregular stem forms, such as deciduous trees.

Smaller value of cell size should be chosen if trees are grown close to each other. Besides, we used the −1 cm and +2 cm thresholds for inscribing the circle in 2D. They are chosen based on the fact that the mean bark thickness is ∼18 mm in our study region according to empirical functions (see Section 5.2). The step length L should be determined according to the stem shapes. For forests where trees are in general vertical (e.g., managed forests), L can be longer to simplify the modeling; while a smaller value should be used for complex stem shapes in order to represent the curvature of stems. The cylinder shifting and rotating parameters are also related to stem shapes. They can be explicitly refined if smaller shifting and rotating intervals for cylinder tests are required. However, this will increase the computation time accordingly.

5.2. Error Analysis

Terrain model estimation is a prerequisite to determine the normalized height of a tree. The accuracy of DTM, theoretically, could introduce a bias of DBH and stem curve, because of the false estimation of stem heights. However, we do not anticipate that the stem diameter would change sharply within a short distance. In addition, the field measurements also contain uncertainties arising from the human-determined absolute height of the tree [53]. Practical tests show that a change of terrain height by ±5 cm will only introduce a small change of the results of DBH determination by $0.2\%$–$0.6\%$; thus, this does not affect the results significantly.

Further, one of the main challenges in high mountain forest is the terrain accessibility. Therefore, the scan locations are often disjointed. Point coverage and density for every individual tree vary according to the scanning visibility, as well as the occlusions. To account for the effects of such coverage inconsistency, we plot the results of DBH against the TLS coverage percentage for every tree, as demonstrated in Figure 14. In general, the higher TLS coverage results in smaller absolute bias. The RMSE of DBH estimation deceases from 1.80 cm down to 1.50 cm for trees with TLS coverage larger than 80%, comparing with the results from all trees (Table 3). Nonetheless, the improvements are not significant, because our algorithm models the partially covered stems with a robust approach. In fact, the largest RMSE is smaller than 4 cm, demonstrating the capability of our approach for suppressing the deficiency of disjointed TLS scans.

We have shown in Figure 7 that the asymmetric thresholds are critical when fitting the circle or cylinder to points that cover only an arc. The fitted primitive could be heavily biased by the unbalanced weights introduced by the point distribution, if symmetric thresholds are given. To evaluate the effect of the asymmetric thresholds, we fit a circle to a stem whose cross-section is fully covered (i.e., full TLS coverage), by applying asymmetric and symmetric thresholds, respectively. We obtained an underestimation of ∼15 mm of diameter for the asymmetric thresholds. However, such a disparity is on the same scale with natural ambiguities, like tree bark roughness or stem texture. The asymmetric thresholds enforce the fitted primitive inside the stem points (i.e., inscribing), which is a plausible assumption in practice as laser pulses cannot penetrate into stems. However, the bark or surface texture does introduce a certain thickness of laser points around the stem (i.e., effectively forming a toroid). In other words, the inside primitive determines the diameter of the stem without bark. On the other hand, field-measured DBH was based on the peripheral surface of the bark, thus leading to an underestimated result compared to our algorithm. We calculated the bark thicknesses in our study site according to Equation (4) in [54]. The mean bark thickness is ∼18 mm, which is similar to the underestimation of 15 mm in our study. Therefore, the determined diameter is simply compensated by +15 mm in this study.

Asymmetric thresholds could result in significant underestimation of the stem diameter under one circumstance. For example, Figure 15 shows the cross-section at the height of DBH of the tree that is marked as reconstruction failure in our study site. The irregular shape significantly differs from a circle. Thereby, the fitted circle (or cylinder) from our algorithms will have a much smaller diameter. In the landslide region, such a phenomenon often exists, because of the disturbance from soil movements, which change the wood formation mechanism [55].

In addition, the stem shape becomes more anomalous on the parts close to the root (e.g., similar to Figure 15). Lower segments feature more complex stem shapes, especially when the root systems are partially exposed. Figure 11 and Figure 12 show that the biases of estimated diameter and location at the height of 0.65 m are usually larger than, for example, 1 m–2 m. Advanced curve fitting strategy should be developed for future studies, such as [27,56], to resolve such issues. Figure 11 and Figure 12 also indicate that the results’ accuracy drops along with the height starting from approximately 6 m above the ground. This addresses another challenge of using TLS in the forest. In particular, the multi-layered canopy structure in mountain regions blocks the laser pulse from reaching the top of the tree. A previous study also shows that there is a significant disparity of tree height estimation using TLS in urban heterogeneous forests (e.g., [31]). This issue can be possibly improved by combing data from other sources, such as unmanned aerial vehicle (UAV) laser scanning (e.g., [57]). Nevertheless, our self-adaptive cylinder growing strategy aims to model the stem up to the tree top if point data are available. The highest fitted cylinder reaches 12.3 m, which is consistent with the manual delineation.

5.3. Algorithm

Our algorithm is designed to meet the aforementioned challenges of processing TLS data from high mountain heterogeneous forests. The tree stem recognition method combines the z-normal value and projection density. The stem curve is retrieved by a robust cylinder fitting and self-adaptive growing scheme. The performance of our method is similar to those on a flat environment (e.g., [33,34,35]), demonstrating the effectiveness of precise stem modeling in mountain and landslide regions.

The unique capability of our method is that it is noise free. Many previous publications require fine estimation of stem points and then estimate or model the DBH (e.g., [33,34]). Yet, our RANSAC-based stem reconstruction method does not require a fine delineation of stem points and is robust with points from branches, leaves and other outliers. Figure 16 gives an example of one stem section with branches and outliers. Our method achieves a robust and accurate stem fit even when a significant number of disturbing points persist. This enables us to estimate tree stem locations instead of direct stem point delineation, which can be difficult due to the complex field conditions.

The presented cylinder growing strategy handles the challenge of occlusion. The growing cylinder always takes advantage of previous robustly fitted cylinder. In other words, the rough location of the next cylinder is confident. Therefore, the shifted cylinder connects points that belong to the stem even if they are partially occluded. Moreover, the cylinder length growth step can be adapted to cover the zone on the stem where the tree is occluded. For example, the step length can be elongated to 1 m or even longer, if necessary, although we fix it as 20 cm in this study.

The employed method in this study utilizes a 2D to 3D robust circle and cylinder fitting scheme. Most of the previous studies on using TLS focused on retrieving tree parameters, such as location, DBH and tree height (e.g., [26,28,33,58]). Some of them fit a circle to a slice of stem points (e.g., [33,34]). Errors may arise from the projection of points to 2D, because tree stems are not always strictly vertical in space, and the project points will form an ellipse if the stem is tilted. Thereby, the circle fitting becomes implausible. In our method, the purpose of circle inscribing in 2D is to provide a good initial guess for 3D cylinder fitting. The final model is performed in 3D by cylinder fitting; thus, it is not subject to a specific case if the stem is strictly vertical or tilted. For previous cylinder fitting approaches (e.g., [35,36]), errors could originate from branches and outliers if the standard least square is directly applied. Moreover, the solution of the nonlinear problem will be heavily affected by the initial guess, whereas in our method, it is determined by an RANSAC scheme.

Table 4. Summary of several selected studies on automatic stem detection and modeling from TLS data. Only the first authors are shown to save space. Some of the contents and values are visually inspected or calculated from corresponding publications. S and M in scan mode stand for single scan and multiscan, respectively. M(p) means that multiscan is only performed partially.

Study	Environment	Sample Size	Density (Stem/ha)	Scan Mode	Method	Completeness	Level of Automation	DBH Result (RMSE (cm) or R2)
Thies, M., et al. [ 59]	steep	50	556	S, M	Circle fitting	22% (S)	full	3.48 (S)
						52% (M)		3.22 (M)
Hopkinson, C., et al.[ 35]	flat	138	661	M	Circle fitting	93%	semi	$R^2$ = 0.85
					1.25 m–1.75 m
Watt, P. et al., [ 34]	flat	12	600	M(p)	Circle fitting	100%	semi	$R^2$ = 0.92
	planned
Wezyk, P., et al. [ 36]	flat	-	-	M	Cylinder fitting	63%–90%	semi	$R^2$ > 0.946
					1.28 m–1.32 m
Maas, H., et al. [ 33]	flat	14–29	212–410	S, M	Circle fitting	87%–100%	full	1.48–3.25
Yao, T., et al. [ 60]	flat	-	1017–3281	S	Angular width	-	full	7.0–8.0
Calders, K., et al. [ 61]	flat	65	317–347	M	Circle fitting	-	semi	2.39
Olofsson, K., et al. [ 26]	-	-	358–1042	S	Circle fitting	87%	full	2.0–9.6
					RANSAC	on average		14%
Moskal, L.M., et al. [ 31]	heterogeneous	25	-	S	Cylinder fitting	-	full	9.2
					voxel modeling			$R^2$ = 0.91
Liang, X., et al. [ 37]	flat	28	-	M	Cylinder fitting	-	full	0.82
					robust			4.2%
Brolly, G., et al. [ 32]	flat	213	852	S	Circle fitting	81%	full	4.2–7.0
					cylinder fitting
Our work	steep landslide	27	554	M(p)	Cylinder fitting	93%	full	1.8 (5.5%)
					RANSAC			$R^2$ = 0.99

gives a comparison of the performance of our work with previous studies. It is worth mentioning that Table 4 is not a comprehensive overview of TLS applications in forests, but it rather shows how different methods were applied to various forests and their corresponding performances. From the comparison, we can see that our method performed equally well or even better when compared to previous works that focused on flat regions or planned forests.

In addition, our method is also applicable for other tree species, such as deciduous trees, and the point cloud acquired from other approaches, such as UAV laser scanning. Figure 17 shows an example of our method applied to point clouds of a deciduous tree acquired from TLS and UAV laser scanning, respectively.

5.4. Applicability of TLS in Landslide-Affected Forest

In this paper, “landslide-affected forests” not only refers to the regions where intense landslides occurred. In fact, a wide range of high mountain forests share similar features. Soil movement is a common threat in many mountainous areas, such as the Alpine mountain range. Trees in these regions are interacting with soil as the tree roots try to cling to the soil. On the other hand, the moving soil drags the tree. Therefore, tree growth is disturbed, and stems typically have different shapes compared with flat areas.

Our self-adaptive cylinder growing strategy enables the modeling of the whole stem. Stem parameters are then calculated from the fitted models. This opens up the opportunity for precisely deriving more valuable forest variables, such as taper, stem volume and other information related to tree morphology. In particular, stem shape is of special value in high mountain landslide regions. Biomass estimation on the single tree level can be integrated into physical models to improve the understanding of tree-soil interaction (e.g., [11]). Different biomass components are also of great interest, for example, in estimating the stem, branches and leaf/needle biomass (e.g., [47]) and to quantify the leaf area index. Stem volume provides the basis for single tree biomass determination (e.g., [29]).

6. Conclusions

In this study, we present a full automatic tree stem detection and modeling method. The approach is tested for 27 trees in a landslide-affected forest in the Austrian Alps. A robust cylinder fitting scheme is exploited to reconstruct the whole stem. Stem curves are compared at various heights for the diameter and center location. Results showed good accuracies when compared to manually-measured reference data. Furthermore, we have discussed some challenges of TLS applications in landslide-affected forests. In general, the lack of application of TLS in high mountain forests calls for a specific point cloud processing approach, and our study highlighted the potential of the methodology. In the form demonstrated here, a limitation of our method is the handling of irregular stem cross-sections, which significantly deviate from a circle. The possible solution is to fit free-form curves. Our method can also be exploited for branch modeling. However, a higher point density is needed, especially for the upper parts of the tree. The combination of TLS with other data sources, such as airborne laser scanning, UAV laser scanning, optical images and hyperspectral images, is assumed to advance the research of forest management and other ecosystems research activities. In future studies, we will test our algorithm for other tree species and study areas.