Measuring stem diameters with TLS in boreal forests by complementary fitting procedure

Point clouds generated by terrestrial laser scanners (TLS) have enabled new ways to measure stem diameters. A common method for diameter calculation is to fit cylindrical or circular shapes into the TLS point cloud, which can be based either on a single scan or a co-registered combination of several scans. However, as various defects in the point cloud may affect the final diameter results, we propose an automatized processing chain which takes advantage of complementing steps. Processing consists of two fitting phases and an additional taper curve calculation to define the final diameter measurements. First, stems are detected from co-registered data of several scans using surface normals and cylinder fitting. This provides a robust framework for localizing the stems and estimating diameters at various heights. Then, guided by the cylinders and their indicative diameters, another fitting round is performed by cutting the stems into thin horizontal slices and reassessing their diameters by circular shape. For each slice, the quality of the cylinder-modelled diameter is evaluated first with co-registered data and if it is found to be deficient, potentially due to modelling defects or co-registration errors, diameter is detected through single scans. Finally, slice diameters are applied to construct a spline-based taper curve model for each tree, which is used to calculate the final stem dimensions. This methodology was tested in southern Finland using a set of 505 trees. At the breast height level (1.3 m), the results indicate 5.2mm mean difference (3.2%), −0.4mm bias (-0.3%) and 7.3mm root mean squared error (4.4%) to reference measurements, and at the height of 6.0m, respective values are 6.5mm (3.6%), +1.6mm (0.9%) and 8.4mm (4.8%). These values are smaller compared to most of the corresponding contemporary studies, and outperform the initial cylinder models. This indicates that the applied processing chain is capable of producing relatively accurate diameter measurements, which can, at the cost of computational heaviness, remove various defects and improve the modelling results.


