An algorithm for the automatic parametrization of wood volume equations from Terrestrial Laser Scanning point clouds: application in Pinus pinaster

ABSTRACT

The parametrization of wood volume equations has traditionally been carried out with destructive samplings, which are highly resource-intensive. These equations must be specifically set up for each species and set of conditions, meaning that, in many cases, they are unfeasible or non-existent. Here, we present a nondestructive and fully automated methodology for the parametrization of merchantable volume equations from terrestrial laser scanning (TLS) data, which aims at being applicable to any species and stand typology. It is based on the estimation of diameters along the stem and the height of each tree, including a robust system for the automatic identification and correction of anomalous values. The implementation considers several types of volume equations, the most suitable equation being selected and parameterized using the diameter and height estimations. The methodology was tested in a Pinus pinaster plot with 428 trees, steep slopes, low branches and dense understory. The results showed that 97% of trees were automatically detected, and RMSE of the height and diameter estimations was 1.52 m and 1.14 cm, respectively. A volume ratio equation was automatically selected as the best option for the test dataset. RMSE in automatic volume estimations was 0.0233 m³, and 0.0149 m³ using diameters reviewed by an operator.

KEYWORDS:

• diameters along the stem
• total height
• anomalous sections
• Pinus pinaster

Highlights

1. Fully automated methodology for the parametrization of total and merchantable volume equations (stem taper and/or volume ratio) from TLS.

2. Parametrization derived from automated estimations of tree height and diameter along the stem of all trees in a plot.

3. Robust to the presence of branches and understory, and scalable to different plot sizes.

4. Incorporates an identification and labelling system for anomalous diameter estimations as well as mechanisms for their optional individual supervision.

5. Automatic selection of the most suitable volume equation.

6. Allows for supervision of diameter estimations and selection of the volume equation used.

1. Introduction

Estimating merchantable volume as accurately as possible is essential in forest management in order to know the amount of wood available for the various uses that industry demands. There are basically two different methodologies for determining merchantable volume: stem taper equations and volume-ratio equations (González et al. 2001; Trincado, Klaus, and Sandoval 1997). Stem taper equations rely on the mathematical relationships between diameters along the stem and the height at which they are located (Newnham 1992), while volume ratio equations estimate the volume of the tree up to a certain point in the stem (defined by its diameter or height) as a percentage of its total volume (Cao, Burkhart, and Max 1980; Van Deusen, Sullivan, and Matvey 1981). In both cases, tree stems are divided into logs, which are assimilated as a cylinder, and the top section is modeled as a cone. The total stem volume is obtained by summing the log volumes and the volume of the top section. In the few published studies where the two methodologies are compared, no significant differences between them have been obtained (Parresol, Hotvedt, and Cao 1987; Trincado, Klaus, and Sandoval 1997). In order to develop these types of equations, large and comprehensive mass data is required to ensure that variability in the population, as well as in the geographic areas, is covered, so that the models can be applied in different stand typologies (Picard, Saint-André, and Henry 2012). Traditionally, to obtain such measurements, destructive techniques have been used, according to which a number of suitable trees (not forked nor excessively branched) that provide a representative distribution of diameter and height classes are selected and felled. They are then cut into logs of variable length and two perpendicular diameters are measured for each cross section (Menéndez-Miguélez et al. 2014) up to a top diameter of 7 cm. These techniques are resource-intensive (Corral-Rivas et al. 2007; Crecente-Campo, Alboreca, and Ulises 2009), since they require the felling of trees and consequently, stem-volume models for many tree species are non-existent and/or have only been developed for certain specific geographical areas, thus limiting their applicability.

The use of remote sensing techniques, more specifically, Terrestrial Laser Scanning (TLS), is gaining great importance in the estimation of variables related to stems and branches. This is mainly due to the high resolution and precision of the information it provides at the plot level, making possible the 3D modeling (Zong et al. 2021) and reconstruction of the trees (Liang et al. 2016) without the need for felling. As a result, since the early 2000s the possibilities offered by TLS in the forestry sector are increasingly being recognized for research and commercial ends (Fröhlich and Mettenleiter 2004), particularly in terms of compiling forest inventories (Hollaus, Mokroš, and Wang 2019). Although manual extraction processes are still sometimes used (Holopainen et al. 2011; Kankare et al. 2014) the general trend is to use automatic procedures (Liang et al. 2012; Raumonen et al. 2015). Liang et al. (2018) detailed a wide range of automatic or semi-automatic forest inventory methods focused on stem detection and diameter at breast height (dbh) calculations.

Regarding parameter extraction, most studies focus on the estimation of stem variables such as diameters (d_i) at different height positions, including dbh – the diameter at 1.3 m above the top of the stump – and total height (h) (Cabo et al. 2016; Tiago et al. 2017; Hackenberg et al. 2014; Hauglin et al. 2014; Henning and Radtke 2006; Pfeifer and Winterhalder 2004; Thies et al. 2004; You et al. 2016) but there are some authors who have also studied crown variables (Jung et al. 2011; Mohammed, Majid, and Izah 2018). In the case of dbh, TLS-based estimations have provided RMSE values ranging from 0.71 to 2.64 cm (Hollaus, Mokroš, and Wang 2019) which fulfil the requirements of many practical applications, e.g. national forest inventories (Liang et al. 2018). This high accuracy makes TLS a suitable technique to develop locally adjusted merchantable volume equations (Sun et al. 2016; Trincado and Burkhart 2006).

Although there are many studies where stems are modeled, this number is drastically reduced when it comes to extracting diameter at different heights for merchantable volume equations. There are, in addition, only a few studies evaluating the performance of TLS in contrast to field data in the measurement of stem taper (Henning and Radtke 2006; Liang et al. 2014; Maas et al. 2008; Mengesha, Hawkins, and Nieuwenhuis 2015). That said, some authors have put their efforts into the use of TLS for the estimation of variables of great importance for the wood industry such as branch architecture (Laurin et al. 2016), solid wood volume (Dassot et al. 2012) and taper equation construction (Gabriel 2017; Olofsson and Holmgren 2016; Sun et al. 2016).

