Estimation of breast height diameter and trunk curvature with linear and single-photon LiDARs

Abstract

Key message

New technologies can take us towards real precision forestry: the terrestrial single-photon avalanche diode (SPAD) light detection and ranging (LiDAR) has a great potential to outperform conventional linear mode LiDARs in measuring tree parameters at the stand level.

Context

Precision forestry together with new sensor technologies implies Digital Forest Inventories for estimation of volume and quality of trees in a stand.

Aims

This study compared commercial LiDAR, new prototype SPAD LiDAR, and manual methods for measuring tree quality attributes, i.e., diameter at breast height (DBH) and trunk curvature in the forest stand.

Methods

We measured 7 Scots pine trees (Pinus sylvestris) with commercial LiDAR (Zeb Horizon by GeoSLAM), prototype SPAD LiDAR, and manual devices. We compared manual measurements to the DBH and curvature values estimated based on LiDAR data. We also scanned a densely branched Picea abies to compare penetrability of the LiDARs and detectability of the obstructed trunk.

Results

The DBH values deviated 1–3 cm correlating to the specified accuracies of the employed devices, showing close to acceptable results. The curvature values deviated 1–6 cm implying distorted range measurements from the top part of the trunks and inaccurate manual measurement method, leaving space for improvement. The most important finding was that the SPAD LiDAR outperformed conventional LiDAR in detecting tree stem of the densely branched spruce.

Conclusion

These results represent preliminary but clear evidence that LiDAR technologies are already close to acceptable level in DBH measurements, but not yet satisfactory for curvature measurements. In addition, terrestrial SPAD LiDAR has a great potential to outperform conventional LiDARs in forest measurements of densely branched trees.

1 Introduction

Forestry has been behind most other industries in the adoption of digital technology, but this is about to change. Recent studies are showing productivity increases in general agriculture at rates of 5 to 25 percent annually, with returns on investment of 1 to 2 years for digital technology, and similar gains are also being realized by some pioneers today for forest products. Inspired by advances in agriculture, forestry operators globally have begun pioneering the use of advanced technologies to improve forest management results. Among the most promising technologies and practices in precision forestry is the “Digital Inventory”: measurement of forest standing inventory—volume, species, and sometimes grade mix—by aerial remote sensing and in-forest devices (Choudhry and O’Kelly 2018). Forest inventory produces information for both strategic planning of forest management and for everyday forest operations such as tending and harvesting. With accurate forest inventory data, the value of harvested timber can be maximized, and the efficiency of operations can be increased. Potential wood buyers are also interested in knowing what kind of timber assortments is available in each forest stand to meet their specific supply needs.

Finnish Forest Centre is responsible for producing and updating stand-based forest inventory data of private forests, while the Natural Resources Institute Finland (Luke) conducts the National Forest Inventory (NFI) of all forests regularly in 5–10 years cycles (Luke 2020a; Tomppo et al. 2014). The Finnish Forest Centre employs remote sensing methods such as airborne laser scanning (ALS) with 5 pts/m², aerial photographs, and manual measurements from sample plots (Heikkilä 2017). ALS, based on LiDAR (light detection and ranging) technologies, provides precise 3D information on the structure of the forest, such as tree height, volume, biomass, and terrain as laser pulses reach the ground by partly penetrating canopy (Hyyppä and Inkinen 1999). The aerial photographs are useful in determining tree species. However, information produced by remote sensing and related supportive model computation (like model prediction and growth simulation) still requires field verification (Vähä-Konkka et al. 2020). Manually measured sample plots are necessary for increasing the accuracy of inventory and serve as a reference for modeling at the stand level (Metsäkeskus 2018). Currently operational forest planning relies on stand based forest inventory data, which should be as accurate as possible so that the forest owner can manage forests well and perform all operations according to best practices (Holopainen et al. 2013).