Introduction
Terrestrial laser scanning (TLS) is a relative recent technology in the production and analysis of forest inventory data. Through high-resolution three-dimensional data, TLS acquisitions offer a great potential to improve the quality of information which traditionally has been estimated using relatively limited field data and pre-existing models, or calculated through expensive (i.e. time consuming) measurements from felled trees. In particular, the capabilities of TLS systems to automatically, non-invasively and rapidly produce a large pool of accurate data related to trees' structural and biophysical metrics have been found highly beneficial (Liang et al., 2016;Watt and Donoghue, 2005). Taking advantage of these opportunities, TLS data has been successfully applied to various purposes at single tree, plot and stand levels, such as diameter, stem curve and height measurements, crown width estimation, volume and biomass calculations, and species detection (Åkerblom et al., 2017;Liang et al., 2014;Maas et al., 2008;Olschofsky et al., 2016;Srinivasan et al., 2015). Repetitive TLS measurements can also be applied for change detection purposes such as assessing increased volume or biomass (Kaasalainen et al., 2014;Mengesha et al., 2015).
One of the common TLS-based tasks in forestry has been to measure stem diameters of single trees, which further enables, for example, assessing their volumetric dimensions or timber assortment distribution (Dassot et al., 2012;Kankare et al., 2014). Ideally, good data quality is assured by selecting the specific trees to be scanned, designing scanning strategy based on these predefined targets, and performing fieldwork only in suitable conditions related to e.g. wind and phenology (Dassot et al., 2011;Liang et al., 2016). In reality, however, emphasis is normally given on plot-level study designs, which focus on scanning multiple trees by using a limited number of separate scans. To avoid excessive amount of manual work in the data analysis, automatic procedures capable of detecting single trees without significant omission or commission errors, extracting their stems, and calculating reliable stem diameters at various heights are normally in a key role during the data analysis.
Reliable automatized calculation of stem diameters requires carefully planned fieldwork and data processing procedures. One important initial determinant is the decision whether to use only one scan, based on a point cloud representing a single full field-of-view, or a set of multiple scans, which are normally made both inside and outside of the intended area of interest in order to collect complementary data. In the case of multiple scans, clouds are co-registered and combined with the help of artificial reference targets or natural objects in the scene, or matching is later performed at feature level (Liang et al., 2016;Newnham et al., 2015). Acquisition speed, simplicity of the design and exclusion of operations needed for co-registration are favoring single scan setup (Olofsson and Olsson, 2017). This makes single-scan applications to be also a good option to be used in operationally oriented purposes, which aims at maximizing sample size while simultaneously minimizing time and costs required for field work (Astrup et al., 2014;Kelbe et al., 2015).
The most obvious drawback of single-scan data, however, are viewing limitations deriving from occlusion effects, which may result in lower capabilities to detect individual trees and reconstruct their structures Kankare et al., 2017;Maas et al., 2008;Wilkes et al., 2017). As occlusion is further highly affected by stem density (Liang et al., 2012;Pueschel et al., 2013), single-scan data is likely to be most deficient in denser forests. Due to this reason, acquisition of multi-scan data is favored in many studies despite of its more tedious collection. Furthermore, determining the scanning locations in relation to the scanned trees as well as selecting appropriate scanner settings, such as point density, have also important implications to the final results. These topics however are not covered here in details, and readers are referred to Liang et al. (2016) and Wilkes et al. (2017) for related discussion.
Once TLS data has been collected and single clouds have potentially been merged into a combined data set, single trees and their structures can be extracted. There are various techniques applied for this purpose in different studies, such as voxelization, meshing, or detection of specific geometric shapes, as well as their combinations and further extensions (Dassot et al., 2012;Newnham et al., 2015;Olschofsky et al., 2016). Voxels do not initially expect any given geometric shapes, but rather aim at conceptualizing and representing the point-filled geographical space as a set of volumetric elements (Wu et al., 2013). This has been found to be useful for e.g. estimating tree volume, biomass, leaf area or canopy structure (Béland et al., 2014;Cifuentes et al., 2014;Hackenberg et al., 2015;Hosoi et al., 2013). Mesh-based methods, in turn, focus on constructing a tree structure as a form of mesh, which is further used to extract the needed measurements (van Leeuwen et al., 2013;Wang, Z. et al., 2014). When aiming at to produce accurate diameter calculations along the whole stem, however, fitting procedures based on geometric shapes are often preferred and measurements are typically made by dividing the trunk into longer vertical cylinders Liang et al., 2012;Raumonen et al., 2013), or detecting shorter cross-sections using circular shapes (Olofsson et al., 2014;Pueschel et al., 2013).
If diameter calculation using geometric shapes is automatically applied for the whole stem, the used algorithm must also be able to process trees which are not growing straight, have curved stem, or include occluded regions (Olofsson and Holmgren, 2016;Wang et al. 2016aWang et al. , 2016b. Regardless of carefully planned fieldwork, at least partial occlusion will always remain in the resulting point cloud, which deteriorates the shape fitting. Occlusion near to the scanner is most severe upwards from the level where canopy starts to induce shadowing effects, whereas upper parts of the tree have often better visibility from further away (Liang and Hyyppä, 2013;Watt and Donoghue, 2005). This supports the use of multiple scans, but the final extent of occlusion is affected by factors like branching structure and whether the scanning is made in leaf-on or leaf-off conditions (Olofsson et al., 2014;Pueschel et al., 2013;Srinivasan et al., 2015). Merging multiple scans, however, may also result in additional defects in terms of dislocated tree halves, which fail to join together as a complete shape. These co-registration problems may not only derive from defects in merging the distinct scans, but will also appear due to wind. In windy conditions, trees are moving during the scanning event as well as between the scans, which may considerably hinder the diameter measurements (Åkerblom et al., 2012;Henning and Radtke, 2006;Seidel et al., 2012;Vaaja et al., 2016). The only way to prevent wind-derived inaccuracies in co-registered data is to collect data only by still weather, which may be impractical or even impossible depending on the data needs or location.
In addition, circular or cylindrical shapes can only offer simplifications of the actual stem form, which may not always correspond to the reality. These relatively simple shapes are favored over more complex geometric primitives or irregular shapes mainly due to their computational simplicity and robustness (Åkerblom et al., 2015), but they may not work properly with deformed stems (Kankare et al., 2015;Saarinen et al., 2014). Selecting another symmetrical or even asymmetrical shape would be possible, but such an approach would considerably increase calculation efforts and simultaneously complicate outlier detection, i.e. exclusion of non-stem points such as TLS returns from branches or leaves. Sensitivity to outlier points depends on the selected shape detection method, but their removal is normally tied to a specific shape and regarded as a highly necessary action prior to fitting procedures (Paláncz et al., 2016;Pueschel et al., 2013). Given the reality when working with TLS data, outliers will always exist and this must be considered in the fitting strategy (Olofsson and Holmgren, 2016).
As point clouds are potentially affected by various defects, approaches for detecting stems and calculating their diameters should on one hand be robust enough to handle with deficient data, but on the other hand be capable of improving the resulting accuracy when allowed by the data quality. To respond to this challenge, we developed and tested an automatic processing chain, which consisted of two separate phases. At the first phase, stems were recognized and modelled from the co-registered multi-scan point cloud using a relatively robust surface normal and cylinder-based iterative modelling approach, which was run several times. This provided the initial stem structures and their dimensions by using all the potential data, but also induced certain defects deriving from the scanning conditions, forest structure, and the modelling methodology itself. First, the constructed cylinders must be long enough to provide an acceptable fitting solution, but this will simultaneously decrease their capabilities to follow the stem details. Second, except of knowing that the cylinder represents the best possible fit with certain preconditions, the quality of the solution is largely unknown. Third, given the occlusion effects from ground vegetation, branches and other trees, the modelled stem generally remains too short. And fourth, co-registration errors will hinder the cylinder construction, which are often exaggerated towards the treetop. Furthermore, cylinder modelling is an iterative process which may result in different structural solutions, and therefore performing several modelling attempts is preferred to assure an acceptable outcome.
Due to the potential errors described above, the second modelling phase aimed at confirming the cylinder-derived results and, if necessary, fine tuning the diameter measurements. This included selection of the best cylinder models, extracting thin slices along the modelled stems, and processing these slices using a circle fitting procedure. The goodness of the fit was assessed based on the continuity of the circle edge, referring to a sufficient coverage of TLS points along the expected edge zone, which was first evaluated using the co-registered point cloud. If an accepted solution was not found, points of the corresponding slice were split into single scans, which were processed separately to find the best estimate for the slice diameter. Then, finally, a taper curve model based on the refined slice diameters was constructed to balance out potential over-and underestimates, and used to calculate stem diameters at heights of 1.3 (i.e. breast height) and 6.0 m. These results, along with cylinder-based diameters, were finally compared to field measurements from corresponding heights.

Study area and tree measurements
Field data for this study was collected in summer 2016 from Lapinjärvi, southern Finland (60.7°N, 26.1°E), which belongs to southern boreal forest zone. The study area, however, is a designated research forest established in 1933, including a wide variety of management regimes and therefore differing forest types within a relatively small area. Tree measurements and TLS scans were done from 18 plots ( Fig. 1), which were selected to represent variable forest types in terms of species composition, age and stem density. From each sample plot, all trees having at least 45 mm diameter at breast height (d1.3) within a radius of 9.00 m from the plot center were selected to be measured. For all the selected trees, d1.3 was measured using Masser Sonar electronic caliper and tree species was recorded. D1.3 measurement was performed at a perpendicular direction against the line towards the plot center. If the trunk was significantly deformed at the breast height, measurement point was lowered to the smallest non-deformed diameter below the breast height level. Furthermore, two additional measurements were performed for tally trees, which were included in an angle count plot (Bitterlich plot) with q = 1.5 (Tomppo et al., 2011). These measurements included diameter at 6.0 m height (d6, only for trees taller than 8.0 m) and tree height. D6 was measured using a manual caliper, and height using Vertex hypsometer. The angle count based selecting strategy ensured that larger trees had a higher probability to be measured compared to small ones, while keeping the amount of fieldwork at tolerable limits.

Terrestrial laser scanning
Each of the 18 field-measured plots was scanned using Leica P40 terrestrial laser scanner. Leica P40 is a time-of-flight scanner, which is able to scan up to 1 million points per second, uses < 0.23 mrad beam divergence and ≤3.5 mm beam diameter at front window, and offers a 3D position accuracy of 3 mm at 50 m distance. As tree measurements and laser scanning were performed by different field groups, their timing could not be fully synchronized due to logistical reasons. Time lag between these two varied from one day to a maximum of 35 days, being approximately two weeks on average. Given the relatively short time lag and normal variation in the measurement accuracy, however, this was not regarded to induce major effects on the gained results. Scanning was always performed in dry conditions, having temperatures between 17 and 24 degrees centigrade and average wind speed between 1.6 and 6.5 m/s as recorded during the scanning at the closest windmeasuring weather station, located at 40 km distance from the fieldwork area (Harabacka, Porvoo).
At each plot, one scanning station was placed at the plot center and four stations around the plot edge (9 m distance) at all the cardinal directions (N, E, S and W). Having one station in the middle and additional stations evenly placed over the sample plot is a commonly used setting in TLS studies, offering normally a good visibility on the scanned forest stand without an excessive use of resources (Abegg et al., 2017). To further optimize the scanning results, ensure the proper inclusion of T.P. Pitkänen et al. ISPRS Journal of Photogrammetry and Remote Sensing 147 (2019) 294-306 all the measured trees and minimize occlusion effects, however, the actual station points were allowed to deviate up to 1.5 m from the predesignated locations based on plot conditions. Scans were acquired by using a spacing of 3.1 mm of adjacent points at 10 m distance. Separate scans were co-registered with black-and-white target boards using Leica Cyclone 9.1 software, resulting in plot-wise absolute errors between 1 and 3 mm and having single target errors up to 11 mm. Coregistered point clouds were then restricted to a maximum of 15 m radius, thus enabling the detection of all the plot-wise trees including their branches, and exported from Cyclone at their initial accuracy.

Overview of the diameter modelling
Diameter modelling is presented as a simplified flowchart in Fig. 2. Processing starts from a co-registered and preprocessed point cloud generated by terrestrial laser scanner, which is used to locate the stems and model their structure with cylinders, until all the stems on the plot are detected and their respective diameters are calculated. Further, dimensions of these cylinders are utilized in the complementary steps, which include extracting stem points from the cloud, dividing them into thin horizontal slices, re-evaluating slice diameters by circle fitting, and constructing a taper curve model. Finally, the diameters of these two approaches are compared to the field-measured values at the heights of 1.3 and 6.0 m, later referred as d1.3 and d6. The complete processing chain is described in details in Sections 2.4 and 2.5.

Stem extraction and cylinder-based stem modelling
Stem extraction and cylinder modelling was the first phase in the stem diameter processing, and performed using MATLAB 2016b software. The processing of cylinders was based solely on the points' location (x-y-z) without utilizing any additional information such as scan numbers, intensity values or color information.

Point cloud filtering
Before extracting the stems from the TLS point cloud, the initial coregistered point cloud was filtered. The aim of the filtering was to decrease the unnecessary points by removing points from the ground, reduce the point density and exclude low-density regions. To remove points from the ground and understory layers but enable measuring tree height at later stage, ground level was defined using a grid of minimum z-coordinates in 0.5 m × 0.5 m squares. The lowest 20 cm layer was expected to be consisted of ground returns and low vegetation, and removed from each square. Furthermore, if squares contained points up to two meters height but having at least one meter empty space above them, thus indicating only the existence of herbaceous plants, shrubs or saplings, such points were removed as well. Next, the density of the point cloud was decimated by partitioning it into small cubes, having an edge length of 5 mm, and keeping only one point per each cube. Finally, filtering based on point density was applied to the decimated point cloud by partitioning it into 5 cm cubes, and totally removing points inside those cubes that contained less than three points. As a result, point cloud was markedly reduced, thus enabling faster processing while not significantly affecting to its stem-related information content.

Stem extraction from TLS data
Surface patch generation: Initially, filtered point cloud was first partitioned into small sets called surface patches. The patches were generated using cubical partition of the point cloud and constant-size neighborhood balls with random but distance-restricted centers, which is explained in more details in Raumonen et al. (2013). The neighbors, surface normals and heights of the patches were also determined ( Fig. 3). Diameters of these patches varied randomly between 15 and 30 cm, and their neighbor relations were used to determine connected patch clusters.
Locating potential stem sections: Stems were located using the patches by an iterative process, which was based on detecting nearly horizontal surface normals. Iterations were performed using different elevations and degrees of horizontality of the normals to locate as many stems as possible, and to bypass stem sections affected by occlusion effects. This process started by first selecting one-meter layers of patches between 1.5 and 2.5 m above the ground, and continued by increasing the elevation in 0.5 m steps up to 8 m. For each one-meter layer, patches with normal angle of less than 15 degrees from the horizontal were selected and connected components of at least eight patches were determined as potential stem sections, which were extended into whole stems as explained in the next two paragraphs below. At this step, only patches not yet assigned to any of the previously detected stem were used. After the components with less than 15 degrees normals were tested, angles of less than 20 degrees, and lastly 25 degrees, were tested as well. Similar idea of using nearly horizontal surface normals to locate potential stems from point clouds has been used various earlier studies such as Raumonen et al. (2015), Holopainen et al. (2013) and Liang et al. (2012).
Acceptance of stem sections and their cylinder modelling: Next, at each iteration step the selected stem components were checked if they are part of a stem, and further modelled as cylinders. purpose, vertical point density, i.e., number of points in vertical columns within 2 cm × 2 cm squares was calculated. This density was used to weight the cylinder fitting based on the least squares fitting, assuming that points deriving from tree stems will have more points directly below or above them compared to branches and outliers. The first modelled cylinder fit was used to filter out the furthest 30% of the points, expecting that they may not be part of the stem. Then, second density-weighted cylinder fitting was applied using the filtered point data to gain a better estimate of the stem axis and radius. Further, this refined cylinder was used to check if there are points remaining inside the expected stem, assuming that the number of such points should be low. To finally accept the cylinder as a part of a true stem, the number of points inside the given cylinder excluding 3 cm tolerance distance (2 cm for cylinders having a radius of less than 80 mm) was defined as not exceeding 10. Expanding the stem sections into whole stems: The accepted stem components were expanded into whole stems, or at least as to be as complete as allowed by the data. The first fitted cylinder between 1.5 and 8 m height above the ground was used to select all the patches within a 3 m horizontal distance from the cylinder and expected to contain the whole stem. After that the points, which were within 4 cm distance from the cylinder in horizontal direction and within 30 cm distance above or below the cylinder ends, were selected as parts of the stem point cloud. Processing was continued by defining the whole stem, points and cylinders, first towards the bottom, and thereafter towards the top. In the downwards processing, the one-meter section below the previous cylinder/section was selected, and points within the upper cylinder radius plus 5 cm were used for cylinder fitting to define the stem part of the section as described above. In the upwards processing, the procedure was similar except that selection of the points was based on the cylinder radius below plus 2 cm. When approaching higher parts of the stem, occlusion started to deteriorate the quality of the cylinder fitting. When no new upwards sections were found anymore, or the detected axis direction between the two sections to be compared deviated more than 25 degrees, the point cloud quality was regarded as being too bad to continue any further. Finally, the base of the stem was defined as being the lowermost 3 cm section.
Improved stem extraction based on branch segmentation: The surface normal and cylinder-based stem detection explained above defined the stems only approximately. To define them more accurately, separate stems from the branches and make more accurate cylinder fitting, stems were defined again based on a branch segmentation method. During this step, the extracted stem point clouds were first expanded to include a few layers of neighboring patches, thus including both stem surface and beginnings of the nearby branches. Then, branch segmentation and cylinder modelling for each stem was performed based on the modified quantitative structure model reconstruction method as presented in Calders et al., (2015). This included covering the point cloud again with surface patches by using small patch diameter (1-2 cm) to allow accurate branch separation from the stem. After this, the neighbor relations of the sets were modified to make the stem connected from the base to the top. Finally, the stem segments were divided into short sections, and final cylinders were fitted with approximate relative length (length/radius) of 20, which was assessed to provide a suitable balance between robustness and accuracy.

Extraction of the cylinder-based variables
Each plot was modelled at least five times, as it was known that there is a small amount of random variation between the models, which affects stem detection, definition of their components, and further measurements of the cylinder diameters. Stems less than 25 cm apart in the distinct modelling attempts of the same plot, as measured at the height of 1.3 m, were interpreted to represent the same individual. This distance was found as a good compromise to distinguish single trees without unintentionally merging separate stems together. Furthermore, to make the extracted stems to correspond the field data, those stems modelled to have d1.3 of less than 45 mm were removed. Cylinder models were used to extract cylinder-based d1.3, d6 and tree (stem) height to be compared with later circle-based evaluations. Of these variables, d1.3 and d6 were however not extracted directly from the initial model, but by fitting separate cylinders specifically between 1.1-1.5 m and 5.8-6.2 m heights. This was considered as providing a better estimate of the potential accuracy compared to the initial and often longer cylinders. The height was calculated as the total length of the cylinder structure. Final cylinder-based d1.3, d6 and height values for a given tree were calculated as an averaged value of all the stems representing the same individual. Cylinder structure (i.e. dimensions, inclinations and coordinates) to be used in the further cloud-based measurements, however, was based on the longest stem among the models of the same individual. If a stem was detected only once but missing in all the remaining modelling results, it was discarded and considered a probable false detection of a non-tree structure.

Complementary modelling steps
After construction of the cylinder models, additional fitting procedures were performed using the R statistical software version 3.3.1. Similarly to MATLAB modelling, the decimated point cloud, containing only one point per 5 mm cube, was used to reduce redundant information and facilitate faster processing. Regardless of the co-registration, however, the initial scan numbers as well as backscattered intensity values were stored as point attributes, thus enabling their use in the analysis. It should be noticed that the parameter values and limits presented in the following subchapters were primarily decided based on multiple tests, and found to be suitable for the used TLS material. They may not, however, be optimal for other types of forests, field setups, or scanner settings.

Evaluating tree dimensions and slicing the stem
Ground level and treetop evaluation: For each stem, TLS points in the vicinity of the modelled stem were first extracted. As the cylinder model was known to be often too short due to occlusion effects of branches as well as ground vegetation, both ground level and treetop height were re-evaluated by vertical dimensions of the point cloud. Ground level estimation was based on extending the lowermost cylinder downwards, extracting TLS points within 50 mm horizontal distance from the cylinder edge, and defining the new ground level by the uppermost 10 mm vertical section with no points. Maximum allowed downwards extension was decided to be 30 cm, corresponding to the maximum height of the herbaceous plants commonly found on the area, and interpreting any potential stem extensions below that to be likely originating from non-stem structures. The treetop height was correspondingly evaluated by selecting the potential TLS points above the modelled stem, and observing their vertical continuity. The applied selection distance was based on the cylinder-derived d1.3 instead of any fixed distance, given that smaller (i.e. shorter) trees were more likely to be located at the understory layer and having branches of other trees above them. More precisely, TLS points were extracted using a horizontal distance of 2 • d1.3 on top of the uppermost cylinder, and the lowest vertical section of at least 3 • d1.3 with no points was used to define the treetop.
Stem slicing: Next, the extracted point cloud around the respective stem was cut into thin slices, which were used for circle fitting as planar projections (x,y). Distinct slices were initialized at every 20 cm height along the whole cylinder model by selecting TLS points within 25 cm horizontal distance from the respective cylinder edge, and straightening up the stem if the z-axis indicated an inclination of 10 degrees or more.
Then, centered at the respective height, the actual slice was constructed by either reaching a limit of 5,000 TLS points, or alternatively finalizing the slice at a vertical maximum distance of ± 10 cm. This approach resulted in shorter slices along the well-scanned stem parts but prevented too long vertical variation to exist even in the lack of points. Furthermore, point clouds of single stems were also plotted on the screen from two different sides (according to x,z and y,z axes) and saved as image files for later inspection.

Slice processing
Diameter evaluation with co-registered points: For each slice, applicability of the cylinder diameter was initially evaluated against the co-registered points. Instead of approving the fit by any given number of associated points, which may depend on various reasons related to nearby obstacles, scanner settings, preprocessing procedures or surface properties, a method calculating the continuity of the circle edge was applied for the evaluation. The edge was defined in proportion to the respective cylinder diameter (d c ) and set to have a tolerance of ± 0.0375 • d c , respecting the fact that larger trees tend to have more uneven surface and often slightly less circular shape. Then, the extracted circumference was divided into 100 equal sections (i.e. each covering 1% of the total edge length), or into 50 sections for diameters less than 100 mm due to relatively shorter circumference length and proportionally longer point-to-point distances. Continuous edge was defined as at least three consecutive sections having edge points, thus filtering out the effects of single outliers. Further, the maximum allowed number of points inside the circle (minus edge) was set to be 100 • d c , where d c was expressed in meters, as only few points were expected to be found inside the modelled stem. Finally, if continuous edge covered more than half of the circumference, and points inside did not exceed the maximum limit, the initial evaluation was approved. Additionally, given that cylinder was at correct location but slightly under-or overestimated the stem diameter, 5% and 10% larger as well as 5% and 10% smaller diameters were tested with similar approach.
Circle fitting (co-registered points): If the initial evaluation was approved, the final slice diameter was modelled using the co-registered point cloud and random sample consensus (RANSAC) algorithm (Fischler and Bolles, 1981). RANSAC tests a selected geometric shape against point data by picking a small sample from the full set of points, fitting the focused shape using the sample, and evaluating the result with the initial data. A sample of three points is sufficient to initiate a solution for a circle arc in two-dimensional space, which can be evaluated e.g. by the number of points in the vicinity of the arc. If the result is accepted, smoothing technique such as least squares fitting can be used to get an improved estimate of the circle parameters. RANSAC requires multiple repetitions and is therefore computationally heavy, but it has been found useful in TLS-based data analysis given its low tendency to be affected by outliers (Olofsson et al., 2014, Schnabel et al., 2007, Wang et al., 2016a, 2016b. RANSAC also allows for estimating a sufficient number of iterations (k) with Eq. (1) (Fischler and Bolles, 1981): where z is the desired probability that at least one of the random selections is error free, w is the probability that a randomly picked point belongs to the modelled shape or locates within a defined tolerance limit from it, and n is the number of sampled data points used in the shape construction. Prior to RANSAC iterations, 25% of the lowest intensity points were filtered out. In general, intensity values depend on various factors such as distance, wavelength, incidence angle and the material itself, but woody surfaces have generally higher spectral response compared to e.g. leaves of conifer needles (Côté et al., 2009, Kaasalainen et al., 2011. Number of RANSAC iterations was set to be 500, decided based on Eq. (1) by using values z = 0.999, w ≈ 0.25 and n = 3. Iterated circles deviating more than 80 mm from d c were not accepted, as cylinder error compared to the field-measured values did not generally exceed this limit. Minimum diameter of 45 mm was applied, corresponding to the smallest diameter measured in the field, as evaluation of very small radii was expected to be uncertain. In addition, circles located further away than 50 mm from the initial cylinder center, or having too many points inside (greater than100 • d c ), were rejected. Among the accepted circles, the final slice diameter was decided based on the highest number of edge points, and the proportion of continuous edge was calculated for later taper curve modelling. Against the normal RANSAC procedure, however, least squares fitting was not applied as the solution without additional smoothing was found to be more robust against potential outliers.
Circle fitting (single-scan points): If the initial cylinder-based evaluation using co-registered points was not approved, reasons included stem deformations, occlusion effects or co-registration flaws. As however partial circle edges were expected to be found from the distinct scans, the analyzed slice was split into data from single scans and a RANSAC-assisted fitting procedure similar to above was performed for each of them. As the actual stem was potentially further away from the initial cylinder, approved circle distance from the respective cylinder center was increased to 200 mm, and the number of RANSAC iterations to 2000 (i.e. z = 0.999, w ≈ 0.15 and n = 3). However, as single scans were only able to see the stem from one side, the allowed number of points inside the stem was decreased to 50 • d c , and the total proportion of the continuous edge length was multiplied by two to better correspond to the co-registered data. Resulting circle iterations with less than 25% proportion of continuous edge length were discarded due to their expected poor quality. As the number of edge points was not comparable between the distinct scans, the final slice diameter and continuous edge proportion were calculated as a mean value weighted by edge proportions of the single measurements. If no accepted diameter solution was found after the RANSAC evaluation of single scans, the analyzed slice was discarded and no diameter was calculated.

Calculating final d1.3 and d6 values
As occlusions and outliers affect the single circle fitting results, d1.3 and d6 values were calculated using a taper curve model. A taper curve indicates the relationship between the stem diameter and height, which is often modelled using spline functions and potentially calibrated according to the focused species or stand-related stem form correlations (Kilkki and Lappi, 1987;Koskela et al., 2006;Laasasenaho, 1982;Lappi, 2006). Both quadratic and cubic spline functions for taper curves are used. Of these, cubic estimation allows creating more flexible curve but on the other hand is more susceptible for under-and overestimates if the used data is not consistent, and does not necessarily preserve monotonicity (Lahtinen and Laasasenaho, 1979;Lahtinen, 1988).
Prior to taper curve calculation, potentially overestimated slice diameters were removed as they may have had significant impacts on the resulting curve. Similarly to Liang et al. (2014), a strategy of comparing a diameter to the mean value of three preceding (i.e lower in the stem) measurements was applied. If a diameter larger than the calculated reference value was found, and it simultaneously had a smaller continuous edge proportion than any of the reference diameters, it was replaced by the previous diameter measurement. Furthermore, treetop height was added to the curve calculation with a diameter of 0 mm and continuous edge proportion of 75%, therefore guiding but not forcing the diameter calculation in the upper part of the stem. After these procedures, a taper curve was calculated using cobs package in R (Ng and Maechler, 2017). The calculated curve was estimated using 2 nd degree spline function, which was constrained to be decreasing, weighted by the observed continuous edge proportions, and set to automatically select a suitable lambda value i.e. the penalty parameter. This resulting spline model was then applied to calculate the final d1.3 and d6 values.

Field data and TLS tree detection
Altogether 505 trees were measured in the field, ranging from 8 to 52 individuals per field plot. Their d1.3 varied between 45 mm and 468 mm, having average diameter of 154.1 mm and median diameter of 134 mm. In total, d1.3 was recorded from all the trees, d6 from 209 trees and height from 242 trees, according to the angle count based measurement strategy as specified before. In terms of species composition, 42.8% of the trees were Norway spruce (Picea abies), 29.3% Scots pine (Pinus sylvestris), 18.4% Silver birch (Betula pendula), 8.3% Downy birch (Betula pubescens) and 1.2% other broadleaved species, respectively.
When the field-measured 505 trees were checked against the cylinder models and point cloud data, one tree appeared to be fallen down prior to laser scanning and one could not be matched with any detectable tree near to the given location, therefore decreasing the number of confirmed field-measured trees to be 503. Of these trees, 442 (87.9%) could be matched with the final TLS-modelled stems, which were then used to compare the d1.3 (n = 442), d6 (n = 203) and h (n = 234) values. Majority of the non-detected 61 trees were relatively small, having an average d1.3 of 69.6 mm and only 8 of them exceeding a d1.3 of 100 mm. Most of them were young spruces with dense branches near to the ground, small trees which were growing in groups, or trees with relatively high occlusion rate due to various obstacles. In addition, a few non-detected trees had a broken trunk or they were at a close proximity of another tree, which potentially caused modelling attempts to either discard them or merge as a single detection. Majority of the matched stems (88.9%) were extended downwards from the modelled cylinders prior to circle fitting procedures with an average value of 7.7 cm. As calculated using the final circle-based taper curve, the effect of this extension on the calculated diameter value was on average 0.78 mm at 1.3 m and 0.66 mm at 6.0 m height levels.
In addition to the field-detected trees, 32 trees were detected from the TLS data, but were not recorded in the field. Of these, 18 were existing trees which could be confirmed from the point cloud, but had been excluded in the fieldwork mainly because of having measured as too thin, or being slightly out of the plot radius. The remaining 14 TLSdetected trees were no actual trees but rather included other vertical objects (e.g. a partially broken branch of another tree, twisted into a vertical position), markedly erroneous measurements (e.g. two small trees detected as a single structure) or various non-tree structures, which could not be verified from the point cloud. Most of these could be immediately detected as non-trees based on their previously stored stem profile images, and they contained rarely a single slice, which would have had a continuous edge proportion exceeding 60%. Compared to any real and even partially occluded tree, having normally several slices with over 80% edge detection, this was a notable difference and offering a relatively accurate automatic exclusion method if needed.

Modelling results
In general, diameter calculations are relatively accurate for both using cylinder modelling as well as when processed through further taper curve assisted circle fitting procedures (Fig. 4). However, differences between the two are distinctive. In terms of mean absolute error of all trees (Fig. 4 a), circle fitting results indicate an error of 5.25 mm (3.2%) for d1.3 and 6.45 mm (3.6%) for d6, whereas the corresponding errors for cylinders are 6.18 (3.7%) and 8.52 mm (4.8%). Regarding to bias, both circle-calculated d1.3 (−0.43 mm/−0.3%) and d6 (+1.63 mm/+0.9%) are statistically unbiased at the confidence level of 0.95. Respective biases for cylinders indicate 3.97 mm (2.4%) underestimation for d1.3 and 4.14 mm (2.3%) underestimate for d6, which are considerably larger compared to circle fitting results, but not T.P. Pitkänen et al. ISPRS Journal of Photogrammetry and Remote Sensing 147 (2019) 294-306 exceeding statistically significant limit. For RMS errors, values for circle fitting are 7.34 mm (4.4%) for d1.3 and 8.44 mm (4.8%) for d6, and for cylinders 9.76 mm (5.9%) and 13.09 mm (7.4%), respectively. It should however be noticed that also the field measurements have their own sources of errors as discussed later, therefore rather providing comparatively accurate reference estimates than truly error-free diameter values. When observing trees of different sizes (Fig. 4 b-e), mean absolute and RMS errors of the d1.3 increase according to the size of the tree, but relative errors are rather decreasing. Circle fitting produces smaller Fig. 4. Errors related to taper curve assisted circle fitting method (Circle) and cylinders fitted at a specific height (Cylinder). Errors have been calculated for d1.3 and d6 using mean absolute error (MAE), bias and root mean squared error (RMSE), and number of included trees to calculate the given statistic has been indicated on top of the respective figure. Relative errors, as compared to the mean value of the corresponding field measurements, are indicated in brackets. For bias, an asterisk (*) has been added to statistically significant deviations at the confidence level of 0.95, i.e. values larger than ± 1.96 * σ M . errors and lower bias in almost all the cases, whereas cylinder model tends to provide larger underestimates when the size of the tree increases. In terms of d6, circle fitting results outperform cylinder-based results with exception of the largest size class where circles produce larger mean absolute error and bias, the latter however not being statistically significant. Difference between the cylinder and circle fitting models is particularly obvious with d6 measurements of small trees (Fig. 4 b), however, it should be noticed that number of measurements in this category is only six, therefore not providing representative figures. Regarding to the different tree species (Fig. 4 f-h), the most notable deviation is the poor performance of circle fitting results for spruces as compared to cylinders, given particularly their high although not statistically significant d6 bias. This bias is heavily influenced by a small number of large individuals at a single plot, for which the circle fitting is producing d6 overestimates between 15 and 18 mm, but cylinder errors remain between 7.1 and 11.3 mm. Circle fitting results, however, are comparatively better for pine and birch. Further, results were checked against different wind speeds (Fig. 4 i-j) to observe potential wind-induced registration errors. In this context, the circle fitting seems to perform slightly more consistently, but no major differences are to be found. However the highest wind speed measured during the scanning was 6.5 m/s, which may not yet be enough to show major differences at or below the height level of 6.0 m.
At the level of single slices, most circle diameter measurements were relatively accurate and could not be considerably improved by minor methodological changes. Fig. 5a (spruce, slice taken at a height of 8.0 m, circle diameter 153 mm) shows a typical situation when the measured tree is well scanned and round in shape, thus leading to a successful measurement regardless of outliers. Dividing the slice into single scans assisted in measuring diameters like Fig. 5b (pine, h = 15.0 m, d = 181 mm), where the data was affected by co-registration faults due to wind. Diameters were generally following the tree surface even in case of slight stem deformations or bark roughness (Fig. 5c, pine, h = 4.0 m, d = 298 mm), but this often resulted in a larger bias depending on the actual stem form and applied measurement direction.
Larger diameter errors were resulting from obvious stem deformations (Fig. 5d, birch, h = 0.8 m, d = 296 mm) or lack of continuous circular shapes to be found (Fig. 5e, spruce, h = 6.6 m, d = 53 mm), which were primarily connected to occlusion effects, but often further impaired by co-registration defects and poor initial cylinder-based measurement. Some less frequent reasons for unsuccessful measurements were found as well, such as Fig. 5f (grey alder, Alnus incana, h = 1.0 m, d = 128 mm). In that particular case, the stem had short slanted sections, which were not indicated by the cylinder inclination, resulting in thicker point concentrations along the stem edge and difficulties in detecting a continuous circle within the set tolerance limits. This problem could have been avoided by increasing tolerance limits, but such an action appeared to have negative impacts to cases when a slice had a considerable amount of outlier points just outside the stem edge.
When the modelled tree heights were compared to the field-measured values, cloud-adjusted heights prior to circle fitting procedures resulted in an absolute error of 1.23 m (6.3%), bias of -0.89 m (-4.5%) and RMSE of 2.1 m (10.7%). The respective values for cylinder model without height adjustment were 4.19 m (21.2%), −4.16 m (21.1%) and 5.09 m (25.8%). Relatively high and almost equal mean absolute error and bias for the cylinders indicate that the height was practically always underestimated, which derives from the fact that a reliable cylinder model can only rarely be constructed until the very top of the modelled tree. Further cloud-based calculations, however, provide a relatively accurate estimate that corresponds well with the field-measured values. In addition, it should be noticed that measuring heights in Fig. 5. Examples of circle fitting results (black dashed line) against TLS points (colored by separate scans) of a single slice, including a typical stem slice and successfully fitted circle (a), stem section affected by wind (b), stem section with considerable bark roughness (c), strongly deformed and non-circular stem (d), lack of TLS points and deficient fitting result (e), and an erroneous result due to slanted stem (f). Coordinates indicate meters relative to the plot center (0,0). T.P. Pitkänen et al. ISPRS Journal of Photogrammetry and Remote Sensing 147 (2019) 294-306 the field is not always too straightforward, and tends to result in slight underestimates as well due to occlusion effects and consequent difficulties in detecting the respective treetop.

Accuracy of the modelling results
Success and reliability of detecting stem diameters using TLS data has been tested and discussed in several recent papers. In one of the earlier studies, Henning and Radtke (2006) gained an average error of less than 10 mm at the lower parts of loblolly pines (Pinus taeda), which were planted in rows with regular spacing. De Conto et al. (2017) reported RMS errors of 10.6-32.1 mm for Scots pine (Pinus sylvestris) and 11.9-35.8 mm for Norway spruce (Picea abies) using a small set of trees, performing scans specifically designated for single trees, making field measurements from various heights and testing several different modelling methods. Bias in their study was generally between 7 and 13 mm. Olofsson et al. (2014) used circle-based RANSAC method and a plotwise single scan setup (r = 20 m) to detect stem diameters of various tree species, resulting in a breast height level RMSE of 28 mm and bias of 6 mm with an outlier-excluded set of trees. Liang et al. (2014) collected data from nine sample plots (r = 10 m) with seven separate scans, of which two to three nearest ones to a particular stem were used for cylinder modelling and further calculation of the stem curve. In their study, bias was 1.5 mm and RMSE 11.3 mm for a combined set of Scots pines and Norway spruces. Srinivasan et al. (2015) reported of RMS errors between 18.3 and 21.3 mm when scanning was performed on two study sites dominated by post oak (Quercus stellata) and loblolly pine (Pinus taeda) at different times including both leaf-on and leaf-off, and using different setups (both single and multi-scan). Maas et al. (2008) were also comparing several plots with different scanning procedures, equipment, tree species compositions and scanning dates, concluding to RMS errors between 14.8 and 32.5 mm for DBH measurements. Further, Pueschel et al. (2013) compared various earlier TLS studies with circle-or cylinder-based diameter measurements at breast height, reporting from RMS errors of 14.8 to 42 mm and biases of 3.0 to 36 mm.
Although direct comparison between studies is difficult due to numerous differences in e.g. scan setup, scanner settings, scanning conditions, tree species, stem density and shape, plot size, co-registration details, level of automatization and applied methodology to detect the final diameters (Liang et al., 2018), the results presented in this manuscript are among the most accurate ones. Given the limited data used in this study and lack of benchmarking the results on any other data set, however, further tests are needed to ensure the applicability of the used methods in different conditions. Moreover, fieldwork is also subject to inaccuracies both in terms of measuring a correct value as well as selecting a consistent height, which affect the reported modelling errors. In one of the few studies made on this topic, Päivinen (1987) compared the results of 14 experienced fieldworkers who measured the DBH for a set of trees using a manual caliper. In terms of the measurement height, he concluded that the standard deviation in selecting the breast level was 3.6 cm in vertical direction while single deviations were over 15 cm. If the height level was standardized and marked, standard deviation in the diameter measurement in one direction (against the plot radius) was 3.5 mm. In another study, Elzinga et al. (2005) estimated that diameter errors of at least 5% may be expected to occur in at least 5% of the measurements, particularly for smaller trees. Therefore, errors reported for TLS-based measurements are not only deriving from scanning per se or analyses conducted to calculate the diameters, but also including unknown variance between field measurements and true diameters. Further, this leads to the outcome that even if TLS-based methods would be perfectly successful in measuring the diameter, errors would be unlikely to drop very close to zero. Bias, however, should not be too large if field measurements are not expected to have any systematic errors. In addition, it should be noted that the field-measured diameters used in our study were based on a fixed direction against the plot radius, therefore being sensitive to partial stem shape irregularities. In this matter, a better strategy would be to field-measure diameters from at least two perpendicular directions, or to use a diameter tape.
Tree height was not the primary target to be modelled in this study, but the gained cloud-based estimates (mean error 1.23 m, bias −0.89 m and RMSE 2.1 m) are also encouraging. For example, Srinivasan et al. (2015) reported a height bias of +0.3 m and RMS error of 1.51 m for a subset of the studied trees, whereas Olofsson et al. (2014) had very small bias but a RMSE of 4.3-4.9 m. Liang et al. (2016) summarized several earlier studies, which had resulted in a height bias between −0.3 and +0.6 m, and RMSE between 0.8 and 6.5 m. In our study, however, part of the plots were highly occluded and a number of understory trees were measured, therefore making it crucial to separate their tops from the branches of upper trees instead of only recording the uppermost TLS points. This indicates that a measurement strategy based on empty space on the top of the tree, which is adjusted with the detected d13, can be a workable option when detecting trees of various heights and layers.
Differences between the final slice-based circle fitting results and the cylinder measurements are clear, but it should be stressed that these two approaches have certain important dissimilarities. In principal, circles can be regarded as flattened cylinders, which should be comparable to their three-dimensional counterparts at similar heights. In practice, however, vertically shorter circle-based measurements with assistance of a spline-based taper curve have potential to neglect obvious over-and underestimates, bypass partial occlusions with less problems, and balance out local stem deformations deriving from e.g. to branches, their bumps or bark roughness Saarinen et al., 2017), which may cause differences compared to field measurements. Further, selecting the longest stem among several cylinder models and capabilities to use data from single scans instead of relying only on the co-registered point cloud are decreasing sensitivity to defects resulting from co-registration inaccuracies or modelling faults. Potential height adjustment of the stem base gives slight advantage for the circle-based approach in finding the correct measurement levels, although the effect is relatively minor e.g. compared to the improvement in the modelling bias. The measured average vertical difference of 7.7 cm between the cylinder model and the adjusted stem base is at a similar level to observations of Mengesha et al. (2015). These advantages are however gained at a cost of applying a tedious and computationally more intensive approach, which relies on the cylinderbased a priori information, but potentially lacks the robustness of the cylinder model at the level of a single measurement. Additional processing caused by the complementary steps will at least double the time used for modelling, but this depends highly on the number of TLS points, number and dimensions of stems detected on the plot, and the potential need to use single scans for diameter evaluation instead of the co-registered data. The most influential single factor affecting the computational heaviness is the number of selected RANSAC iterations, which further determines the likelihood of finding the best possible fit. Therefore, there is no single optimal solution to suit for all the cases, but rather several alternatives which should be selected based on the data, modelling needs and available resources.

Sources of errors
Errors found in the circle fitting were primarily resulting from two reasons: either the tree was not round enough to be detected using a circular shape, or there was too large amount of outlier points in relation to the actual stem edge. In terms of selecting an appropriate shape, using a circle or circular cylinder provides a robust way of modelling the stem which can also handle outliers and gaps reasonably well. Most other shapes are more sensitive to data quality and should T.P. Pitkänen et al. ISPRS Journal of Photogrammetry and Remote Sensing 147 (2019) 294-306 therefore be used with caution (Åkerblom et al., 2015). Instead of a circle, ellipses have been sometimes found useful e.g. in detecting diameters of slanted trees or modelling trunk parts without obvious occlusion effects (Belton et al., 2013, Bu and. Elliptic shapes increase flexibility in stem modelling, but their usability is significantly deteriorated in conditions when occlusions and outliers are present, thus leading easily to larger errors. Use of ellipses could be justified if a number of modelled trees are known to have a distinctively elliptical shape and the absence of outliers can be confirmed. Circular form however was found the best option to be used in our study, and slanted trees did not pose a particular problem as most of their effects were avoided with assistance of the cylinder model. Outlier points not belonging to the stem, when co-occurring with occlusion effects, appeared to be the main reason to hinder stem detection and produce erroneous measurements. As long as the stem shape was relatively round and its edge was detectable through TLS point concentrations in the slice (see Fig. 5a), the existence of outliers outside of the stem did not significantly impair the calculations. With help of splitting the slice into single scans, the effects of misplaced TLS points due to co-registration faults could also be often neglected (Fig. 5b). Difficulties however started to appear particularly in upper parts of the stem where the point cloud was sparser and outliers occupied a larger proportion of all TLS points (Fig. 5e). According to Henning and Radtke (2006), Olofsson and Holmgren (2016) and Mengesha et al. (2015), occlusion effects are affecting the measurements from 10 to 15 m height upwards, and our results are supporting these observations. Problems were also more distinctive for Norway spruce than the other dominant species, given its tendency to have relatively dense branching structure throughout the whole stem. To some extent and at a cost of increased computational efforts, a larger number of RANSAC iterations could have resulted in better stem detection. One of the underlying problems, however, was in the circle evaluation method related to the co-registered data, which emphasized the number of points. While this method was found to be relatively robust and providing correct measurements in most situations, it was more prone to misdetections when there were denser point concentrations, which did not belong to the stem. Further fine tuning with e.g. intensity values or by setting stricter limits for RANSAC iterations could improve the results in these conditions.
To further improve the results and reduce excessive number of outliers, scanning design and fieldwork conditions are in a key role. Making more scans will help to avoid occlusion effects, but from practical point of view it may not be always possible. Therefore, deciding optimal places for stations with respect to the scanned trees is essential, which would probably benefit from increased flexibility compared to the relatively static scanning design used in this study. Regarding to the fieldwork conditions, both weather and phenology have importance for the quality of the results. While scanning during wet weather is normally avoided due to altered laser pulse transmission and scattering properties (Wilkes et al., 2017), similar obvious restrictions are not applicable for windy conditions. A maximum wind speed of around 5 m/s has been mentioned in a few studies to assure good data quality (Dassot et al., 2011;Seidel et al., 2011;Seidel et al., 2012), but splitting the co-registered cloud into single scans should give more flexibility in this matter. Wind-induced problems are however most prominent in the upper parts of the stem (Vaaja et al., 2016), and suggesting any specific limits to restrict scanning would require analysis from more diverse conditions. Data for this study was collected during leaf-on conditions, which increases occlusion effects of deciduous trees compared to leaf-off time. The period when deciduous trees are leafless but there is no snow cover, is normally however quite short in the boreal zone. Therefore, restricting station-based TLS scanning to leaf-off conditions would pose serious limits for data collection and areal coverage thus not being applicable for extensive studies. To some degree and at a cost of slower operation, this could be improved by using a scanner capable of recording multiple echoes. A more feasible alternative for efficient data collection, however, could be offered by mobile scanning systems, which are based on an intertial measurement unit and typically tied to a satellite navigation system. While they can significantly accelerate and facilitate scanning efforts, their potentially associated problems include lower accuracy compared to stationary TLS data, instability when operated under dense forest canopies, and a limited range to detect the backscattered signals (Bauwens et al., 2016;Liang et al., 2016). Therefore, relying on the mobile scanning devices may not presently offer sufficient accuracy for many applications, but this topic is without doubt subject to further developments.

Conclusions
In this study, a processing chain was tested where tree stems were first modelled as cylinders, which were then used as a priori information for further circle fitting procedures, taper curve estimation and diameter calculation. According to the results, the presented approach was capable of improving the initial cylinder measurements and reducing their diameter errors, including minor co-registration faults and wind disturbances. These adjustments can help to extend the range of suitable scanning conditions both in terms of forest types and weather conditions, and augment the potentialities of TLS-based analyses in forest inventories. Further, the improved taper curve has potential to increase the accuracy of volumetric calculations as well as facilitate the separation of timber grades. Occlusion effects and large quantities of outlier points, however, are still posing problems that require further methodological development and optimization of the scanning design. While these results are encouraging, similar analysis should be extended to a larger geographical area, different forest types and more variable conditions, as well as tested with different TLS technologies, which would assist in defining best practices and applicable preconditions for successful diameter measurements.