The consistent increase in popularity of LiDAR devices, linked to their decreasing cost and rapid developments in their associated computer hardware and scanner technology in recent years, will make 3D point cloud data easily available for a wider range of users. Indeed, robust, time efficient and flexible software for data processing is already being demanded (Tiago et al. 2017), especially in the field of forestry. Within this context, the automated processing of TLS forestry data remains a challenge, particularly in regard to the efficient storage and analysis of voluminous merged scan data, the filtering of noise and unwanted data, and methods for dealing with the partial and complete occlusion that is common in forest plots (Calders et al. 2020; Li et al. 2020; Pitkänen, Raumonen, and Kangas 2019). Another challenge is to reduce the time invested in data collection and processing to make it similar to or less than the time required to parametrize conventional volume equations. This depends largely on the level of automation that can be attained by the software-aided analysis of the scanned plots (Thies et al. 2004).

The aim of this study is to design, implement, and test a fully automated methodology for the parametrization of stem volume equations (stem taper and/or volume ratio equations). This is accomplished using data from TLS point clouds to automatically detect individual trees and measure diameters along the stem (d_i) and total height (h) at the individual tree level.

2. Methodology

Figure 1 shows a workflow of the methodology detailing the steps involved the automatic parametrization of volume equations, including the input/output data needed and the procedures to be followed. In Section 2.1, the procedure for obtaining h and d_i in a fully automatic way from the TLS point cloud is explained in detail, including the identification and individualization of the stems, and the subsequent estimation of their geolocation. Section 2.2 details the procedure followed for the parametrization of the volume equations.

Figure 1. Workflow of the methodology proposed for the volume equations parametrization from TLS point clouds (Input/output data are shown in yellow and procedures in green)

2.1 Automatic estimation of dendrometric variables

2.1.1 Preprocessing of point cloud data

This part of the methodology encompasses all the steps that need to be followed to process the 3D point cloud prior to the automatic estimation of d_i and h: the aim is to automatically detect individual trees in the plot, then isolate and identify all points belonging to each tree and extract those points, which are along the stem. The workflow can be itemized as follows (see Figure 2).

1. Construction of a digital elevation model (DEM) to define the ground.

2. Height normalization of the point cloud to obtain Z coordinate, which refers to the ground (Z_norm).

3. Selection of candidate point groups along a horizontal strip (equal height from the horizontal ground), and where no scrub or branches are expected, on the basis of clustering of points using the DBSCAN algorithm (Birant and Kut 2007). The consistency of the groups in the horizontal strip is checked according to the methods described in Cabo et al. (2018a; 2018b).

4. Location of stem axis through the determination of the centroid and PC1 of detected stems: The methodology assumes that candidate tree stems are essentially linear features. As such, in the principal component (PC) analysis of the XYZ coordinates of each stem-candidate point group, the first PC (PC1) aligns with the direction of the maximum variance from the XYZ coordinates of each group. Then, PC1 is used to define the direction of each stem axis. The other two principal components (PC2 and PC3) are, by definition, perpendicular to each other and to PC1 (Oviedo-de et al. 2021) and can therefore be used as local/tree coordinate axes to calculate the distance of any point from the PC1 axis.

5. Calculation and storage of distances from each point to the closest axis: D_axis This has various applications, but the most obvious is that it allows the points of each tree to be filtered by their distance from the axis, making the separation of the main parts of each individual tree into crown/branches and stem relatively straightforward.

6. Tree individualization. All the points are assigned and linked to their closest axis (Tree Id)

Figure 2. Workflow for the identification of tree stems prior to automatic estimation of d_i and h. 1) DEM generation 2) Point cloud normalization (Z_norm) 3) Strip on the height normalized point cloud where the groups of candidate points likely to be tree stems are identified 4) Tree stem axes resulting from PCA analysis, where the PC1 is represented by a red line 5) Axis distance mapping of each tree (D_axis) 6) Tree individualization procedure, where all the points are labeled with a tree identifier (Tree Id)

2.1.2 Tree height estimation

The estimation of tree height is based on the determination of the highest point in the crown in the proximity of each stem axis. For this, the points belonging to each tree (Figure 3A) are first clustered. From the resulting groups, those that are few in number and far from the stem, which may correspond to noise, artifacts, or canopy from higher neighbor trees, are automatically removed. The remaining points are filtered by distance from the tree axis in order to only consider those closest to the axis for the tree height estimation. This is done based on a threshold distance that delimits a cylinder around the top part of the tree, and can be varied for different species, conditions, or environments. h is finally estimated as the elevation above the ground of the highest point inside the cylinder (Figure 3B), which reduces the probability that h corresponds to the top of another higher tree, and allows appropriate estimations, even for tilted trees (Figure 3C). A default distance threshold value can be set as half of the average distance between neighbor trees in the plot. However, this value can be altered, generally reduced, if the stems are essentially straight and/or there is a high diversity in terms of tree height within the plot.

Figure 3. Tree height estimation procedure, where h is estimated for each individualized tree (A) as the highest point within the area of a delimited cylinder around the tree axis after noise filtering (B) and which is also valid for tilted trees (C)

2.1.3 Diameter estimation along the stem

The tree stem is initially segmented by extracting the points that are within a certain distance of each axis. This distance is set as the maximum stem radius expected in the plot. Then, diameters at different heights along each stem (d_i) are measured at a user defined interval starting from the ground (Figure 4A). To do so, each d_i is measured in a slice, which comprises all the points in a thin horizontal section of stem of uniform thickness. The maximum height with points in these sections strongly depends on the quality of the results obtained in the stem filtering step, which, in turn, is mainly dependent on occlusions, the point cloud density and the presence of artifacts around the stem, such as branches or understory. These three factors are interrelated and their influence on the resulting d_i is expected to be stronger in the upper sections of the stem. As well as the height interval for slice extraction, the thickness of the slices and the upper height limit can be specifically set for different species and/or conditions, environments or studies (Figure 4B).