As presented by Mandlburger et al. (2019), LiDAR technology has three main approaches: (1) the conventional linear mode (or multi-photon or full waveform) LiDARs, (2) Geiger-mode LiDARs, and (3) single-photon LiDARs. The linear mode LiDARs utilize linear mode avalanche photodiodes (APD), whereas more modern airborne single-photon and Geiger-mode LiDARs employ single-photon avalanche diode (SPAD) detectors. The current terrestrial LiDARs are still based on linear mode APD limiting their operation in harsh conditions as the pulse detection threshold is in the order of hundreds of photons (Ullrich and Pfennigbauer 2016, Mandlburger et al. 2019). Because of higher single-photon sensitivity, higher dynamic range, and higher temporal resolution, the SPAD LiDAR can detect partly obstructed targets by means of time-correlated single-photon counting (TCSPC) in which the actual distance to target is based on a histogram defining the temporal distribution of the reflected photons (Henriksson et al. 2016, Fu et al. 2020). However, the low 3D image acquisition rate hinders application of SPAD LiDARs as single-pixel devices typically require data collection time of several minutes (Henriksson et al. 2016) and the latest SPAD array sensors tens of seconds (Morimoto et al. (2020). Nevertheless, higher 3D image acquisition rates, 10 Hz and above, have been demonstrated feasible with SPAD array sensors (Matsubara et al. 2018, Ruokamo 2019).

However, despite extensive research and huge development steps from traditional manual, often ocular forest inventory methods, the full economic benefits of digital forest inventory remain unobtainable. The current terrestrial LiDARs, if applied in harvesters, are too inaccurate compared to what the Finnish wood measurement law requires. In addition, measuring trunks of the trees in dense forests gets more unreliable as range measurement relies on detecting the most intense returning light pulses. Higher-density LiDAR data, with a UAV-LiDAR, for example, may increase the density of pulses/m² which can get between the branches, which would to some extent mitigate the problem of trunk obscuration, by producing enough data from the trunk surface. The accuracy of present portable LiDARs is limited to the level of several centimeters (Xie et al. 2020; Ghimire et al. 2017), whereas the Finnish wood measurement law requires accuracy of ± 4% for trunk diameter measured in the harvester head corresponding to ± 1 cm accuracy with nominal trunk diameter of 25 cm (Heikurainen et al. 2018). It should be noted, however, that Saarinen et al. (2020) have presented promising results of terrestrial laser scanning (TLS) with a ground-based geosystem LiDAR, placed on a tripod, which clearly limits the mobility and flexibility of measurements. Compared to movable or portable devices, accuracy of measuring, e.g., DBH, was significantly better with an average accuracy of 0.21 cm.

The majority of the previous research has focused on off-line processing of the point clouds, whereas improving the efficiency of the Cut-To-Length method would require measuring tree dimensions in real-time before cutting trees in the forest as proposed by Miettinen et al. (2010). Currently, the harvester operator selects the trees and cutting lengths based on a visual estimate and a given cutting matrix, which is error prone, because obstructing branches, weather and lighting conditions limit visibility from harvester cabin to the trunks of the trees. The sensors of the harvester head measure the actual diameter and length of the cut trees but ignore curvature of the trunk, although sawmills have varying directives for the maximum allowed log curvature (Suuriniemi and Marjomaa 1998).

In this paper, we introduce results in acquiring and estimating DBH and curvature for stand-based digital forestry inventories. We concentrate on key parameters reflecting the volume and quality of the trees in the stand: stem diameters and curvature of the trunks. We apply the promising new LiDAR technology, the SPAD technology, enabling more reliable range measurements for limited views, and also demanding environmental conditions such as rain, fog, canopy, and dense branching. In our proof-of-concept study, the primary aim was to compare the accuracy of the prototype SPAD LiDAR (VTT Technical Research Centre of Finland Ltd, Oulu, Finland) to conventional LiDAR technology, i.e., the commercial 3D LiDAR (ZEB Horizon by GeoSLAM) and manual measurement methods. Second, we aimed at demonstrating the penetrability of SPAD LiDAR by measuring obstructed trunks through dense branching. Keeping in mind recent developments in tripod-based scanning solutions, we rather focused on solutions fitting in harvesters. To our knowledge, this is the first study applying terrestrial SPAD LiDAR for measuring trees at the stand level, and we show preliminary evidence that SPAD LiDAR has a great potential for harvesters, outperforming traditional LiDAR technologies in conditions where the economically most valuable objects, the tree trunks, are partially obscured by branches, leaves, and foliage.

2 Materials and methods

2.1 Experimental setup

We measured a group of sample trees, including 7 Scots pine trees (Pinus sylvestris), with commercial hand-held 3D LiDAR ZEB Horizon (GeoSLAM, Nottingham, UK) and with a prototype SPAD LiDAR available at the VTT Technical Research Centre of Finland in Oulu. Stand characteristics are presented in Table 1. ZEB Horizon (Fig. 1) applies Simultaneous Localization and Mapping (SLAM) algorithm for localizing the sensor with respect to the environment and reconstructing the 3D point cloud in a common coordinate frame. While scanning the sensor measures 300 k points per seconds with maximum range of 100 m and with relative accuracy of 1–3 cm (GeoSLAM 2020). Both LiDAR methods studied in this work produce high point densities, up to thousands of points/m², scanned from a direction as close as possible to orthogonal to the trunk. Thus, compared to typical point density of ALS, 5 pts/m², much higher measurement points per trunk can be obtained.

Table 1 Forest inventory attributes of the sample stand per hectare. The actual sample plot forest stand was limited, with 0.18 ha (Metsäkeskus (Finnish Forest Centre) 2021)

Fig. 1

ZEB Horizon hand-held scanner head and data acquisition unit

Reconstruction of the point cloud with ZEB Horizon included two process phases. Firstly, the sample area was scanned, and the raw data was stored on the device’s local disk. Secondly, the raw data was uploaded to the GeoSLAM Hub cloud service that reconstructed the final point cloud of the sample forest (Fig. 2). Because of utilizing SLAM, the raw point cloud was in a common coordinate system. The raw point cloud measured with ZEB Horizon included ca. 93 million points so the points belonging to individual sample trees were cropped and saved to separate Point Cloud Data (PCD) files.

Fig. 2

The point cloud measured with ZEB Horizon including the whole sample area including 93 million 3D points

Figure 3 shows the prototype SPAD LiDAR comprising laser emitter, single-photon detector, and actuated mirrors for deflecting the coaxial laser and receiver beams. The maximum measurement range is ca. 20 m, and relative accuracy for distance is approximately 1 mm. The horizontal and vertical resolution of scanning is adjustable, but the mechanics limit the angle range of the deflecting mirrors between ± 20 degrees. The prototype SPAD LiDAR lacked the SLAM algorithm so only one side of the sample tree was measured from stationary location.

Fig. 3

The SPAD LiDAR prototype (VTT Technical Research Centre of Finland Ltd)

The measurements with VTT’s SPAD LiDAR were carried out sample tree by sample tree, and Fig. 4 shows the SPAD LiDAR targeted on the sample tree at a distance near 8–10 m. Because of mechanical limitations of the deflecting mirrors, the view was limited so that the root end of the trunk was scanned up to the height of 4–5 m. Duration of scanning depended on the selected angle resolution and with the used settings scanning of one tree took 5–10 min. The scanning raw data were stored on the disk on the host PC, and the off-line MATLAB application reconstructed the final point cloud. Figure 5a shows the cropped points of the sample tree measured with ZEB Horizon, and Fig. 5b shows the final point cloud including one sample tree acquired with the SPAD LiDAR.

Fig. 4

The VTT SPAD LiDAR equipment on a carriage and targeted at the sample tree

Fig. 5

The 3D points of the sample tree measured with ZEB Horizon (a) and VTT SPAD LiDAR (b)

To emulate measurement of densely branched trunks in forest, we brought a sample spruce (Picea abies) into the research facility allowing us to conduct experiments without a risk of exposing non-protected SPAD LiDAR equipment to harsh weather. To acquire comparable results, both ZEB Horizon and SPAD LiDAR scanned the sample spruce from a single direction according to the aspect shown in Fig. 6. The aim of limiting scanning to a single direction was to exclude effect of the SLAM algorithm, which could include measurements from sparse locations of the stem if ZEB Horizon moved around the tree. It should be noted here that, in practice, the actual scanning directions of a LiDAR integrated to a forest harvester would also be limited and dependent on the location of trees with respect to the logging tracks.

Fig. 6

The sample spruce (P. abies) with dense branching

2.2 Point cloud processing

Figure 7 shows the point cloud processing pipelines for both data measured with the ZEB Horizon and VTT SPAD LiDAR. The raw point clouds measured with ZEB Horizon and VTT SPAD LiDAR were pre-processed with the open-source software CloudCompare (https://cloudcompare.org/) in which the points belonging to each sample tree were cropped based on manual visual analysis and saved to separate PCD files. Each PCD file was imported to VTT’s proprietary 3D vision software for determination of DBH and centerline points of the sample trees by filtering with a median filter, segmenting and fitting cylinder models to the cylindrical point segments measured from the trunk. Finally, the curvature was analyzed in MATLAB based on the centerline points.

Fig. 7

The point cloud processing pipeline for determining breast height diameters and curvature of the sample trees

2.3 Determining diameter at breast height

The trunk DBH was determined based on point cluster cropped at the height of 1.2–1.4 m above ground level or neck of the root swelling. Then a cylinder model was fitted iteratively to the 3D points so that the distance between points and cylinder surface was minimized. The equations for iterative cylinder fitting are given in Appendix 1. The fitting algorithm has been implemented in C +  + language in VTT’s proprietary 3D vision software.

Figure 8a and b show the fitted cylinder rendered on top of the point cloud measured with ZEB Horizon, and Fig. 8c and d show the fitting result for SPAD LiDAR data. The DBH values were determined for each sample tree from data measured with ZEB Horizon and SPAD LiDAR.

Fig. 8

The cylinder model fitted to the 3D points measured with ZEB Horizon (a,b) and with SPAD LiDAR (c,d)

Figure 9a shows the residual errors of cylinder fitting to the 3D point cluster measured from the trunk at breast height with ZEB Horizon and Fig. 9b with SPAD LiDAR, respectively. The denser ZEB Horizon data contained 6820 points, whereas the sparser point cluster measured with SPAD LiDAR included 151 points. The mean absolute error (MAE) of cylinder fitting was 0.36 cm for the SPAD LiDAR data and 1.44 cm for the ZEB Horizon data.

Fig. 9

The residual errors of fitting cylinder to 3D point cluster for estimation of diameter at breast height a with ZEB Horizon data and b with SPAD LiDAR data

The mean absolute error was calculated as

$$MAE=\frac{\sum_{i=1}^{n}\left|{e}_{i}\right|}{n}$$

(1)

in which e_i is the residual error, i.e., distance to fitted cylinder surface, of i:th point and n is the sample group size.

As a reference for 3D measurements, we measured DBH values with a manual caliper (Masser 500). One DBH was measured per tree as an average of minimum and maximum values of two orthogonal measured diameters.

2.4 Determining curvature of a trunk

The curvature of a trunk of sample trees was measured using the LiDAR data from the trunk surface and then compared to measured caliper data. The curvatures were determined from root neck level up to the height of 4.2 m. From 3D data, the curvature was determined with respect to the centerline of the trunk by fitting cylinder models to the 3D point clusters cropped in 50 cm intervals resulting in 9 centerline points per trunk. The curvature was the maximum distance of the intermediate points from the line connecting centers of the butt and top end of the trunk. As the 3D line is defined by point p_line and direction vector v_line, the closest point on 3D line with respect to point p can be calculated:

$${p}_{closest}={p}_{line}-\frac{({p}_{line}-p){v}_{line}^{T}}{{v}_{line}^{T}{v}_{line}}{v}_{line}$$

(2)

Then the distance d between point p and 3D line is

$$d=\sqrt{{(p-{p}_{closest})}^{T}\left(p-{p}_{closest}\right)}$$

(3)

The reference for trunk curvature was manually determined by attaching a 4.2 m wire on side of the trunk and by measuring the maximum distance between the wire and the surface of the trunk. The lower end of the wire was attached to the root neck of the trunk and the top of the wire up to the trunk surface at 4.2 m distance. In practice, the curvature of a 4.2-m-long butt log was measured. As with DBH, one measurement per tree was taken. The curvature estimation done with respect to the centerline or surface of the trunk will produce slightly differing curvature values, but with typical dimensions of the sample trees, the difference is small. This is explained in the following, showing how to derive the curvature for the trunk along the centerline and curvature for the trunk along the surface.

Figure 10 shows the geometrical analysis of the cone with a known curvature B₁ with respect to the centerline s, known end points c and d, and known diameters D_t and D_b, reflecting the conic shape of the trunk. Here we assume that the bending is regular and the centerline of the trunk forms a circular segment between points c and d. The cf is half of the chord of the segment, and the curvature B₁ is perpendicular to the chord, defining also the point s, c, s, and d determine the point o. Now a and b can be calculated from points a and d, and the diameters D_t and D_b. Then the curvature B2 with respect to the surface of the cone is the difference of the radius os, the average diameter, and the cevian oe of the triangle oab:

$${B}_{2}=\left|os\right|-\frac{{D}_{b}+{D}_{t}}{2}-\left|oe\right|$$

(4)

in which |os| is

$$\left|os\right|=\frac{{\left|cf\right|}^{2}+{B}_{1}^{2}}{2{B}_{1}}$$

(5)

and |oe| is

Fig. 10

The geometrical analysis of curvature based on the centerline and the surface of the bent cone

$$\left|oe\right|=\sqrt{\frac{{\left|oa\right|}^{2}\left|be\right|+{\left|ob\right|}^{2}\left|ae\right|}{\left|be\right|+\left|ae\right|}-\left|ae\right|\left|be\right|}$$

(6)

Based on the geometrical analysis, the curvature value given by the methods deviate of order 0.014 cm (0.28%) if the trunk is an ideal bent cone with nominal length of |cf|= 2.1 m, centerline bent B₁ = 5 cm, base diameter D_b = 35 cm, and top diameter D_t = 20 cm. The difference between the methods is negligible as compared to the other error sources such as range measurement noise of the LiDARs, local irregular shapes of the trunk, and accuracy of the measuring tape used in manual measurement.

3 Results

3.1 DBH of the sample trees

DBH values of the sample trees are shown in Table 2. The DBH values varied from 25 to 40 cm, and values determined with different devices deviated typically 1–3 cm. Compared to manual caliper measurements, RMSE was 0.94 cm, and MAE (or bias) was − 0.47 cm for the ZEB Horizon, and RMSE was 1.77 cm, and MAE (or bias) was – 0.24 cm for the SPAD LiDAR.

Table 2 Measured breast height diameters

The measured curvature values of the sample trees are shown in Table 3 for the trunk length of 4.2 m, determined based on 3dD LiDAR data and manual measurements. The curvature values varied from 4 to 10 cm. The curvature values given by different methods deviated from 1 to 6 cm. Compared to manual tape measurements, RMSE was 2.59 cm, and bias was − 2.51 cm for the ZEB Horizon, and RMSE was 2.90 cm, and bias was – 1.69 cm for the SPAD LiDAR.

Table 3 Measured curvature from root neck to 4.2 m height

4 Detection of the trunk of densely branched spruce

Fig. 11a shows the raw point cloud of the densely branched spruce measured with ZEB Horizon. We segmented the point cloud with region growing algorithm available in Point Cloud Library; it was seen that even with various parameters, the result was unsuccessful in separating the points belonging to the trunk. Fig. 11b shows a characteristic segmentation result showing only clusters belonging to the branches that obstructed the trunk. The scanning direction was from right to left so the points belonging to trunk; if they exist, it should be located on the left in the images. Because of the dense branches, there were hardly any signals reflected from the stem, and filtering did not help.

Fig. 11

a The raw point cloud measured with ZEB Horizon from densely branched spruce (left view) and b the point cloud segmented to clusters. The black arrow indicates the direction of scanning

Fig. 12a shows the raw point cloud measured with VTT SPAD LiDAR. The points belonging to the trunk are distinguishable although relatively sparse points from the branches enclose the trunk. Some parts of branches were cut out from the measured point cloud, because the time gating excluded the photons reflected from the tips of the branches. In other words, the time gating influenced like a pass-through filter that determined the low and high limit distances from which the reflected photons may trigger the SPAD.

Fig. 12

a The raw point cloud measured with SPAD LiDAR from densely branched spruce and b the point cloud segmented to clusters

Points of the trunk surfaces were determined by a segmentation algorithm based on surface normals. Fig. 12b shows the segmentation result where points belonging to the trunk are in two cylindrical segments (labeled as 1 and 2) on the left. The two irregularly shaped segments on the right (labeled as 4 and 5) include points from branches, and these could be separated from the trunk segments, e.g., based on residual error of the cylinder fitting. The cluster number 3 located at the bottom of the image represents the base tub in which the spruce was standing.

5 Discussion

5.1 Quality of the derived tree parameters

The DBH values determined from the LiDAR data were at the expected level of precision when considering the range measurement accuracy of the devices and the reference DBH value measured with a caliper. The sample trees were upright oriented trees (P. sylvestris), and the laser was approximately perpendicular to the trunk and thus resulted in reliable range measurements. The MAE was clearly larger for the ZEB Horizon than for the SPAD LiDAR data. Still, even the ZEB Horizon results were close to acceptable level, i.e., ± 1 cm accuracy. This has, however, important consequences if the sensor were used in a harvester, and the view angles are limited, in extreme cases only one single view. Then the SPAD LiDAR would outperform the conventional LiDAR ZEB Horizon.

The deviation between curvature values determined from LiDAR data was higher than the specified accuracy of 1–3 cm of ZEB Horizon. This was an unexpected result as the DBH values determined from LiDAR data were in line with the validated accuracy of the devices. The large difference of curvature values may result from the GeoSLAM Hub SLAM algorithm that averaged range measurements acquired from different directions and locations around the sample trees. In contrast, the SPAD LiDAR performed only stationary measurements from a single direction, and the optical angle between the laser beam and trunk surface increased as the ray deflected from the horizontal plane towards top of the tree. The higher inclination angle may have distorted the SPAD LiDAR range measurements as simulated by Ullrich and Pfennigbauer (2018) and observed by Mandlburger et al. (2019) and Brown et al. (2020). Noticeably, the emitter and the detector of the SPAD LiDAR were located in the carriage (Fig. 4) close to breast height, which is a favorable location for measuring DBH values as the laser beam was near perpendicular to the trunk surface.

The scanning of densely branched spruce revealed the strongest benefits of the SPAD LiDAR over conventional LiDAR,. The points belonging to the trunk were distinguishable even from raw data of the SPAD LiDAR and region-growing algorithm of Point Cloud Library based on the surface normal, successfully segmented the points of the trunk. However, the numerous attempts to segment the trunk points from raw data of the ZEB Horizon were unsuccessful, indicating the point cloud lacked reliable range measurements off the trunk. An obvious reason for lacking points was that obstructed trunk reflected too few photons and was insufficient for triggering the pulse detection threshold of the ZEB Horizon.

6 Compatibility with standards

In Finland, the purpose the Act of Measuring Timber (FINLEX 414/2013) is to ensure reliability of the methods, equipment, and results of measuring unprocessed timber. Methods and accuracy of measuring are stated in the corresponding regulation (1323/14/2013). The Ministry of Agriculture and Forestry is responsible for supervising this law. In practice, this task is assigned to official measurers of the Institute of Natural Resources Institute of Finland (Luke 2020b). However, for measuring DBH, a standard is set only for measuring with a harvester head while cutting the logs. This is explained in the introduction of this paper. No official standards have been set on the accuracy of forest inventory, but the Finnish Forest Centre has set a target for the standard error ranging from 10% in tree height to 30% in stand volume (Heikkilä 2017).

7 Conclusions

In this paper, we introduced results in acquiring and estimating detailed forestry data for stand-based digital forestry inventories. We concentrated on key parameters reflecting the volume and quality of the trees in the forest stand: diameters at breast height (DBH) and curvature of the trunks. We applied new LiDAR technology employing SPAD detector and time-correlated single-photon counting, which demonstrated more reliable range measurements through dense branching than conventional LiDAR technology.

We compared the accuracy of the prototype SPAD LiDAR to conventional LiDAR technology, i.e., commercial 3D LiDAR ZEB Horizon and also manual measurement methods for determining DBH and curvature of trunk. The deviation of the DBH values from both LiDAR technologies were at the expected level and close to being acceptable, typically 1–3 cm, and they match to the specified accuracy of the measurement devices.

However, the deviation of the curvature values was unexpectedly high at a level of 1–6 cm. This result implies the inclined vertical surface at the top part of the trunk might have distorted SPAD LiDAR range measurements. In addition, the manual curvature measurement method was sensitive to the selection of locations for determining the maximum curvature of the trunk, which also reflects challenges in verifying the sensor-based measurements.

To our knowledge, this was the first study applying single-photon terrestrial 3D LiDAR for measuring trees in a forest stand, in this case at edge of the stand, and estimating tree parameters. We also showed preliminary evidence that SPAD LiDAR has the potential to outperform conventional LiDAR technologies in heavy branching conditions. However, further research and technological innovations are needed to increase the measurement speed and to develop a portable device that can be used in real forest conditions. Measurement speed can be increased by using a linear or matrix SPAD detector arrays (Ruokamo 2019) instead of single pixel SPAD used in this study.

Appendix 1. The equations for fitting cylinder model to 3D points

In cylinder fitting, the estimated parameters were position, orientation, and radius of the cylinder. Rotation matrix R and translation vector T define the position and orientation of the cylinder with respect to the point cloud coordinate system. The cylinder coordinate system’s origin was located on the cylinder axis and the Z-axis was collinear with cylinder axis.

Point p_ci in cylinder coordinate system

$${p}_{ci}={R}^{T}\left({p}_{L}-T\right)$$

(7)

in which R is 3 × 3 rotation matrix, T is 3 × 1 translation vector from point cloud coordinate system to the cylinder coordinate system, and p_L is a point in the point cloud coordinate system, i.e., LiDAR coordinate system.

Position, orientation, and radius of the cylinder were corrected iteratively as follows:

$$T=R\Delta T+T$$

(8)

$$R=R\Delta R$$

(9)

$$r=r+\Delta r$$

(10)

The position and orientations are with respect to the cylinder coordinate system, so four pose parameters were estimated together with correction of the radius Δr.

ΔT is translational correction vector

$$\Delta T=\left[{\Delta t}_{x}{\Delta t}_{y}0\right]$$

(11)

and ΔR the rotational correction matrix composed of rotations around x- and y-axes.

$$\Delta R=Rot\left(x,\Delta {\varphi }_{x}\right)Rot\left(y,\Delta {\varphi }_{y}\right)$$

(12)

$$Rot\left(x,\Delta {\varphi }_{x}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{cos}(\Delta {\varphi }_{x})& -\mathrm{sin}(\Delta {\varphi }_{x})\\ 0& \mathrm{sin}(\Delta {\varphi }_{x})& \mathrm{cos}(\Delta {\varphi }_{x})\end{array}\right]$$

(13)

$$Rot\left(y,\Delta {\varphi }_{y}\right)=\left[\begin{array}{ccc}\mathrm{cos}(\Delta {\varphi }_{y})& 0& \mathrm{sin}(\Delta {\varphi }_{y})\\ 0& 1& 0\\ -\mathrm{sin}(\Delta {\varphi }_{y})& 0& \mathrm{cos}(\Delta {\varphi }_{y})\end{array}\right]$$

(14)

The correction vector Δ in cylinder coordinate system

$$\Delta =\left[\begin{array}{c}\begin{array}{c}\Delta {\varphi }_{x}\\ \Delta {\varphi }_{y}\\ \Delta {t}_{x}\end{array}\\ \begin{array}{c}\Delta {t}_{y}\\ \Delta r\end{array}\end{array}\right]$$

(15)

The pose correction vector Δ was calculated as

$$\Delta ={-({J}^{T}J)}^{-1}{J}^{T}E$$

(16)

in which J is the Jacobian matrix and E is error vector.

$$E=\left[\begin{array}{c}{e}_{1}\\ \vdots \\ {e}_{n}\end{array}\right]$$

(17)

Error e_i is the distance between point p_ci and surface of the cylinder

$${e}_{i}=\sqrt{{p}_{cix}^{2}+{p}_{ciy}^{2}}-r$$

(18)

in which r is the radius of the cylinder. Jacobian matrix J includes partial derivatives of e_i with respect to Δ

$$J=\left[\begin{array}{c}\frac{\partial {e}_{1}}{\partial \Delta }\\ \vdots \\ \frac{\partial {e}_{1}}{\partial \Delta }\end{array}\right]$$

(19)

$$\frac{\partial {e}_{i}}{\partial \Delta }=\frac{\partial {e}_{i}}{\partial {p}_{ci}}\frac{\partial {p}_{ci}}{\partial \Delta }+\frac{\partial {e}_{i}}{\partial r}\frac{\partial r}{\partial }$$

(20)

$$\frac{\partial {e}_{i}}{\partial {p}_{ci}}=\left[\begin{array}{cc}\frac{{p}_{cix}}{\sqrt{{p}_{cix}^{2}+{p}_{ciy}^{2}}}& \frac{{p}_{ciy}}{\sqrt{{p}_{cix}^{2}+{p}_{ciy}^{2}}}\end{array}\right]$$

(21)

$$\frac{\partial {p}_{ci}}{\partial \Delta }=\left[\begin{array}{cc}\begin{array}{cc}0& {-p}_{ciz}\end{array}& \begin{array}{cc}-1& \begin{array}{cc}0& 0\end{array}\end{array}\\ \begin{array}{cc}{p}_{ciz}& 0\end{array}& \begin{array}{cc}0& \begin{array}{cc}-1& 0\end{array}\end{array}\end{array}\right]$$

(22)

$$\frac{\partial {e}_{i}}{\partial r}=-1$$

(23)

$$\frac{\partial r}{\partial \Delta }=\left[\begin{array}{ccc}0& 0& \begin{array}{ccc}0& 0& 1\end{array}\end{array}\right]$$

(24)

The iterative cylinder fitting algorithm requires initial values for R, T, and r (radius of the cylinder). The computation is most sensitive to the initial orientation, which needs to be approximately within ± 20° from the actual orientation of the trunk; otherwise, the estimation might not converge or converges to a local minimum. With processed LiDAR data, the initial cylinder axis was set collinear with the Z-axis of the point cloud coordinate system as the data was from standing sample trees. A reasonable initial value for T was computed as the mean of the fitted point cluster and initial radius was computed based on the 2nd eigenvalue of the point cluster.