Figure 4. A) Cross-sections along the stem B) Detail of stem slices showing height interval for their extraction, and their thickness C) Circle fitting to the points of a stem slice

For the initial diameter estimation, the algorithm computes the diameters by fitting a circle to the points in each slice in each tree stem (Figure 4C). The circle fitting method is described in detail in Cabo et al. (2018a;b). The algorithm calculates the diameter and the X and Y coordinates of each section center, thus their position along the stem is defined and stored.

2.1.3.1 Automatic identification and labeling of anomalous sections

Frequently, the horizontal sections at different heights used to fit the circles contain points that do not belong to the stem but to other parts of the tree, such as branches, leaves or other artifacts. In those cases, the circle fitting may result in an erroneous estimation of the stem radius (r). The algorithm includes two complementary functions that allow, on the one hand, the identification and labeling of the anomalous sections, and, on the other, automatic correction in those cases where it is possible. Both functions are explained in detail below:

1. Inner circle (IC): this function was presented for the first time in Cabo et al. (2018a). Here, a concentric inner circle with radius r’ = r/2 is set up in the fitted section (Figure 5). If the position and size of the circle fitted is correct, no points are expected inside the inner circle.

2. Sector occupancy (SO): this function is based on the assumption that the points within a section follow a close-to-circular distribution and, except in cases of very severe occlusions, almost half of each tree section should be covered by points if it is measured from only one TLS setup, but more if it its measured from several setups. In order to verify that the circle fitting is compatible with this point coverage around the stem sections, an annulus divided into sectors is created around the position of each initially fitted circle (Figure 5). The width of the annulus can be varied if irregularities on the stem/bark are expected. If the position and size of the circle fit is correct, at least half minus one of the sectors should be occupied by points: this also allows circle fittings in stem sections that are only reached from one scan position, while the shape of the point group is continuously checked.

Figure 5. Checking functions of the circle fitting step in each section (IC +SO) showing an example of an ideal case, where the section would be labeled as correct

If in the first circle fitting of a section there are no points in the inner circle, and there are sufficient sectors occupied, the diameter estimation is considered correct. However, if there are points within the inner circle, or there are not sufficient sectors occupied in the annulus, it is highly probable that the diameter fitting is not correct. To resolve this, all the points within the section are clustered based on the distance between them. Then, the circle fitting is carried out again using only the cluster with the largest number of points, which is more likely to contain a higher percentage of points on the stem (Cabo et al. 2018a). This gives a new estimation for the diameter and position of the section of the stem. The new circle fitting is then checked again (i.e. IC and SO).

If a section passes the validation functions (either in the first or second round), it is labeled as “correct (C).” If both the first and second round of validation functions fail, the section is labeled as “flagged (F).” If in the first or the second round of validation functions there are not enough points to perform the circle fitting, the section is labeled as “not enough data (ND).”

Figure 6 shows 10 different hypothetical section scenarios and the behavior of the circle fitting and validation functions in each case. They are classified into four groups: True positives (TP), True negatives (TN), False positives (FP) and False negatives (FN). TP are correct estimations (from the first or second fitting round) and they are automatically stored as correct diameters (Figure 6. (A, B, C, D)). TN (i.e. those sections correctly labeled as F or ND) can be either left to be inspected visually by an operator or eliminated straightaway (Figure 6. (E, F, G, H)). FP and FN correspond to rare situations in which the algorithm would wrongly classify the results. For example, Figure 6I corresponds to an irregular stem section, while Figure 6J represents a branch with a quasi-circular section. However, these specific misassignments could be avoided by changing the parameters of the checking functions/tests (diameter of the inner circle, annulus width, and/or minimum number of occupied sectors).

Figure 6. Algorithm for diameter estimation performance when labeling sections. Sections can be classified into 4 categories, TP) True positives: labeled as C when they indeed are; TN) True negatives: labeled as F when they indeed are; FP) False positives: labeled as C when they should be F; and FN) False Negative: labeled as F when they should be C

2.2 Parametrization of volume equations

The merchantable volume can be calculated using different types of volume equations: mainly stem taper and volume-ratio equations, which are similar in terms of accuracy with respect to merchantable volume. In order to parametrize either of them, it is necessary to have a longitudinal data structure, that is, multiple diameter measurements for each individual at different heights along the stem (d_i, h_i) as well as the total height (h). The estimated dbh, d_i and h obtained from TLS data following the methodology described are used as input values to estimate merchantable volume using volume equations. For each type of volume equation, a large number of models (i.e. different specific expressions of the general volume equation) and parameter estimation approaches have been proposed in recent decades (Westfall and Scott 2010).

The first type of equation, a system of taper and volume equations, describe stem taper by providing (i) a taper equation with the diameter at any point along the stem or the height of the stem for a fixed diameter (eq.1), (ii) merchantable volume to any top diameter and from any height or individual volume for logs of any length at any height above the ground (eq.2) and (iii) total volume by integration of taper equation (eq.3). Taper equations can be classified into: simple taper equations, segmented taper equations and variable exponent taper equations. Of these, segmented taper equations have the greatest flexibility, as well as the capacity to describe the total volume efficiently (Pang et al. 2016). Ideally, a taper equation should be compatible, meaning that the volume computed by integration of the taper equation should be equal to that calculated by a total volume equation (Clutter 1980; Demaerschalk 1972). The specific model included in this work, that proposed by Fang, Borders, and Bailey (2000), meets this requirement as it includes a segmented taper equation (specific to f_1, eq.1), and complementary formulations to directly calculate total and merchantable volumes (specific to f_2, eq.2; specific to f_3, eq.3).

(eq. 1) $d_i=f_1\left(dbh,h,h_i\right)$(eq. 1)

(eq. 2) $v_i=f_2\left(dbh,h,d_i\right)$(eq. 2)

(eq. 3) $v=f_3\left(dbh,h\right)$(eq. 3)

It is important to consider that, in spite of the versatility and accuracy of stem taper equations, their practical use is limited unless they are embedded within a computer program (Larsen 2017). Moreover, when the data source is TLS, it is highly probable that the convergence in the iterations needed to calibrate the taper equation could be compromised due to the difficulty of d_i extraction in the upper stem because of branch occlusions in the crown, which leads to poor delimitation of the stem profile.

The second type of equations, volume ratio equations are, in contrast to taper equations, directly targeted at estimating wood volume up to a predefined diameter depending on its merchantable volume. They are very easy to both develop and use (Crecente-Campo, Alboreca, and Ulises 2009), and a hypothetical lack of upper stem data would not have a strong influence on their parametrization. Volume estimations from volume ratio equations, whose general expression is shown in eqs. 4–6, are based on the combined use of two equations – a total volume equation ($v$) (eq. 4) and a ratio equation, R_i (eq. 5) – to estimate the portion of volume up to a certain point (limit diameter d_i or height h_i) with values from 0 to 1, and, as such, the merchantable volume ($v_i$) is their product (eq. 6):

(eq. 4) $v=g_1\left(dbh,h\right)$(eq. 4)

(eq. 5) $R_i=g_2\left(dbh,d_i\right)$(eq. 5)

(eq. 6) $v_i=v.R_i$(eq. 6)

The algorithm developed in the current work, includes five total volume models (v; specific to g_1, eq.4) and four ratio models (R_i; specific to g_2, eq.5). The five specific total volume models (v, g1) are the linearized allometric model and four other models proposed by Spurr (1952): combined variable, generalized combined variable, quadratic, and cubic polynomial. The four ratio models (R_i; g_2, eq.5) included are the ones described in Burkhart (1977); Clark and Thomas (1984); Reed and Green (1984); Van Deusen, Sullivan, and Matvey (1981). Finally, the algorithm processes all the models and selects and integrates one of the total volume models (eq.4) and one of the ratio models (eq.5) into a volume ratio equation (eq.6), based on statistical criteria.

Considering all the above, the steps implemented in the algorithm are synthesized as follows:

1. In the first step, outliers in the stem analysis data (d_i, h_i) are detected following the systematic procedure proposed by Bi and Long (2001). A nonparametric taper curve is fitted using local regression with a typical smoothing factor value (0.25). Then, those values, which fall below a point equivalent to the lower quartile minus twice the interquartile range, or above a point corresponding to the upper quartile plus twice the interquartile range were removed owing to the fact that most of these data points correspond to stem deformations, and volume equations are not intended for deformed trees (Rodríguez, Torre, and Oviedo 2015).

2. The models implemented in the methodology (stem taper equation and the volume ratio models) are fitted and the algorithm calculates the values and significance of parameters (fitting by generalized least squares for non-linear models), the goodness-of-fit statistics of every model (coefficient of determination R², root mean square error RMSE, Akaike´s information criterion, AIC) are used to select the best model and measure the quality of the fittings, which can also be checked by graphical analysis (residual versus predicted values and predicted versus observed values and a scatterplot of relative diameter (d_i/dbh) vs relative height (h_i/h)).

3. After the debugging of the data set, total volume is calculated using a cone equation for the top section and, for each log, the Smalian equation (a formula that calculates the volume for each log by multiplying the average of the areas of the two end cross-sections by the section’s length). This total volume then becomes the dependent variable in the parametrization of the different equations.

4. The algorithm automatically selects the model with the lowest AICd (that is, the smallest AIC value in the set of models) where all parameters are significant as its candidate for the final merchantable volume equation. The user can either accept the algorithm choice or reject it and select any of the other proposed models according to their own preferences or specific requirements on the basis of all the goodness-of-fit statistics and graphical behavior of each model.

3. Case study

3.1 Study area and data

The methodology explained above was tested in a coniferous stand (Pinus pinaster) located in the autonomous region of Asturias (Figure 7A), where Pinus pinaster is the third species in terms of wood volume logging, representing 21% of total coniferous wood felled in Asturias (SADEI 2018).

Figure 7. A) Location of the study area B) Schema of the real distribution of the scans and the spheres in the test case plot C) Pinus pinaster stand and understory within the test case plot D) 3D representation of size and slope of the test case plot

The study plot, which has an approximate area of 5700 m² (Figure 7B), is part of a breeding program where the growth of different families and provenances of Pinus pinaster is being evaluated (RTA2017-00063-C04-02, 2017). Initially, in 2005, 900 trees were planted following a gridded distribution. A thinning treatment was made in 2018, reducing the number of trees to 428. Regarding the stand characteristics, density is 806 trees/ha with an average dbh and h of 18 cm and 10.5 m, respectively, while the average canopy cover fraction is 44%. The terrain is irregular with a slope of 60% and dense cover of different species of scrub such as Ulex, Ericas, Pteridium among others (average height is 50 cm, canopy cover is 80%) (Figure 7C).

Data acquisition was carried out in November 2018 with a TLS model FAROFocus3D (Faro 2018). In total, 24 scans were necessary to ensure the full coverage of the study area and minimize the effect of occlusions. With the aim of merging the point clouds of individual scans into a unified coordinate system, polystyrene spheres of 25 cm diameter fixed to surveying rods were used. Also, the position of the spheres was measured in the field using GNSS (Global Navigation Satellite System), which has an accuracy of 1 cm, thereby ensuring that the unified point cloud had absolute coordinates. The RMSE of the registration was 4 mm, and the final matched point cloud obtained for all scans was approximately 150 million points, after removing duplicates within 6 mm. The scanning distribution is shown in Figure 7B. Finally, for all trees within the plot, dbh was measured in the field with a caliper (dbh_f) to the nearest 1 mm and total height (h_f) to the nearest 10 cm, using a digital hypsometer (Vertex IV 360º).

3.2 Assessing algorithm performance

The performance of the algorithm was assessed from two perspectives: tree detection rates and deviation in the estimation of diameters and tree heights. Tree detection rates were evaluated in terms of completeness and correctness with respect to the tree locations measured with the GNSS+Total station. Completeness was computed as the proportion of the trees within the plot that were detected by the algorithm (number of correctly detected trees divided by the total number of trees in the dataset) and Correctness was calculated as the proportion of elements identified as trees by the algorithm that were actual trees (number of correctly detected trees divided by the total number of elements detected as trees in the dataset). The deviation in the estimation of both diameter measurements, dbh and d_i, and height, h, was evaluated by comparing the results obtained from the algorithm and the values obtained using (i) traditional inventory methods and (ii) manual/visual measurements of the point cloud data.

In order to obtain the manual/visual measurements, a graphic interface was designed and implemented.

1. For manual tree height estimations (h_m) each tree is always shown from two different perspectives and the user can identify the highest point from both perspectives before selecting which to use, ensuring misassignments are avoided (e.g. if ‘Perspective 1ʹ in Figure 8A was the only available perspective of that specific tree, the h_m estimation would be wrong). The elevation of the selected point above the ground is stored as the total height of each individual.

2. For manual diameter (d_m) estimations, the horizontal distribution of the points included in each section is graphically displayed on the screen. Over each section, the user can draw 3 points (A, B, C) along the section contour, so that a circumference is automatically fitted to them (Figure 8B), after which the diameter and the coordinates of the fitted circle are stored. In the case of dbh, the deviation in the diameter estimated by the algorithm is also calculated by comparing it with the direct measurements in the field made with a caliper (428 sections: one per tree). In the case of the diameters along the stem, the algorithm estimations are compared with the manual measurements from the point cloud (9844 sections: 428 trees X 23 sections).

Figure 8. A) Manual total height estimation in the point cloud by marking the highest point in the cloud from two different perspectives B) Manual estimation diameter at different heights in the point cloud, by drawing 3 points (A, B, C) around the section contour

With the aim of minimizing the influence of the operator on the estimations, all height measurement and diameter fitting processes were conducted twice, each time by a different operator (Op 1 and Op 2).

3.3 Experimental results for automatic estimation of dendrometric variables

3.3.1 Stem detection rate

All the features detected as trees were actual trees with no false positives present, therefore, the correctness of the algorithm is 100% for the dataset. In regards to completeness, only 14 trees of the total 428 remained undetected, giving a value of 97%.

3.3.2 Total height

The algorithm automatically measured the total height of 98% of the trees. The average values obtained by the algorithm (h), by traditional inventory methods (h_f) and by manual measurement in the point cloud (h_m) as well as their standard deviation, are shown in Table 1. RMSE obtained when comparing h with h_f was 1.52 m.

Table 1. Mean and standard deviation (in meters) of the total height obtained from: algorithm (h), traditional inventory methods (h_f) and manual measurements of the point cloud data (h_m)

	h	hf	hm
Mean	8.78	10.29	8.72
Standard deviation	1.44	1.48	1.43

Figure 9A shows the differences analyzed tree by tree. More precisely, it shows probability density functions for the differences between h, h_f, and h_m. The differences between h and the h_m are less than 0.5 m for 83% of the trees, and less than 1 m for 94%, meaning that the operator and the algorithm each estimate a very similar height in the vast majority of cases. Although differences between h_f and h_m are higher, they are less than 2 m in 80% of cases.

Figure 9. A) Probability density functions for height differences comparing h_m versus h (blue line) and h_f versus h_m (orange line)

3.3.3 Diameters along the stem

In this dataset, due to the complexity of the forest stand structure (young trees with a great number of branches, along with high stand density and steep terrain; Figure 10A), the automatic estimation of diameters was carried out from 0.5 to 4.9 m in height (i.e. from the average height of the dense understory to approximately half the average height of the trees in the plot; Figure 10B). The circle fittings were computed for sections with a thickness of 5 cm and spaced 20 cm apart vertically (i.e. 23 sections per tree; Figure 10C)

Figure 10. Visual results obtained by the algorithm for the diameter fitting in the study plot: (A) initial point cloud (B) clean stems (C) detail of circle fitting along the stem

Regarding dbh, the algorithm measured it fully automatically in 93% of the trees within the plot, while in the remaining 7%, diameter fitting was not possible or was labeled as F or ND by the validation functions. The dbh average value was 17.27 (±3.2 cm) and the RMSE was 1.14 cm when compared with dbh_f. In absolute terms, 56% of the trees have an error value lower than 1 cm, and this increases to 82% for an error value of less than 2 cm. The distribution of the size of the differences in diameter estimation is shown in Figure 11 (pink line).

Figure 11. Probability density functions for diameter at breast height differences comparing dbh_m versus dbh (blue line) and dbh_f versus dbh (pink line)

Figure 12 shows the comparison of the diameters along the stem estimations obtained by the algorithm and the manual measurements carried out by the two operators. This comparison has a double aim: (i) comparing operator 1 and operator 2 provides information about the disparity of the different reference measurements and allows the identification and elimination of outliers (in cases where the differences between operators are high) and (ii) the comparison between the algorithm and the operators provides an estimation of the accuracy of the method while allowing the subjectivity associated with the operators’ measurements to be considered.

Figure 12. Boxplots of the differences between A) manual estimation of operator 1 (Op 1) and operator 2 (Op 2); B) and C) comparison of the manual estimation of Op1 and Op2, respectively, with algorithm estimation excluding sections flagged as candidates for review. Q1 and Q2 represent the percentiles, and min and max refer to 1.5 times the interquartile range in mm from the upper and lower quartiles respectively. On the right of each boxplot the blue line represents the number of sections (as a percentage) that were used to build the boxes for each section

Specifically, Figure 12A shows the differences between the estimations of the two operators in those sections where both could see and draw a circle. These differences are approximately of the same order and remain constant along the whole stem: they show apparent symmetric distributions and no bias at any specific height.

Figure 12B and 12C show the differences between the diameter estimations along the stems as calculated by the operators and by the algorithm in the sections that were not flagged for review or elimination by the algorithm. The differences between the two operators are very small according to the average interquartile range (8.2 mm) and remain constant along the stem. However, when the estimations of each operator are compared to those from the algorithm, the interquartile range is slightly higher (Op1-alg: 11.2 mm and Op2-alg: 24.0 mm) and the difference increases with the height of the section. Specifically in the case of dbh, RMSE was 1.02 cm when compared with dbh_m and the distribution of the size of the differences in diameter estimation is shown in Figure 11 (blue line). On the right of each boxplot the blue line represents the proportion of sections that were used to build the boxplots for each section. The general tendency in each of the 4 boxplots is for the number of sections to decrease with height. This is to be expected because of the presence of low branches and the high density of trees, a combination, which makes diameter fittings more difficult higher up the stem, particularly when the branches of the crown begin to appear (at around 2.5 m in this test case).

3.3.3.1 Anomaly detection in stem sections

Figure 13 shows the distribution of the sections along the stem in terms of how the algorithm classifies them: correct sections (C: gray), flagged (F: magenta), and those where there was not enough data to compute the circle fitting (ND: blue). Also, to assess the performance of the algorithm when identifying ND, and to compare the results with the manual estimations, ND sections have been subdivided into NDa (not enough points for either circle fitting or visual/manual estimation: dark blue) and NDb (not enough points for circle fitting, but sufficient for visual/manual estimation by the operators: pale blue).

Figure 13. Labeling results obtained by the algorithm expressed in terms of percentage of sections classified into each category

As demonstrated in Figure 13, in general terms, the proportion of both F and ND sections is lowest around breast height (1.2–1.6 m) where stems are less affected by occlusions from understory or tree branches, and, in addition, the incidence angle of the laser beams is in general more favorable. In these sections, more than 90% of the diameter estimations are labeled as C, and less than 2% as ND.

The proportion of F sections increases steadily on either side of breast height (i.e. below 1.1. and above 1.5 m), and remains stable above approximately 2.5 m, the average height of the first branches. Conversely, the proportion of ND increases consistently with height from breast height (same increase rate for NDa and NDb up to the last two sections, where NDa becomes much larger than NDb)

3.4 Experimental results for parametrization of volume equations

Selected taper and volume ratio equations were fitted using two datasets that include the following variables for each tree: d_i, h_i and h. The data sets differ in the way in which these three variables have been calculated. In the first (hereafter, “Automatic”), all the variables were estimated automatically by the algorithm. The second (hereafter, “Supervised”), contains the same data as the Automatic, but those sections it labeled as F or ND were substituted by the diameter manually estimated by an operator, i.e. the data from Op1 and Op2 that was also used to evaluate the automatic diameter estimations. There were a small number of measurements, which were very different between the operators, and these were considered outliers, and removed from the dataset. However, in general, measurements from the two operators were very close (<5 mm) so the use of data from either of them would provide virtually the same result. On this basis, the diameter from either Op1 or Op2 was selected randomly for each section.

Figure 14 shows the scatter plot of the relative diameters (d_i /dbh) vs relative heights (h_i/h) that were used for the parameterization of the merchantable volume equations. In both cases there is a reduction in the number of sections that were available for the model parametrization in the upper part of the stems, especially in the top third of the stem. For the supervised data the situation is slightly better than in the Automatic, because there are more sections in the band of relative height (h_i/h), corresponding to 0.60 to 0.80 whereas in the Automatic data there are virtually no sections in that interval. Moreover, there are more sections providing results in the Supervised data and a smaller proportion of outliers than in the Automatic data (7384 vs 5362 sections and a proportion of outliers of 9.7% vs 11.6.%, respectively).

Figure 14. Scatterplots of relative diameter (d_i/dbh) vs relative height (h_i/h), with a local regression smoothing curve (red line) for the two datasets used in fitting merchantable volume equations Note: point (1,0) represents the final section of the stem, where d_i = 0 (d_i/d = 0) and h_i = h (h_i/h = 1)

For both datasets, the algorithm offered the following outputs: parameter estimation, goodness-of-fit statistics and graphical analysis (residual and predicted vs. observed data) of the merchantable volume equations. Looking first at the taper equation system, it did not converge in either of the two data sets due to the lack of sections in the upper third of the stem. From among the volume ratio equations, the algorithm selected for each data set was that with the lowest AICd and where all parameters were significant, which in this particular case was the allometric volume equation (eq. 7) and the ratio equation proposed by Clark and Thomas (1984) (eq. 8).

(eq. 7) $v=b_{0}db{h}^{b1}{h}^{b2}$(eq. 7)

(eq. 8) $R_i=exp(b_3{d}_i^{b4}db{h}^{b5})$(eq. 8)

where v is total over bark stem volume (m³), calculated as the sum of the volume of all the logs within the stem, and the top of the tree; R_i is the proportion of total volume accounted for at diameter (d_i); dbh is diameter at breast height (cm); d_i over bark diameter at height h_i (cm); h is total tree height in m, and b₀, b₁, b₂, b₃, b₄ and b₅ are the parameters to be estimated.

Table 2 shows the number of sections used and the comparison of the goodness-of-fit statistics for the volume ratio equation selected by the algorithm, for each of the datasets.

Table 2. Goodness-of-fit statistics of the model fitted in the different data sets

	Data set	n	P	RMSE	R^2
Total volume equation	Supervised	7382	3	0.0149	0.8480
	Automatic	5362	3	0.0233	0.6714
Ratio equation	Supervised	7384	3	0.1630	0.6759
	Automatic	5362	3	0.1747	0.6382

Note: n is the number of sections used for the equation fitting, P is the number of parameters in the model, RMSE root mean square error, and R² is the coefficient of determination.

The Supervised dataset gives better results than the Automatic, especially in the case of the total volume equation. The results demonstrate that a low number of sections, and more specifically their absence in the upper third of the tree, is what is responsible for the drop in the values of the model statistics. Finally, Table 3 shows the estimation of the model parameters, which were significant at the 5% level in the two datasets.

Table 3. Parameter estimated for the models fitted in the different data sets

Data set	Total volume equation			Ratio equation
	b0	b1	b2	b3	b4	b5
Supervised	0.00008	1.833	0.7875	−0.5664	4.044	−3.797
Automatic	0.000091	1.659	0.9396	−0.7997	3.946	−3.814

Note: b₀, b₁, b₂, b₃, b₄ and b₅ are the estimated parameters in the volume ratio equations.

4. Discussion

4.1 Comparison with previous work

Comparing the performance of the methodology employed in this work with previous studies is not straightforward as plot size, tree density, forest structure, etc. vary between trials (Cabo et al. 2018a) and all have a strong influence on the results. Moreover, the test case plot has the additional difficulties of high slope, young, and therefore small trees, a tendency of the trees to tortuosity and, particularly, the high density of branches and understory. Despite this, it is possible to establish a range of values obtained by other authors for each dendrometric variable and for volume equations that can serve as a reference for discussion. Regarding tree detection, it is strongly influenced by stem density (Liang and Hyyppä 2013; Thies et al. 2004): the lower the tree density in the sample plot, the higher the detection rate of the stems (Chen et al. 2009). According to Liang et al. (2016), in sparse forests (100–200 stems/ha), the detection rate may reach 100%, while in a forest with a stem density of over 1000 stems/ha, the stem detection rate is generally around 70%. The results obtained in this study for stem individualization are extremely good, with 97% success in a plot with 806 stems/ha. The approach of selecting only the points close to the projection of the stem, based on the estimation of a linear axis for each tree, indeed outperforms previous works by some of the authors (Cabo et al. 2018a; 2018b) as it reduces the probability of assigning a tree the height of a higher neighbor tree and eliminates points from branches, leaves, or other features around the stem for a more effective diameter estimation.

4.2 Tree height estimation

Regarding the results obtained for the calculation of h, most of the differences between the algorithm estimations and the field measurements are between 1 and 2 m. Tree height underestimation using TLS datasets is a common issue (Olofsson and Holmgren 2016), explained by the fact that the light beam of the sensor is not capable of reaching the top of the crowns, due to foliage obstruction, so the treetops are very often not well defined in the resulting point cloud. Nevertheless, when comparing h with h_m, the differences are barely noticeable (6 cm of difference between average values), which indicates that the limitation does not lie with the algorithm, but rather in the nature and quality of the data itself. In cases where h needs to be estimated with high accuracy TLS data could be complemented with aerial data, such as LiDAR (Hadas et al. 2017; LaRue et al. 2020; Puletti et al. 2020) or photogrammetric data (Aicardi et al. 2016; Tian et al. 2019) or with a subsample of field measurements in those cases where no remote sensing data is available.

4.3 Diameter estimation

Concerning the dbh versus dbh_f comparison, RMSE was 1.14 cm, which is compatible with the values of between 0.71 and 2.64 cm reported in previous studies (Hollaus, Mokroš, and Wang 2019). The difference between Automatic and manual procedures as diameters along the stem increases is also consistent with previous research (Henning and Radtke 2006; Maas et al. 2008) and can be attributed to low point density in higher parts of the stem, the smaller object surface area and obstructions caused by branches and foliage (Liang and Hyyppä 2013). This fact also explains why the most accurate estimations occur in the intermediate part of the tree (between 0.7 and 2.9 m, approximately), which is the interval considered to be “occlusion free.” The distribution of anomalous sections (F and ND) follows the same pattern, that is, they are more frequent in the lowest and highest parts of the stem. This problem with diameter estimation in the upper and lower sections could possibly be improved by scanning at higher resolutions, which could lead to better penetration of the foliage, and more laser points. However, any improvement would probably be limited because higher resolutions would still not reduce the occlusion effect (Liang et al. 2016). Stem profile measurement could also be improved by combining diameter measurements retrieved from different sources, or, as shown in the supervised method used in the test case, by incorporating visual checking or manual fitting in challenging sections if very high accuracy is needed. This might include the possibility of selecting the trees to be used in the parametrization of the equations.

4.4 Parametrization of volume equations

For the parametrization of the merchantable volume equations, the number of sections available, especially in the upper third of the stem, clearly influences the convergence of the stem taper equation and the accuracy of volume estimations. In this particular test case, the use of supervised data improves on the results obtained with the Automatic data, resulting in an increase of R² both in the stem volume equations (from 0.67 to 0.85) and in the ratio equations (from 0.64 to 0.68), as well as the reduction of RMSE (from 0.0233 to 0.0149 m³). When these results are compared with those obtained by destructive techniques in Pinus pinaster, where R² for stem volume equations ranges between 0.91 and 0.98 (Alegria and Margarida 2011; Yousefpour et al. 2012) and for volume-ratio equations is between 0.80 and 0.89 (Teshome 2005), the values for the same parameter obtained in this study might seem to be low. However, a comparison in absolute terms is not realistic because: (i) in destructive techniques only healthy and standard-shaped trees for which there is a full range of measurements along the stem are considered and (ii) in contrast to this study, where all the trees in the plot were considered, in destructive sampling, young trees with erratic growth are generally excluded. Consequently, it is only to be expected that results would be better in study plots with adult trees, low-density understory, and the reduced presence of branches.

4.5 Algorithm implementation performance

In terms of time, destructive sampling is by far the most time-consuming step involved in the parametrization of volume equations by means of traditional techniques (Yu et al. 2013). In the case of TLS sampling, only one day would be necessary for a sample of 100–150 trees, which is the number of trees commonly used to adjust merchantable volume equations (Crecente-Campo, Alboreca, and Ulises 2009). The authors have estimated that only 3 minutes per tree (16GB RAM, i7 4-core 8 G processor) is needed to obtain the parametrized volume-ratio equation (including scanning field work), which is a substantial improvement in time compared to destructive methods. Furthermore, in the case of the supervised dataset, only 20 seconds was necessary to review a complete tree (i.e. 2.5 hours of intense work for all the trees in the test plot). Consequently, the supervised option is cost-effective in terms of time invested vs improvement in the merchantable volume equations.

4.6 Transferability of the method

This study is eminently methodological, and it has been applied and evaluated in a specific plot of Pinus pinaster to assess the performance of the method. The results from a single study plot can be transferred to others as long as they have similar characteristics to the test case plot. However, the strong point of the methodology presented here is that it provides the possibility of parametrizing the volume equations of any tree species or stand typology. The only requirement, is to scan an adequate number of trees (in one or several plots, depending on the scope of the application of the volume equations) so that the algorithm has the necessary input diameters and tree heights (from TLS point clouds) to parametrize the volume equations, after which the algorithm selects the most suitable equation for each case.

From the point cloud processing perspective, the basic principles of tree detection and diameter estimation have already been tested in plots with different characteristics in benchmark studies (Liang et al. 2018; Hollaus, Mokroš, and Wang 2019). These studies describe how more than 20 different algorithms work in a wide variety of scenarios (i.e. different point cloud and vegetation densities, plot configurations, and point cloud technologies) and provide very good results when compared with traditional field measurements. However, it is possible that fully automatic approaches for diameter estimations give poorer results in stands with very complex structural or orographic conditions (Liang et al. 2018; Cabo et al. 2018a, 2018b; Hollaus, Mokroš, and Wang 2019). In this sense, the methodology here offers contingency alternatives for such complex scenarios: it automatically analyses the coherence of the tree measurements (diameters along the stem), flags anomalous estimations and includes the possibility of the visual checking and correction of flagged measurements. Furthermore, this can be also applied in not-so-complex scenarios, when, for instance, there are few trees available in the TLS datasets and thus problematic or anomalous sections cannot be discarded. What is more, while the methodology selects the most appropriate option (based on statistical indicators: AICd), it also offers the possibility of selecting a different model, when, for instance, the reliability of a specific model has been demonstrated to be better for a particular species and/or plot configurations (Westfall and Scott 2010; Menéndez-Miguélez et al. 2014; Saarinen et al. 2019).

4.7 Influence of measurement errors

Regarding the management of possible measurement errors, in forest inventory modeling it is common to assume that the independent variables are free of measurement error (Morales-Hidalgo, Kleinn, and Scott 2017). Also, some studies have compared the impact of measurement errors on the parametrization of taper functions when using direct traditional methods and indirect nondestructive methods like TLS (Rodríguez, Torre, and Oviedo 2015; Marchi et al. 2020; Torresan et al. 2021). Most of these studies concluded that the errors are similar or that they are smaller when using indirect (i.e. remote and nondestructive) measurements. Nevertheless, considering an error model for the independent variables could reinforce the analysis, so this may be something worth exploring in future studies. In the specific case of this study, using error models for the independent variables would, however, affect the automation of the method, as it would involve the need for a specific error model for each scanner and/or stand characteristic, which would potentially impact the simplicity and usability of the method.

5. Conclusions

In this study, a fully automatic methodology for the parametrization of stem merchantable volume equations from TLS data was developed and tested. The methodology, applicable to any species or stand typology, automatically detects the trees within a plot, and measures diameters along the stem (d_i, including dbh) as well as total height (h) at the individual tree level. These variables are then used by the algorithm as direct inputs for the parametrization of a set of merchantable volume equations (one stem taper equation and several volume-ratio equations) from which the one with the lowest AICd is automatically selected. The methodology incorporates an automatic system for identifying and labeling stem sections, which present anomalous diameter estimations so that they can either be discarded or manually reviewed.

The methodology was tested in a dense Pinus pinaster plot with 428 trees, steep slopes, low branches and dense understory. The automatic estimations of tree height and diameter were compared with field data and manual measurements made on the point cloud data. The results show good tree height estimation, but a clear misrepresentation of treetops. As for diameters along the stem, the best estimations are in the intermediate part of the stem and have no apparent bias. Discrepancies with reference data, however, increase with height, as does the number of sections that are flagged for review or elimination. Consequently, the number of sections available for the merchantable volume equation fitting is reduced, especially in the upper third of the stem. This fact has a strong influence on both the type of merchantable volume equation that can be fitted, and the degree of uncertainty in the merchantable volume estimation.

In the test case, the manual review of sections flagged for review or elimination improved the volume estimation. Nevertheless, the results obtained using the fully automatic approach are still suitable for the parametrization of the equations when there is no need for a high degree of accuracy in the merchantable volume estimation Also, the nature of the methodology suggests that better results could be obtained with less challenging plots, but in order to achieve a well-founded analysis of the performance of the methodology, further analysis covering a wider range of test cases is needed.

Finally, it is important to note that TLS sampling overcomes the traditional drawbacks of destructive methods, providing a 3D reconstruction of the forest at the moment of the scan. This has many advantages such as multitemporal studies with various objectives, enabling new variables of interest to be studied retrospectively, and the possible estimation of other stem shape variables (closely related to wood quality).

Data and Codes Availability Statement

The data that support the findings of this study are available with the identifier 10.6084/m9.figshare.14046578 at the private link: https://figshare.com/s/50286f1ab59f5da937bd.

Acknowledgments

The authors wish to thank CETEMAS field workers (Manuel Alonso-Graña López and Ernesto Menéndez Álvarez) for carrying out the scanning of the whole plot. Thanks also to Ronnie Lendrum, scientific editor and proof-reader, for correcting the English of the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Instituto Nacional de Investigación y Tecnología Agraria y Alimentaria (INIA) in the framework of the RTA2017-00063-C04-02 (2017) project entitled : “Evaluation of relevant characters for the sustainable management ofPinus pinaster Ait. and their interaction under new climatic scenarios„ and by the UK Environment Research Council (NE/T001194/1).