- Email: [email protected]
Performance of stem denoising and stem modelling algorithms on single tree point clouds from terrestrial laser scanning
Computers and Electronics in Agriculture 143 (2017) 165–176
Contents lists available at ScienceDirect
Computers and Electronics in Agriculture journal homepage: www.elsevier.com/locate/compag
Original papers
Performance of stem denoising and stem modelling algorithms on single tree point clouds from terrestrial laser scanning
MARK
⁎
Tiago de Contoa, , Kenneth Olofssonb, Eric Bastos Görgensc, Luiz Carlos Estraviz Rodriguezd, Gustavo Almeidad a
Swedish University of Agricultural Sciences (SLU), Sundsvägen 6, 230 53 Alnarp, Sweden Swedish University of Agricultural Sciences (SLU), Skogsmarksgränd, 907 36 Umeå, Sweden c Universidade Federal dos Vales do Jequitinhonha e Mucuri (UFVJM), Campus JK, Alto da Jacuba, 39100-000 Diamantina, MG, Brazil d Escola Superior de Agricultura Luiz de Queiroz – Universidade de São Paulo (ESALQ-USP), Avenida Pádua Dias, 11, 13418-900 Piracicaba, SP, Brazil b
A R T I C L E I N F O
A B S T R A C T
Keywords: LiDAR Cylinder/circle fit Tropical and boreal tree species Robust estimation RANSAC
The present study assessed the performance of three different methods of stem denoising and three different methods of stem modelling on terrestrial laser scanner (TLS) point clouds containing single trees – thus validating all tested methods, which were made available as an open source software package in the R language. The methods were adapted from common TLS stem detection techniques and rely on finding one main trunk in a point cloud by denoising the data to precisely extract only stem points, followed by a circle or cylinder fitting procedure on stem segments. The combination of the Hough transformation stem denoising method and the iteratively reweighted total least squares modelling method had best overall performance – achieving 2.15 cm of RMSE and 1.09 cm of bias when estimating diameters along the stems, detecting 80% of all stem segments measured on field surveys. All algorithms performed better on point clouds of boreal species, in comparison to tropical Eucalypt. The point clouds underwent reduction of point density, which increased processing speed on the stem denoising algorithms, with little effect on diameter estimation quality.
1. Introduction Terrestrial Laser Scanning (TLS) sensors are increasing in popularity in the Forestry sector, due to the high resolution and precision of the three-dimensional information provided at the sample plot level, making further 3D modelling and tree reconstruction possible and allowing foresters to extract dendrometric variables with high accuracy from point clouds (Liang et al., 2016). In order to apply the TLS technology in forest inventory at its full potential, robust tool sets to extract useful information from TLS point clouds, and that perform well on a variety of circumstances are required – especially stem denoising and stem modelling, since stem volume is the most targeted variable in commercial forest management. TLS is a technique that could be used to develop locally adjusted stem taper functions (Trincado and Burkhart, 2006; Sun et al., 2016). Routines for processing TLS point clouds containing tree data have been tested, with applications ranging from isolation of single trees in scanned plots (Liang et al., 2012; Olofsson et al., 2014) to segmentation of tree components (e.g. branches, stem and canopy) (Raumonen et al., 2013; Hackenberg et al., 2014).
⁎
In order to automate the process of isolating tree components in point clouds, general assumptions need to be made about a tree’s geometrical structure. A reliable algorithm for stem isolation, thus, requires a level of generalization that allows it to capture the essential information of a variety of tree shapes in a point cloud. Some commonly adopted assumptions are the vertical orientation of the bole (approximately orthogonal angle with the ground), approximately circular stem cross sections and approximately cylindrical stem segments. Those assumptions are the core of many already implemented algorithms, but a key issue must be tackled beforehand: fitting circles and cylinders works well on point clouds with the stem clearly visible, but can’t be applied directly upon noisy point clouds. Good methods for noise filtering are paramount for a robust stem isolation algorithm. There is a handful of methods to denoise stem point cloud data. One early technique is the Hough transformation, used by Simonse et al. (2003), Rabbani and van den Heuvel (2005) and Olofsson et al. (2014). The Hough transformation is a technique in which every data point “votes” for a cylinder center. Another method is to use saliency features of a point dataset to classify different parts of a TLS dataset of a forest
Corresponding author. E-mail address: [email protected] (T. de Conto).
http://dx.doi.org/10.1016/j.compag.2017.10.019 Received 21 January 2017; Received in revised form 17 October 2017; Accepted 22 October 2017 0168-1699/ © 2017 Elsevier B.V. All rights reserved.
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
described in detail by Liang et al. (2012) and Olofsson et al. (2014), for instance. Twenty-five scanned trees – 10 spruces, 9 pines and 6 eucalypts were manually calipered – at 1.3 m above ground (dbh) and every meter from the ground level up to 15 m above ground on the Norway spruce and Scots pine trees, and up to top height on the eucalypt trees. Those trees were processed using our implementations of the following methods. The callipered diameters were used as ground truth for validation, thus only trees with both laser data and field measurements were analyzed. On the present study, the effect of distance from the trees to the LiDAR sensor was not investigated – all studied trees were between 0 and 10 meters from the laser scanners on multi-scan setups.
plot into, for instance, stems and canopy. Several authors use this technique, such as Lalonde et al. (2006), Liang et al. (2012) and Xia et al. (2015). Another way to denoise stem data is to store the TLS data in a 3D voxel space and use morphological operations, as in the studies by Gorte and Pfeifer (2004), Gorte and Winterhalder (2004), Vonderach et al. (2012) and Raumonen et al. (2013). The most common way to model the volume or biomass of a tree’s stem is to fit circles or cylinders to TLS data points by non-linear, or linearized, least square adjustment (Simonse et al., 2003; Aschoff and Spiecker 2004; Hopkinson et al., 2004; Pfeifer and Winterhalder, 2004; Thies et al., 2004; Forsman and Halme 2005; Watt and Donoghue 2005, Henning and Radtke, 2006, Lalonde et al., 2006; Brolly and Király, 2009; Lindberg et al., 2012, Pueschel et al., 2013, Raumonen et al., 2013). Some studies use robust shape fitting of cylinders, like Liang et al. (2012, 2014b), who used iteratively reweighted least squares and m estimators. Another option is to use the RANSAC method (Chum 2005; Choi et al., 2009; Schnabel et al., 2007; Olofsson et al., 2014) which is a robust iterative method that can be used to find circles or cylinders directly in point cloud data. Many studies only evaluate one algorithm, often focusing on tree detection and estimating diameter at breast height (dbh). A few studies evaluated the performance of TLS on measurements of stem taper contrasted to field data, Henning and Radtke (2006), Maas et al. (2008), Liang et al. (2014b), and Mengesha et al. (2015). One feature of interest not yet investigated is to test various combinations of stem denoising and stem modelling algorithms to evaluate which combinations work best. The emphasis of this study is not to develop a tree detection algorithm that operates on TLS data, as other studies have already shown ways to do that. In this study we compare alternatives to extract stems from noisy data. Three different algorithms for stem denoising and three algorithms for stem modelling were chosen and all combinations of those were tested and evaluated. Also, the algorithms were tested on single tree point clouds with different densities, aiming to investigate the performance under lower density surveys. The code during the study is available as an open source package in R language, containing functions for dealing with TLS forest data.
2.2. Denoising techniques 2.2.1. Hough transformation (D1) The Hough transformation is a technique that searches for primitive shapes (e.g. lines, circles, etc.) on raster datasets (Illingworth and Kittler, 1987). On our study, we adapted the algorithm to find circular shapes on two dimensional horizontal layers of the point clouds. Taking a tree point cloud, it is first sliced from bottom to top – every 0.5 m. Every slice contains a stem segment assumed orthogonal to the ground – i.e. parallel to z = [0, 0, 1]. A baseline segment is defined at an arbitrary height interval (1–2 m above ground) and its approximate circle parameters are estimated – center coordinates (x, y) and radius, by applying the Hough transformation (Fig. 1). Points below the upper height limit of the baseline segment, which were far from its center (baseline’s radius plus 5 cm), were disregarded. Above the baseline segment, the search for circle centers in each segment is constrained to the x and y range found for the previous segment. For those segments, points further than its radius plus 1 cm from its estimated center are disregarded. Every pixel displays its local point density (ratio between its number of points and the pixel containing the most points). All pixels with more than 30% density are tested as center coordinates, and several radii between arbitrary limits (from 2.5 to 30 cm) are tested on each valid pixel. The radius between iterations increases by the pixel size – here we adopted 2.5 cm. Another, initially empty, raster layer with the same extent is used to count the number of circle intersections per pixel (votes). The pixel with most votes from circles with same radius is the estimated center – for that radius. The image displays 3 iterations (3 radii) of the Hough transformation. The same pixels are used as centers in all 3 iterations, the first iteration tests a small radius (left), the second iteration tests the true radius (center) and the third iteration tests a large radius (right). Note the central image contains the pixel with most votes, and thus best estimates the circle parameters.
2. Methods 2.1. Data acquisition The TLS data was acquired through multiple scan setups (Hilker et al., 2012). Mixed plots of Norway spruce (mean dbh = 29.57 cm, mean Ht = 22.4 m) and Scots pine (mean dbh = 25.04 cm, mean Ht = 20.7 m) of varying ages were scanned in Northern Sweden. Seven years old eucalypt trees were taken from a scanned plot in SouthEastern Brazil (mean dbh = 12.80 cm, mean Ht = 21.28 m) – dbh: diameter at breast height, Ht: total height. Table 1 summarizes the laser scanner settings of both surveys. All point clouds contained single trees only – 1.5 m radius from the tree’s approximate center. Tree detection at the plot level is not addressed in our study – such techniques were
2.2.2. Eigen decomposition of flat surfaces (D2) This method starts investigating the cloud pointwise. Every point has flatness (Eq. (1)) and verticality (angle between its normal vector and z = [0, 0, 1]) values assigned to it by performing eigen decomposition on its immediate point neighborhood of 30 points (Fig. 2). All (approximately) non-flat (FL < 0.8), non-vertical (angle with z > 10°) points are disregarded, then points which are close in the 3D space are clustered (up to 10 cm distant), and overlapping projections of those clusters on the xy plane are merged. This ensures that points belonging to the trunk are grouped together, even if gaps occur in the point cloud. Small clusters are assumed to be noise and also disregarded.
Table 1 Laser scanner settings of both surveys.
Tree species Scanning period Laser scanner Wave length (nm) Beam divergence (mrad) Pulse rate (kHz) Field of view (degrees)*
Sweden
Brazil
Pinus sylvestris L. Picea abies L. Karst. August 2013 Leica ScanStation C10 532 0.1 50 360h 270v
Eucalyptus sp.
FL = 1−
October 2015 RIEGL VZ-400 V-Line 1550 0.35 122 360h 100v
λ3 λ1 + λ2 + λ3
where FL = flatness; λi = eigenvalues, where λ1 ≥ λ2 ≥ λ3.
* h and v refer to the horizontal and vertical planes, respectively.
166
(1)
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Fig. 1. Hough transformation: a 3D cylinder cross-section is represented by a raster layer.
et al., 2007). One can consider it as a hybrid method between the RANSAC and least squares circle fit. Further parameter optimization was done by applying the Nelder-Mead algorithm (Nelder and Mead, 1965) on the best point samples, in order to achieve finer accuracy in the final estimations. The confidence level adopted at the RANSAC was 99%, with p = 0.8 – expected value of a denoised stem point cloud. The procedure is summarized in Fig. 5.
2.2.3. Voxel space neighborhoods (D3) On this method, at first a rough noise filter is applied twice to remove isolated points from the cloud – removing points with less than 2 neighbors within a radius of 5 and 10 cm, in two “sweeps”, respectively. A voxel space, comprised of 3 × 3 × 3 cm voxels, is assembled around the filtered point cloud. Immediate neighborhoods of 27 (33) voxels are merged (Fig. 3), and all points inside a voxel neighborhood become a cover set (Fig. 4). Cover sets are characterized by flatness (Eq. (1)) and verticality, likewise in method D2, and points belonging exclusively to (approximately) non-flat and non-vertical cover sets are disregarded. The remaining points are reorganized into larger voxel neighborhoods – 64 (43) voxels of 5 × 5 × 5 cm each. The largest cover set is used to estimate a global position for the trunk’s axis, assuming it is orthogonal to the ground. Cover sets which are too small or too far from the estimated global axis are also disregarded.
N=
log(1−P ) log(1−pn )
(2)
where N = number of iterations required to find at least one sample without outliers with a probability P; P = confidence level; n = number of observations per sample; p = probability that one point belongs to the model, which is estimated by M , with M being the number of inliers in the dataset and D D its total number of points.
2.3. Stem modelling 2.3.1. RANSAC circle fit (M1) The RANdom SAmple Consensus (RANSAC) algorithm considers the probability that all points in a random sample belong to a model (Eq. (2)) – i.e. the probability that a sample contains no outliers (Choi et al., 2009). The denoised point clouds were sliced every 0.5 m from ground to top, having their contained stem segments processed in sequence. On this method, we combined the RANSAC algorithm with least squares circle fit, sampling 15 points on every iteration, minimizing the distance from the sampled points to their fitted circle on each iteration and, finally picking the circle parameters that corresponded to the smallest distance. This circle fit step is performed independently for every tree segment. The sample size here adopted relied that the sample points suffice to estimate the optimal distance to the circle, without further calculations. For this reason, and also because this technique is applied after denoising, which leaves few remaining outliers, the sample size is higher than usual RANSAC implementations (Schnabel
2.3.2. IRTLS cylinder fit (M2) The denoised point cloud is segmented from bottom to top, every 0.5 m from the ground. On each segment a cylinder is modelled through Iterated Reweighted Total Least Squares (IRTLS). The first iteration of the ITRLS algorithm bases the initial parameters estimation on the segment’s gravity center. On every iteration, the sum of weighted squared distances from the points to the cylinder’s surface is minimized, employing the Nelder-Mead simplex algorithm (Nelder and Mead, 1965), applying weights calculated using Tukey’s biweight function (Eq. (3)) (Liang et al., 2012). The algorithm runs until convergence – i.e. when the sum of weighted squared distances stabilizes or a maximum number of iterations is reached. The distance from a point to a cylinder is defined in Eq. (4), and its vector components are shown in Fig. 6. Cylinder parameterization is described in detail by Lukács et al.
Fig. 2. The eigen decomposition extracts the eigenvalues and eigenvectors of a dataset. The figure shows a point’s immediate neighborhood on a cylinder point cloud. L1, L2 and L3 are the neighborhood’s eigenvectors, representing the directions from most to least variability, respectively. L3 is the direction of least variability, which is also normal to the cylinder’s axis.
167
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Fig. 3. Neighborhoods of 27 voxels viewed from 3 different angles. Voxels with same color belong to the same neighborhood. Notice that all neighborhoods have voxels in common – the central voxel of a neighborhood becomes an outer voxel of the next. Overlapping reduces voxel misclassification, since a single voxel may belong to flat and non-flat neighborhoods at the same time.
Fig. 4. Cover sets. (a) single tree point cloud segment (b) construction of voxel space (c) cover sets contained by neighborhoods of 27 voxels (d) cover sets distinguished by color. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
168
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Fig. 5. RANSAC circle fit on sample trunk segment. (a) sample tree with highlighted denoised stem (b) circles fitted on random samples of 15 points (c) inliers and outliers – according to the best fitted circle.
(1998). Fig. 7 shows the stepwise cylinder fitting procedure. 2 y 2 ⎤, c ⎦
()
⎧⎡1− w= ⎣ ⎨ ⎩ 0,
w = weight; d cyl y=
1.4826·MAD(dcyl)
if y > c
;
MAD = median absolute deviation; 1.4826 = scale factor that ensures MAD is asymptotically unbiased if dcyl ∼ N(0, σ2); c = efficiency constant at the normal distribution (5 corresponds to 95% efficiency).
if y ⩽ c (3)
where
Fig. 6. Point cloud and its model cylinder – n starts at the origin and ends at the orthogonal intersection with the cylinder’s axis a (equal to the q point in Eq. (4)). The optimal cylinder has the lowest sum of squared distances (dcyl) from the point cloud to its surface.
169
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Fig. 7. IRTLS cylinder fit of a trunk segment. (a) denoised tree point cloud highlighted (b) segment point cloud with undetected noise (c, d, e) fitted cylinder with noise highlighted – 20% lowest weights.
d cyl = |(p−q) × a|−r
(2014), Liang et al. (2012) and Raumonen et al. (2013). All the implemented algorithms are available on GitHub (github.com/tiagodc/ TreeLS).
(4)
where dcyl = distance from a point to the cylinder’s surface; p = point at the cylinder’s surface; q = point at the axis; a = axis’ direction; r = radius; |…| is the Euclidean norm and × is the vector cross product.
2.4. Performance evaluation Performance of all methods was assessed by comparing TLS estimated diameters with field measurements along the stems. Root mean squared error (Eq. (5)), bias (Eq. (6)), Pearson’s correlation and percentage of modelled stem segments were calculated as diagnostic measures.
2.3.3. RANSAC cylinder fit (M3) This method combines features from the previous two. It models cylinders on point cloud segments of 0.5 m by minimizing the sum of squared distances to a model cylinder (Eq. (4), Fig. 6) iteratively, through the RANSAC algorithm (Eq. (2), n = 20, p = 0.9, P = 0.99). Minimization is done also through the Nelder-Mead simplex algorithm (Nelder and Mead, 1965), taking as initial parameters estimations from the segment’s gravity center. For every segment, the chosen cylinder parameters were the ones which corresponded to the smallest distance between the modelled cylinder’s surface and its respective points. All studied algorithms were implemented in the R language (R Core Team 2016). The packages rLiDAR (Silva et al., 2015) and rgl (Adler and Murdoch, 2016) were used for reading LAS files and 3D plotting, respectively. The denoising and modelling techniques implemented on this study combined features of methods published by Olofsson et al.
1 n
RMSE =
bias =
1 n
n
∑i =1 (di−̂ di)2
(5)
n
∑i =1 (di−̂ di)
(6)
where RMSE = root mean squared error; d ̂ = estimated diameter at the ith observation; d = observed diameter at the ith observation; n = number of observations. Processing started after all individual tree point clouds had their point densities reduced. Three densities were tested, keeping up to 104, 170
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Fig. 8. Estimated vs. field measured diameters (D) for all combinations of denoising and modelling methods. The point clouds’ density was 104 points/m above ground (RMSE: root mean square error, r: Pearson’s correlation coefficient, D%: percentage of detected and modelled stem segments).
5 × 103 and 2 × 103 points/m above ground in each case. Points per meter above ground is a more relatable unit when dealing with single trees, since terrestrial laser scanning setups tend to collect more data at lower heights. The points were randomly sampled from their original point clouds, undergoing timed denoising and stem modelling procedures for all methods combinations.
is the most robust denoising method out of the three. The combination D1 × M1 showed best stem detection accuracy, only missing 11% of the total number of segments – compared to the total number of diameters calipered. Diameter size correlates negatively with tree height, thus smaller diameters are found at higher portions of the stems, further from the laser scanner and likely occluded by the canopy. Larger diameters were generally estimated more precisely and more often detected. Fig. 9 displays the performance of diameter estimation per portion of tree height. Diameters of stem segments below 50% of tree height were more accurately estimated. Above 50% height, both accuracy and stem detection decreased in all cases. Stems denoised by D1 were modelled up to taller heights, as it was the only method capable of detecting stem portions above 75% of tree height, and also presented more than 50% detection between 50% and 75% height – reinforcing its higher robustness. Best performance on stem segments above 50% height was achieved by combination D1 × M1, given the higher percentages of stem segments detection and lower bias (Fig. 9).
3. Results Best overall performance was obtained by denoising method D1 (Hough transformation) followed by modelling method M2 (IRTLS cylinder fit), as this combination produced the lowest RMSE and bias values, keeping a high stem segment detection rate (Fig. 8). Combination D3 (voxel space neighborhoods) × M3 (RANSAC cylinder fit) closely followed, as second best. Combinations denoised through methods D2 and D3 were generally less precise. Trees denoised using method D1 presented higher proportion of detected stem segments, also presenting more accurate modelling, on average. We can affirm, therefore, that D1 171
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Fig. 9. Residuals distributed by relative height intervals. D residuals = (estimated D – measured D), where D is diameter in meters. Percentage values shown beneath the boxplots are the proportion of detected stem segments in the respective height interval.
ground) were briefly tested, and the overall performance was very poor using our default parameters, as the algorithms often found false features that did not correspond to the trees in most cases. It is possible that even lower point densities can yield reliable results for some tree species, although optimal parameter setting must be investigated. For better assessment on similarities, Table 2 RMSE, bias, r and stem segments detection rate values went through a hierarchical cluster analysis (Fig. 10) (Kaufman and Rousseeuw, 1990). The Eucalypt sets were clustered together, isolated from spruce and pine, reflecting the less accurate estimations due to lower data quality. Most pine and spruce sets performed similarly, which is reflected on most them being aggregated in a single cluster. The point cloud density does not show any patterns in the clustering, which reinforces it has little effect on estimation quality.
Diagnostic measures were summarized and split in accordance to tree species, denoising method, stem modelling method and point cloud density, respectively (Table 2). Best performance was achieved on pine trees, followed by spruce and then eucalypt. During scanning of the eucalypt stand there was heavy wind effect, affecting the data quality. The same patterns as shown for all species combined in Fig. 8 can also be observed for the species separated. Point cloud density reduction had little influence on bias and RMSE values, with the advantage of speeding up the denoising methods. Percentage of stem segments detection tended to decrease slightly on reduced point cloud densities, mainly for spruce trees. Spruce trees have dense canopies starting from low heights, the laser pulses from the TLS system reach the stem less often naturally, thus further reducing the point cloud size makes stem points more sparse and harder to detect and model. The parameters for all methods were unchanged for the different point cloud densities, with little impact in performance. The parameters involving random sampling are unlikely to account for any differences in the results according to point density, but the ones related to point neighborhood size can yield some effect, and are therefore optimizable, since the area occupied by point neighborhoods of same size increases as point density decreases. Very low point densities (below 1000 points/m above
4. Discussion Liang et al. (2014b) compiled results from circle and cylinder fitting methods for stem isolation in TLS point clouds published between 2008 and 2012. They found that circle based techniques (RMSE: 1.48–5.69 cm) were generally more precise to diameter determination 172
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Table 2 Diagnostic measures for all methods combinations per tree species. Cloud density unit is points/m above ground. Denoising
Stem modelling
Eucalyptus sp. D1 M1
M2
M3
D2
M1
M2
M3
D3
M1
M2
M3
Pinus sylvestris D1 M1
M2
M3
D2
M1
M2
M3
D3
M1
M2
M3
Picea abies D1
M1
M2
M3
Cloud Density (×1000)
RMSE (cm)
*
Bias (cm)
*
Segments detected (%)
r
10 5 2 10 5 2 10 5 2
4.46 4.55 4.77 3.01 3.23 3.34 4.38 3.76 4.56
37 36 40 24 26 27 36 30 35
1.47 2.09 1.70 0.59 1.19 1.04 1.45 1.42 1.68
12 16 14 5 10 9 12 12 13
88 84 85 62 70 73 73 70 62
0.59 0.63 0.47 0.60 0.56 0.55 0.43 0.44 0.45
10 5 2 10 5 2 10 5 2
3.56 3.81 4.13 8.16 7.63 7.35 7.79 4.38 5.76
28 29 31 53 52 49 53 34 42
1.45 1.67 1.86 3.64 3.12 2.93 2.92 1.47 2.35
11 13 14 24 21 20 20 11 17
66 63 59 45 45 52 58 63 60
0.55 0.54 0.53 0.40 0.38 0.30 0.37 0.24 0.46
10 5 2 10 5 2 10 5 2
4.55 4.55 4.37 8.49 7.96 4.56 3.19 5.97 3.37
32 34 33 55 52 35 25 42 28
2.27 2.28 2.50 3.30 3.85 2.58 1.00 2.78 1.45
16 17 19 22 25 20 8 20 12
59 65 69 50 59 66 58 59 65
0.46 0.35 0.37 0.24 0.36 0.44 0.40 0.44 0.36
10 5 2 10 5 2 10 5 2
1.92 1.66 2.79 1.73 1.47 1.60 1.79 1.63 3.21
9 8 13 8 7 8 9 8 15
0.92 1.01 1.24 1.07 1.01 0.93 1.04 0.97 1.38
5 5 6 5 5 4 5 5 6
91 90 86 87 90 87 89 88 90
0.98 0.99 0.95 0.99 1.00 0.99 0.99 0.99 0.95
10 5 2 10 5 2 10 5 2
2.10 2.01 1.24 1.87 1.86 1.66 2.05 1.75 1.18
10 9 6 9 9 7 10 8 5
0.87 0.96 0.76 0.98 0.83 0.69 1.00 0.83 0.68
4 4 3 5 4 3 5 4 3
83 78 79 86 85 82 82 84 80
0.98 0.97 0.99 0.99 0.98 0.99 0.98 0.99 0.99
10 5 2 10 5 2 10 5 2
1.82 1.52 1.49 1.18 1.17 1.06 1.40 1.36 1.27
8 7 7 5 5 5 6 6 6
0.90 0.89 0.84 0.80 0.71 0.69 0.83 0.79 0.74
4 4 4 4 3 3 4 4 3
78 79 73 83 83 78 83 84 77
0.98 0.99 0.99 1.00 1.00 1.00 0.99 0.99 0.99
10 5 2 10 5 2 10 5 2
2.28 2.34 2.89 1.87 1.47 1.63 1.98 1.71 1.41
10 10 13 8 6 6 8 7 6
1.30 1.04 1.37 1.38 1.10 1.00 1.31 1.22 1.00
6 5 6 6 4 4 6 5 4
88 85 82 89 79 78 86 78 74
0.98 0.98 0.97 0.99 0.99 0.99 0.99 0.99 0.99
Relative RMSE (%)
Relative Bias (%)
(continued on next page)
173
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Table 2 (continued) Denoising
Stem modelling
Cloud Density (×1000)
RMSE (cm)
*
Bias (cm)
*
Segments detected (%)
r
D2
M1
10 5 2 10 5 2 10 5 2
3.06 3.58 2.82 1.88 2.47 2.51 2.96 1.63 2.02
12 14 11 8 10 9 12 7 8
1.50 1.66 1.46 1.28 1.50 1.54 1.50 1.10 1.34
6 7 6 5 6 6 6 4 5
72 68 64 81 71 66 82 71 73
0.95 0.93 0.96 0.99 0.97 0.97 0.96 0.99 0.98
10 5 2 10 5 2 10 5 2
2.17 2.16 3.05 1.70 1.84 1.26 2.22 1.67 1.19
9 8 12 7 7 5 9 7 5
1.30 1.32 1.10 1.31 1.19 0.98 1.48 0.96 0.73
5 5 4 5 5 4 6 4 3
65 64 58 71 66 58 77 64 58
0.98 0.98 0.95 0.99 0.99 1.00 0.98 0.99 0.99
M2
M3
D3
M1
M2
M3
Relative RMSE (%)
Relative Bias (%)
* Relative values were divided by the average diameter of all field measurements – for all segments. Fig. 10. Cluster analysis of diagnostic variables for all different denoising and modelling combinations, point cloud densities and tree species (c = cophenetic correlation).
denoising and modelling methods, close attention must be paid during the scanning to retrieve good quality data. In extreme cases the tree may appear deformed in the point cloud, yielding highly overestimated diameters on some of the automated methods. This effect can be seen in Fig. 8 (D2-M2, D2-M3 and D3-M2) as a cluster of points pertaining to a single Eucalypt tree. Liang et al. (2016) made a compilation of studies and results available for dealing with TLS data from forest inventories – mostly from temperate or boreal regions. Studies of TLS applications in the tropics are scarce, but with data becoming more easily available, the amount of studies is expected to rise in these areas shortly. The most robust denoising method was the Hough transformation (D1), as it preserved stem features further up in the canopies, also prompting better accuracy on the subsequent stem modelling phase, due to more well preserved stem point clouds. The best combination was denoising by Hough transformation (D1), followed by the IRTLS cylinder fit algorithm (M2) for stem modelling – although combinations involving M2 and M3 performed similarly, as seen in Table 2. The cylinder modelling approaches have better performance approximating stem features, as they can estimate the angles at which stem segments lean from a straight, vertical bole, while the circle based approach assumes perfect verticality of the stem, estimating only the segments’ radii on a horizontal plane. Combinations involving M2 and M3 performed similarly. Those slight differences may disappear with finer parameter tuning. For every
than cylinder based methods (RMSE: 3.7–9.17 cm). Olofsson et al. (2014) found diameter RMSE values in the range of 2.4–7.5 cm. They tested their adapted Hough transformation denoising and circle fit RANSAC algorithm on boreal tree species point clouds, obtaining best results for Scots pine (RMSE = 2.38 cm, bias = -0.6 cm), followed by Norway spruce (RMSE = 6.18 cm, bias = 1.54 cm) and deciduous species (RMSE = 7.48 cm, bias = 2.52 cm). Liang et al. (2014a, 2014b) also tested their algorithms (eigen decomposition and iterated reweighted least squares cylinder fit) on boreal species, finding a bias of 0.15 cm and RMSE values of 1.17 cm for Pine and 1.08 cm for Spruce. Furthermore, they also reported that automated estimations were as precise as visual assessment of the point cloud, and around 3 times less biased. The reported diagnostic measures were also obtained by contrasting their methods with manual field measurements of diameters along the stems. Most methods implemented on this study performed better on pine trees, followed by spruce and then eucalypt. Pine trees have less dense lower canopy, while spruce stems are often occluded by branches at low heights. Performance on eucalypt trees was poor in comparison, mainly due to wind effect during scanning given the sideways movement of the trees, which provoked the effect of “ghost stems” in the point cloud. In addition, the eucalypt stand was only seven years old, and the trees thin and slender, which explains the high relative values of RMSE and bias observed in Table 2. Our implemented algorithms were not able to properly model those trees, so in order to get good accuracy from stem
174
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
clusters. Testing nine different algorithm combinations we observed that estimations made on stem segments below 10 cm of diameter were generally much less accurate or not even detected. Without density reduction at higher portions of the canopy in the point cloud, high accuracy is likely achievable for even smaller diameters, decreasing this threshold. Since the most valuable, commercial wood is of larger diameters, concentrated on the first few meters above the ground, the modelling techniques have proven to be effective for measuring commercial wood on standing trees. From a biomass and carbon stock estimation perspective, more robust and refined stem extraction methods and reliable taper models are required to attain precise and accurate estimations of upper diameters from TLS point clouds, as well as increased level of detail on the tree modelling, including branches and possibly foliage. Some methods for full tree modelling are described by Dassot et al. (2012), Raumonen et al. (2013) and Hackenberg et al. (2015). Even though the data used here was obtained through terrestrial laser scanning, the algorithms tested can be applied on point clouds obtained using other technologies, such as mobile scanning and photogrammetric methods (Liang et al., 2014a). As the algorithms rely solely on geometrical features, requiring only 3D coordinates as input, any data source that provides this type of data can be targeted. The constant increase in popularity of LiDAR devices, which is associated to decreasing market prices, will make 3D point cloud data easily available and cheaper for a wider range of customers, therefore availability of robust, time efficient and flexible software for data processing is already a demand.
method implemented, both denoising and modelling, a set of general parameters was chosen based on empirical tests, on which the algorithms were capable to properly detect and subsequently model stems from varying tree architectures. There is no best global parameter combination for a given method and default parameters should work well on a wide variety of tree species and ages. Nevertheless, the parameter choice is likely optimizable for distinct target tree species/ architectures for all combinations, in order to improve robustness and accuracy for specific conditions, such as dense canopy or scattered branching. In this study, we opted for a default general parameter setting per method, regardless of tree species, in order to investigate the performance of fully automated procedures. Parameter tuning assessment under different forest conditions should be addressed in a simulation study. The RANSAC circle fit (M1) was the most stable stem modelling method, as its RMSE and bias values fluctuated the least. By taking a closer look at Table 2 it becomes clear that this happens due to the eucalypt trees. This method could process more robustly point clouds of worse quality, especially when applied in sequence from the Hough transformation denoising algorithm (D1). Diameter estimations above half of tree total height were less precise and accurate due to the longer distance from the laser sensor and occlusion by branches and the canopy layer. Larger laser footprint, as well as a sharper angle between laser beams and stem surfaces closer to the tree tops are also a possible source of error in the estimations. Some possible solutions to get further estimations of missing stem segments are regression, interpolation or species-specific taper equations based on the estimated values from the TLS data. Regarding some specifics of the denoising techniques, D1 is more robust but its performance is hindered in case of tilted trees and forking stems. The Hough transformation works well on slightly tilted trees, but in order to adapt the method to detect forked stems, extra parameters and conditions must be added to the algorithm. D2 and D3 overcome those issues by detecting plane-like features whose points can deviate from a strictly vertical distribution, regardless of stem forking. Those methods need denser point clouds at the canopy region though, as they are more prone to omit stem segments with scarce inliers on the stem surface. Regarding the modelling algorithms, M1 requires least amount of points, being less affected by absence of points at the stem surface. Its main issue relates to bent or tilted trees, as it will often overestimate the diameter for any stem segment whose section is not exactly horizontal. For generally straight trees this method is precise and robust, but for other cases, methods that extract a cylinder axis’ direction are better suited, such as M2 or M3. Those methods require connectivity among stem surface patches and a larger amount of points to yield good estimations of cylinder parameters, thus being more sensitive to point cloud density at the canopy layer. Fitting rings on a two-dimensional plane is more robust in this regard, since it is less demanding in terms of data amounts and patch connectivity, thus possibly a better choice for processing tree point clouds acquired by low cost or low resolution scanners. All methods were applied on reduced point density clouds, with no more than 3x105 points per tree, having little to no loss in terms of robustness and accuracy estimating diameters along the stems. Kelbe et al. (2015) estimated dbh from single scan, low density TLS point clouds, achieving a RMSE of 6.4 cm. According to Liang et al. (2014b), further improvement on TLS-based estimation of tree attributes could be achieved by using higher point density. Our study shows that such algorithms can provide accurate estimates on lower point densities, and rather than increasing point cloud density globally, adopting different point densities per height strata may suffice the goal of improved global precision/accuracy – e.g. by maintaining the original point density only at the canopy layer, avoiding missing stem segments. Fitting circles or cylinders on narrow features can be tricky, since the algorithms tend to falsely identify those shapes in small point
5. Conclusions All tested methods are accurate and precise on single tree point clouds of good quality. They perform generally better at lower portions of the stem. The Hough transformation (D1) followed by iteratively reweighted total least squares cylinder fit method (M2) presented best combined performance of denoising and stem modelling, respectively. D1 consists in a more robust denoising method as it is less sensitive to outliers, being able to capture higher percentages of stem features at upper tree heights. The cylinder based modelling methods (M2 and M3) have better modelling accuracy as they consider leaning angles in addition to radii estimations, thus being less prone to diameter overestimation. For roughly straight trees, a combination of D1 and M1 is recommended for robustness and stability, respectively. For varying tree architectures a cylinder modelling approach is recommended (M2, M3), and in case of forking stems, denoising methods D2 and D3 are better suited. All functions developed during this study are available as an open source R package on GitHub: github.com/tiagodc/TreeLS. Currently the package provides a collection of functions for manipulation of TLS point clouds, stem mapping at plot level, automatic stem isolation from single trees and taper extraction. Acknowledgments We would like to thank the Sven and Hildur Wingquist Foundation For Forest Research, for financing the Swedish data acquisition, managed by the forest remote sensing section at SLU Umeå. The Erasmus Mundus programme in Sustainable Forest and Nature Management, SUFONAMA. The Southern Swedish Forest Research Centre, SLU Alnarp. References Adler, D., Murdoch, D., 2016. rgl: 3D Visualization Using OpenGL. R Package Version 0. 96.0. < https://CRAN.R-project.org/package=rgl > . Aschoff, T., Spiecker, H., 2004. Algorithms for the automatic detection of trees in laser scanner data. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 36, 71–74.
175
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Lukács, G., Martin, R., Marshall, D., 1998. Faithful least-squares fitting of spheres, cylinders, cones and tori forreliable segmentation. Proc. Eur. Conf. Comput. Vis. 1, 671–686. Maas, H., Bienert, A., Scheller, S., Keane, E., 2008. Automatic forest inventory parameter determination from terrestrial laser scanner data. Int. J. Remote Sens. 29, 1579–1593. Mengesha, T., Hawkins, M., Nieuwenhuis, M., 2015. Validation of terrestrial laser scanning data using conventional forest inventory methods. Eur. J. For. Res. 134, 211–222. Nelder, J.A., Mead, R., 1965. A simplex method for function minimization. Comput. J. 7 (4), 308–313. Olofsson, K., Holmgren, J., Olsson, H., 2014. Tree stem and height measurements using terrestrial laser scanning and the RANSAC algorithm. Remote Sens. 6 (5), 4323–4344. Pfeifer, N., Winterhalder, D., 2004. Modelling of tree cross sections from terrestrial laser scanning data with free-form curves. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 36, 76–81. Pueschel, P., Newnham, G., Rock, G., Udelhoven, T., Werner, W., Hill, J., 2013. The influence of scan mode and circle fitting on tree stem detection, stem diameter and volume extraction from terrestrial laser scans. ISPRS J. Photogramm. Remote Sens. 77, 44–56. R Core Team, 2016. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
Brolly, G., Király, G., 2009. Algorithms for stem mapping by means of terrestrial laser scanning. Acta Silv. Lignaria Hung. 5, 119–130. Choi, S., Kim, T., Yu, W., 2009. Performance evaluation of RANSAC family. Br. Mach. Vis. Assoc. 81. Chum, O., 2005. Two-view geometry estimation by random sample and consensus. Ph.D. Thesis. Czech Technical University, Prague, Czech. Dassot, M., Colin, A., Santenoise, P., Fournier, M., Constant, T., 2012. Terrestrial laser scanning for measuring the solid wood volume, including branches, of adult standing trees in the forest environment. Comput. Electron. Agric. 89, 86–93. Forsman, P., Halme, A., 2005. 3-D mapping of natural environments with trees by means of mobile perception. IEEE Trans. Robot. 21, 482–490. Gorte, B., Pfeifer, N., 2004. Structuring laser-scanned trees using 3D mathematical morphology. Int. Arch. Photogramm. Remote Sens. 35, 929–933. Gorte, B., Winterhalder, D., 2004. Reconstruction of laser-scanned trees using filter operations in the 3D raster domain. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 36, 39–44. Hackenberg, J., Morhart, C., Sheppard, J., Spiecker, H., Disney, M., 2014. Highly accurate tree models derived from terrestrial laser scan data: a method description. Forests 5 (5), 1069–1105. Hackenberg, J., Morhart, C., Sheppard, J., Spiecker, H., Disney, M., 2015. SimpleTree – an efficient open source tool to build tree models from TLS clouds. Forests 6 (12), 4245–4294. Henning, J.G., Radtke, P.J., 2006. Detailed stem measurements of standing trees from ground-based scanning LiDAR. For. Sci. 52, 67–80. Hilker, T., Coops, N.C., Culvenor, D.S., Newnham, G., Wulders, M.A., Bater, C.W., Siggins, A., 2012. A simple technique for co-registration of terrestrial LiDAR observations for forestry applications. Remote Sens. Lett. 3 (3), 239–247. Hopkinson, C., Chasmer, L., Young-Pow, C., Treitz, P., 2004. Assessing forest metrics with a ground-based scanning LiDAR. Can. J. For. Res. 34, 573–583. Illingworth, J., Kittler, J., 1987. The adaptive Hough transform. IEEE Trans. Pattern Anal. Mach. Intell. PAMI-9 (5), 690–698. Kaufman, L., Rousseeuw, P., 1990. Finding Groups in Data: An Introduction to Cluster Analysis. John Wiley & Sons Inc, New York. Kelbe, D., van Aardt, J., Romanczyk, P., van Leeuwen, M., Cawse-Nicholson, K., 2015. Single-scan stem reconstruction using low-resolution terrestrial laser scanner data. IEEE J. Sel. Top. Appl. Earth Observ. Remote Sens. 8 (7), 3414–3427. Lalonde, J.F., Vandapel, N., Huber, D., Hebert, M., 2006. Natural terrain classification using three-dimensional ladar data for ground robot mobility. J. Field Robot. 23, 839–861. Liang, X., Litkey, P., Hyyppä, J., Kaartinen, H., Vastaranta, M., Holopainen, M., 2012. Automatic stem mapping using single-scan terrestrial laser scanning. IEEE Trans. Geosci. Remote Sens. 50 (2), 661–670. Liang, X., Jaakkola, A., Wang, Y., Hyyppä, J., Honkavaara, E., Liu, J., Kaartinen, H., 2014a. The use of a hand-held camera for individual tree 3D mapping in forest sample plots. Remote Sens. 6, 6587–6603. Liang, X., Kankare, V., Yu, X., Hyyppä, J., Holopainen, M., 2014b. Automated stem curve measurement using terrestrial laser scanning. IEEE Trans. Geosci. Remote Sens. 52 (3), 1739–1748. Liang, X., Kankare, V., Hyyppä, J., Wang, Y., Kukko, A., Holopainen, M., Vastaranta, M., 2016. Terrestrial laser scanning in forest inventories. ISPRS J. Photogramm. Remote Sens. 115, 63–77. Lindberg, E., Holmgren, J., Olofsson, K., Olsson, H., 2012. Estimation of stem attributes using a combination of terrestrial and airborne laser scanning. Eur. J. For. Res. 131, 1917–1931.
176
Contents lists available at ScienceDirect
Computers and Electronics in Agriculture journal homepage: www.elsevier.com/locate/compag
Original papers
Performance of stem denoising and stem modelling algorithms on single tree point clouds from terrestrial laser scanning
MARK
⁎
Tiago de Contoa, , Kenneth Olofssonb, Eric Bastos Görgensc, Luiz Carlos Estraviz Rodriguezd, Gustavo Almeidad a
Swedish University of Agricultural Sciences (SLU), Sundsvägen 6, 230 53 Alnarp, Sweden Swedish University of Agricultural Sciences (SLU), Skogsmarksgränd, 907 36 Umeå, Sweden c Universidade Federal dos Vales do Jequitinhonha e Mucuri (UFVJM), Campus JK, Alto da Jacuba, 39100-000 Diamantina, MG, Brazil d Escola Superior de Agricultura Luiz de Queiroz – Universidade de São Paulo (ESALQ-USP), Avenida Pádua Dias, 11, 13418-900 Piracicaba, SP, Brazil b
A R T I C L E I N F O
A B S T R A C T
Keywords: LiDAR Cylinder/circle fit Tropical and boreal tree species Robust estimation RANSAC
The present study assessed the performance of three different methods of stem denoising and three different methods of stem modelling on terrestrial laser scanner (TLS) point clouds containing single trees – thus validating all tested methods, which were made available as an open source software package in the R language. The methods were adapted from common TLS stem detection techniques and rely on finding one main trunk in a point cloud by denoising the data to precisely extract only stem points, followed by a circle or cylinder fitting procedure on stem segments. The combination of the Hough transformation stem denoising method and the iteratively reweighted total least squares modelling method had best overall performance – achieving 2.15 cm of RMSE and 1.09 cm of bias when estimating diameters along the stems, detecting 80% of all stem segments measured on field surveys. All algorithms performed better on point clouds of boreal species, in comparison to tropical Eucalypt. The point clouds underwent reduction of point density, which increased processing speed on the stem denoising algorithms, with little effect on diameter estimation quality.
1. Introduction Terrestrial Laser Scanning (TLS) sensors are increasing in popularity in the Forestry sector, due to the high resolution and precision of the three-dimensional information provided at the sample plot level, making further 3D modelling and tree reconstruction possible and allowing foresters to extract dendrometric variables with high accuracy from point clouds (Liang et al., 2016). In order to apply the TLS technology in forest inventory at its full potential, robust tool sets to extract useful information from TLS point clouds, and that perform well on a variety of circumstances are required – especially stem denoising and stem modelling, since stem volume is the most targeted variable in commercial forest management. TLS is a technique that could be used to develop locally adjusted stem taper functions (Trincado and Burkhart, 2006; Sun et al., 2016). Routines for processing TLS point clouds containing tree data have been tested, with applications ranging from isolation of single trees in scanned plots (Liang et al., 2012; Olofsson et al., 2014) to segmentation of tree components (e.g. branches, stem and canopy) (Raumonen et al., 2013; Hackenberg et al., 2014).
⁎
In order to automate the process of isolating tree components in point clouds, general assumptions need to be made about a tree’s geometrical structure. A reliable algorithm for stem isolation, thus, requires a level of generalization that allows it to capture the essential information of a variety of tree shapes in a point cloud. Some commonly adopted assumptions are the vertical orientation of the bole (approximately orthogonal angle with the ground), approximately circular stem cross sections and approximately cylindrical stem segments. Those assumptions are the core of many already implemented algorithms, but a key issue must be tackled beforehand: fitting circles and cylinders works well on point clouds with the stem clearly visible, but can’t be applied directly upon noisy point clouds. Good methods for noise filtering are paramount for a robust stem isolation algorithm. There is a handful of methods to denoise stem point cloud data. One early technique is the Hough transformation, used by Simonse et al. (2003), Rabbani and van den Heuvel (2005) and Olofsson et al. (2014). The Hough transformation is a technique in which every data point “votes” for a cylinder center. Another method is to use saliency features of a point dataset to classify different parts of a TLS dataset of a forest
Corresponding author. E-mail address: [email protected] (T. de Conto).
http://dx.doi.org/10.1016/j.compag.2017.10.019 Received 21 January 2017; Received in revised form 17 October 2017; Accepted 22 October 2017 0168-1699/ © 2017 Elsevier B.V. All rights reserved.
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
described in detail by Liang et al. (2012) and Olofsson et al. (2014), for instance. Twenty-five scanned trees – 10 spruces, 9 pines and 6 eucalypts were manually calipered – at 1.3 m above ground (dbh) and every meter from the ground level up to 15 m above ground on the Norway spruce and Scots pine trees, and up to top height on the eucalypt trees. Those trees were processed using our implementations of the following methods. The callipered diameters were used as ground truth for validation, thus only trees with both laser data and field measurements were analyzed. On the present study, the effect of distance from the trees to the LiDAR sensor was not investigated – all studied trees were between 0 and 10 meters from the laser scanners on multi-scan setups.
plot into, for instance, stems and canopy. Several authors use this technique, such as Lalonde et al. (2006), Liang et al. (2012) and Xia et al. (2015). Another way to denoise stem data is to store the TLS data in a 3D voxel space and use morphological operations, as in the studies by Gorte and Pfeifer (2004), Gorte and Winterhalder (2004), Vonderach et al. (2012) and Raumonen et al. (2013). The most common way to model the volume or biomass of a tree’s stem is to fit circles or cylinders to TLS data points by non-linear, or linearized, least square adjustment (Simonse et al., 2003; Aschoff and Spiecker 2004; Hopkinson et al., 2004; Pfeifer and Winterhalder, 2004; Thies et al., 2004; Forsman and Halme 2005; Watt and Donoghue 2005, Henning and Radtke, 2006, Lalonde et al., 2006; Brolly and Király, 2009; Lindberg et al., 2012, Pueschel et al., 2013, Raumonen et al., 2013). Some studies use robust shape fitting of cylinders, like Liang et al. (2012, 2014b), who used iteratively reweighted least squares and m estimators. Another option is to use the RANSAC method (Chum 2005; Choi et al., 2009; Schnabel et al., 2007; Olofsson et al., 2014) which is a robust iterative method that can be used to find circles or cylinders directly in point cloud data. Many studies only evaluate one algorithm, often focusing on tree detection and estimating diameter at breast height (dbh). A few studies evaluated the performance of TLS on measurements of stem taper contrasted to field data, Henning and Radtke (2006), Maas et al. (2008), Liang et al. (2014b), and Mengesha et al. (2015). One feature of interest not yet investigated is to test various combinations of stem denoising and stem modelling algorithms to evaluate which combinations work best. The emphasis of this study is not to develop a tree detection algorithm that operates on TLS data, as other studies have already shown ways to do that. In this study we compare alternatives to extract stems from noisy data. Three different algorithms for stem denoising and three algorithms for stem modelling were chosen and all combinations of those were tested and evaluated. Also, the algorithms were tested on single tree point clouds with different densities, aiming to investigate the performance under lower density surveys. The code during the study is available as an open source package in R language, containing functions for dealing with TLS forest data.
2.2. Denoising techniques 2.2.1. Hough transformation (D1) The Hough transformation is a technique that searches for primitive shapes (e.g. lines, circles, etc.) on raster datasets (Illingworth and Kittler, 1987). On our study, we adapted the algorithm to find circular shapes on two dimensional horizontal layers of the point clouds. Taking a tree point cloud, it is first sliced from bottom to top – every 0.5 m. Every slice contains a stem segment assumed orthogonal to the ground – i.e. parallel to z = [0, 0, 1]. A baseline segment is defined at an arbitrary height interval (1–2 m above ground) and its approximate circle parameters are estimated – center coordinates (x, y) and radius, by applying the Hough transformation (Fig. 1). Points below the upper height limit of the baseline segment, which were far from its center (baseline’s radius plus 5 cm), were disregarded. Above the baseline segment, the search for circle centers in each segment is constrained to the x and y range found for the previous segment. For those segments, points further than its radius plus 1 cm from its estimated center are disregarded. Every pixel displays its local point density (ratio between its number of points and the pixel containing the most points). All pixels with more than 30% density are tested as center coordinates, and several radii between arbitrary limits (from 2.5 to 30 cm) are tested on each valid pixel. The radius between iterations increases by the pixel size – here we adopted 2.5 cm. Another, initially empty, raster layer with the same extent is used to count the number of circle intersections per pixel (votes). The pixel with most votes from circles with same radius is the estimated center – for that radius. The image displays 3 iterations (3 radii) of the Hough transformation. The same pixels are used as centers in all 3 iterations, the first iteration tests a small radius (left), the second iteration tests the true radius (center) and the third iteration tests a large radius (right). Note the central image contains the pixel with most votes, and thus best estimates the circle parameters.
2. Methods 2.1. Data acquisition The TLS data was acquired through multiple scan setups (Hilker et al., 2012). Mixed plots of Norway spruce (mean dbh = 29.57 cm, mean Ht = 22.4 m) and Scots pine (mean dbh = 25.04 cm, mean Ht = 20.7 m) of varying ages were scanned in Northern Sweden. Seven years old eucalypt trees were taken from a scanned plot in SouthEastern Brazil (mean dbh = 12.80 cm, mean Ht = 21.28 m) – dbh: diameter at breast height, Ht: total height. Table 1 summarizes the laser scanner settings of both surveys. All point clouds contained single trees only – 1.5 m radius from the tree’s approximate center. Tree detection at the plot level is not addressed in our study – such techniques were
2.2.2. Eigen decomposition of flat surfaces (D2) This method starts investigating the cloud pointwise. Every point has flatness (Eq. (1)) and verticality (angle between its normal vector and z = [0, 0, 1]) values assigned to it by performing eigen decomposition on its immediate point neighborhood of 30 points (Fig. 2). All (approximately) non-flat (FL < 0.8), non-vertical (angle with z > 10°) points are disregarded, then points which are close in the 3D space are clustered (up to 10 cm distant), and overlapping projections of those clusters on the xy plane are merged. This ensures that points belonging to the trunk are grouped together, even if gaps occur in the point cloud. Small clusters are assumed to be noise and also disregarded.
Table 1 Laser scanner settings of both surveys.
Tree species Scanning period Laser scanner Wave length (nm) Beam divergence (mrad) Pulse rate (kHz) Field of view (degrees)*
Sweden
Brazil
Pinus sylvestris L. Picea abies L. Karst. August 2013 Leica ScanStation C10 532 0.1 50 360h 270v
Eucalyptus sp.
FL = 1−
October 2015 RIEGL VZ-400 V-Line 1550 0.35 122 360h 100v
λ3 λ1 + λ2 + λ3
where FL = flatness; λi = eigenvalues, where λ1 ≥ λ2 ≥ λ3.
* h and v refer to the horizontal and vertical planes, respectively.
166
(1)
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Fig. 1. Hough transformation: a 3D cylinder cross-section is represented by a raster layer.
et al., 2007). One can consider it as a hybrid method between the RANSAC and least squares circle fit. Further parameter optimization was done by applying the Nelder-Mead algorithm (Nelder and Mead, 1965) on the best point samples, in order to achieve finer accuracy in the final estimations. The confidence level adopted at the RANSAC was 99%, with p = 0.8 – expected value of a denoised stem point cloud. The procedure is summarized in Fig. 5.
2.2.3. Voxel space neighborhoods (D3) On this method, at first a rough noise filter is applied twice to remove isolated points from the cloud – removing points with less than 2 neighbors within a radius of 5 and 10 cm, in two “sweeps”, respectively. A voxel space, comprised of 3 × 3 × 3 cm voxels, is assembled around the filtered point cloud. Immediate neighborhoods of 27 (33) voxels are merged (Fig. 3), and all points inside a voxel neighborhood become a cover set (Fig. 4). Cover sets are characterized by flatness (Eq. (1)) and verticality, likewise in method D2, and points belonging exclusively to (approximately) non-flat and non-vertical cover sets are disregarded. The remaining points are reorganized into larger voxel neighborhoods – 64 (43) voxels of 5 × 5 × 5 cm each. The largest cover set is used to estimate a global position for the trunk’s axis, assuming it is orthogonal to the ground. Cover sets which are too small or too far from the estimated global axis are also disregarded.
N=
log(1−P ) log(1−pn )
(2)
where N = number of iterations required to find at least one sample without outliers with a probability P; P = confidence level; n = number of observations per sample; p = probability that one point belongs to the model, which is estimated by M , with M being the number of inliers in the dataset and D D its total number of points.
2.3. Stem modelling 2.3.1. RANSAC circle fit (M1) The RANdom SAmple Consensus (RANSAC) algorithm considers the probability that all points in a random sample belong to a model (Eq. (2)) – i.e. the probability that a sample contains no outliers (Choi et al., 2009). The denoised point clouds were sliced every 0.5 m from ground to top, having their contained stem segments processed in sequence. On this method, we combined the RANSAC algorithm with least squares circle fit, sampling 15 points on every iteration, minimizing the distance from the sampled points to their fitted circle on each iteration and, finally picking the circle parameters that corresponded to the smallest distance. This circle fit step is performed independently for every tree segment. The sample size here adopted relied that the sample points suffice to estimate the optimal distance to the circle, without further calculations. For this reason, and also because this technique is applied after denoising, which leaves few remaining outliers, the sample size is higher than usual RANSAC implementations (Schnabel
2.3.2. IRTLS cylinder fit (M2) The denoised point cloud is segmented from bottom to top, every 0.5 m from the ground. On each segment a cylinder is modelled through Iterated Reweighted Total Least Squares (IRTLS). The first iteration of the ITRLS algorithm bases the initial parameters estimation on the segment’s gravity center. On every iteration, the sum of weighted squared distances from the points to the cylinder’s surface is minimized, employing the Nelder-Mead simplex algorithm (Nelder and Mead, 1965), applying weights calculated using Tukey’s biweight function (Eq. (3)) (Liang et al., 2012). The algorithm runs until convergence – i.e. when the sum of weighted squared distances stabilizes or a maximum number of iterations is reached. The distance from a point to a cylinder is defined in Eq. (4), and its vector components are shown in Fig. 6. Cylinder parameterization is described in detail by Lukács et al.
Fig. 2. The eigen decomposition extracts the eigenvalues and eigenvectors of a dataset. The figure shows a point’s immediate neighborhood on a cylinder point cloud. L1, L2 and L3 are the neighborhood’s eigenvectors, representing the directions from most to least variability, respectively. L3 is the direction of least variability, which is also normal to the cylinder’s axis.
167
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Fig. 3. Neighborhoods of 27 voxels viewed from 3 different angles. Voxels with same color belong to the same neighborhood. Notice that all neighborhoods have voxels in common – the central voxel of a neighborhood becomes an outer voxel of the next. Overlapping reduces voxel misclassification, since a single voxel may belong to flat and non-flat neighborhoods at the same time.
Fig. 4. Cover sets. (a) single tree point cloud segment (b) construction of voxel space (c) cover sets contained by neighborhoods of 27 voxels (d) cover sets distinguished by color. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
168
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Fig. 5. RANSAC circle fit on sample trunk segment. (a) sample tree with highlighted denoised stem (b) circles fitted on random samples of 15 points (c) inliers and outliers – according to the best fitted circle.
(1998). Fig. 7 shows the stepwise cylinder fitting procedure. 2 y 2 ⎤, c ⎦
()
⎧⎡1− w= ⎣ ⎨ ⎩ 0,
w = weight; d cyl y=
1.4826·MAD(dcyl)
if y > c
;
MAD = median absolute deviation; 1.4826 = scale factor that ensures MAD is asymptotically unbiased if dcyl ∼ N(0, σ2); c = efficiency constant at the normal distribution (5 corresponds to 95% efficiency).
if y ⩽ c (3)
where
Fig. 6. Point cloud and its model cylinder – n starts at the origin and ends at the orthogonal intersection with the cylinder’s axis a (equal to the q point in Eq. (4)). The optimal cylinder has the lowest sum of squared distances (dcyl) from the point cloud to its surface.
169
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Fig. 7. IRTLS cylinder fit of a trunk segment. (a) denoised tree point cloud highlighted (b) segment point cloud with undetected noise (c, d, e) fitted cylinder with noise highlighted – 20% lowest weights.
d cyl = |(p−q) × a|−r
(2014), Liang et al. (2012) and Raumonen et al. (2013). All the implemented algorithms are available on GitHub (github.com/tiagodc/ TreeLS).
(4)
where dcyl = distance from a point to the cylinder’s surface; p = point at the cylinder’s surface; q = point at the axis; a = axis’ direction; r = radius; |…| is the Euclidean norm and × is the vector cross product.
2.4. Performance evaluation Performance of all methods was assessed by comparing TLS estimated diameters with field measurements along the stems. Root mean squared error (Eq. (5)), bias (Eq. (6)), Pearson’s correlation and percentage of modelled stem segments were calculated as diagnostic measures.
2.3.3. RANSAC cylinder fit (M3) This method combines features from the previous two. It models cylinders on point cloud segments of 0.5 m by minimizing the sum of squared distances to a model cylinder (Eq. (4), Fig. 6) iteratively, through the RANSAC algorithm (Eq. (2), n = 20, p = 0.9, P = 0.99). Minimization is done also through the Nelder-Mead simplex algorithm (Nelder and Mead, 1965), taking as initial parameters estimations from the segment’s gravity center. For every segment, the chosen cylinder parameters were the ones which corresponded to the smallest distance between the modelled cylinder’s surface and its respective points. All studied algorithms were implemented in the R language (R Core Team 2016). The packages rLiDAR (Silva et al., 2015) and rgl (Adler and Murdoch, 2016) were used for reading LAS files and 3D plotting, respectively. The denoising and modelling techniques implemented on this study combined features of methods published by Olofsson et al.
1 n
RMSE =
bias =
1 n
n
∑i =1 (di−̂ di)2
(5)
n
∑i =1 (di−̂ di)
(6)
where RMSE = root mean squared error; d ̂ = estimated diameter at the ith observation; d = observed diameter at the ith observation; n = number of observations. Processing started after all individual tree point clouds had their point densities reduced. Three densities were tested, keeping up to 104, 170
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Fig. 8. Estimated vs. field measured diameters (D) for all combinations of denoising and modelling methods. The point clouds’ density was 104 points/m above ground (RMSE: root mean square error, r: Pearson’s correlation coefficient, D%: percentage of detected and modelled stem segments).
5 × 103 and 2 × 103 points/m above ground in each case. Points per meter above ground is a more relatable unit when dealing with single trees, since terrestrial laser scanning setups tend to collect more data at lower heights. The points were randomly sampled from their original point clouds, undergoing timed denoising and stem modelling procedures for all methods combinations.
is the most robust denoising method out of the three. The combination D1 × M1 showed best stem detection accuracy, only missing 11% of the total number of segments – compared to the total number of diameters calipered. Diameter size correlates negatively with tree height, thus smaller diameters are found at higher portions of the stems, further from the laser scanner and likely occluded by the canopy. Larger diameters were generally estimated more precisely and more often detected. Fig. 9 displays the performance of diameter estimation per portion of tree height. Diameters of stem segments below 50% of tree height were more accurately estimated. Above 50% height, both accuracy and stem detection decreased in all cases. Stems denoised by D1 were modelled up to taller heights, as it was the only method capable of detecting stem portions above 75% of tree height, and also presented more than 50% detection between 50% and 75% height – reinforcing its higher robustness. Best performance on stem segments above 50% height was achieved by combination D1 × M1, given the higher percentages of stem segments detection and lower bias (Fig. 9).
3. Results Best overall performance was obtained by denoising method D1 (Hough transformation) followed by modelling method M2 (IRTLS cylinder fit), as this combination produced the lowest RMSE and bias values, keeping a high stem segment detection rate (Fig. 8). Combination D3 (voxel space neighborhoods) × M3 (RANSAC cylinder fit) closely followed, as second best. Combinations denoised through methods D2 and D3 were generally less precise. Trees denoised using method D1 presented higher proportion of detected stem segments, also presenting more accurate modelling, on average. We can affirm, therefore, that D1 171
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Fig. 9. Residuals distributed by relative height intervals. D residuals = (estimated D – measured D), where D is diameter in meters. Percentage values shown beneath the boxplots are the proportion of detected stem segments in the respective height interval.
ground) were briefly tested, and the overall performance was very poor using our default parameters, as the algorithms often found false features that did not correspond to the trees in most cases. It is possible that even lower point densities can yield reliable results for some tree species, although optimal parameter setting must be investigated. For better assessment on similarities, Table 2 RMSE, bias, r and stem segments detection rate values went through a hierarchical cluster analysis (Fig. 10) (Kaufman and Rousseeuw, 1990). The Eucalypt sets were clustered together, isolated from spruce and pine, reflecting the less accurate estimations due to lower data quality. Most pine and spruce sets performed similarly, which is reflected on most them being aggregated in a single cluster. The point cloud density does not show any patterns in the clustering, which reinforces it has little effect on estimation quality.
Diagnostic measures were summarized and split in accordance to tree species, denoising method, stem modelling method and point cloud density, respectively (Table 2). Best performance was achieved on pine trees, followed by spruce and then eucalypt. During scanning of the eucalypt stand there was heavy wind effect, affecting the data quality. The same patterns as shown for all species combined in Fig. 8 can also be observed for the species separated. Point cloud density reduction had little influence on bias and RMSE values, with the advantage of speeding up the denoising methods. Percentage of stem segments detection tended to decrease slightly on reduced point cloud densities, mainly for spruce trees. Spruce trees have dense canopies starting from low heights, the laser pulses from the TLS system reach the stem less often naturally, thus further reducing the point cloud size makes stem points more sparse and harder to detect and model. The parameters for all methods were unchanged for the different point cloud densities, with little impact in performance. The parameters involving random sampling are unlikely to account for any differences in the results according to point density, but the ones related to point neighborhood size can yield some effect, and are therefore optimizable, since the area occupied by point neighborhoods of same size increases as point density decreases. Very low point densities (below 1000 points/m above
4. Discussion Liang et al. (2014b) compiled results from circle and cylinder fitting methods for stem isolation in TLS point clouds published between 2008 and 2012. They found that circle based techniques (RMSE: 1.48–5.69 cm) were generally more precise to diameter determination 172
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Table 2 Diagnostic measures for all methods combinations per tree species. Cloud density unit is points/m above ground. Denoising
Stem modelling
Eucalyptus sp. D1 M1
M2
M3
D2
M1
M2
M3
D3
M1
M2
M3
Pinus sylvestris D1 M1
M2
M3
D2
M1
M2
M3
D3
M1
M2
M3
Picea abies D1
M1
M2
M3
Cloud Density (×1000)
RMSE (cm)
*
Bias (cm)
*
Segments detected (%)
r
10 5 2 10 5 2 10 5 2
4.46 4.55 4.77 3.01 3.23 3.34 4.38 3.76 4.56
37 36 40 24 26 27 36 30 35
1.47 2.09 1.70 0.59 1.19 1.04 1.45 1.42 1.68
12 16 14 5 10 9 12 12 13
88 84 85 62 70 73 73 70 62
0.59 0.63 0.47 0.60 0.56 0.55 0.43 0.44 0.45
10 5 2 10 5 2 10 5 2
3.56 3.81 4.13 8.16 7.63 7.35 7.79 4.38 5.76
28 29 31 53 52 49 53 34 42
1.45 1.67 1.86 3.64 3.12 2.93 2.92 1.47 2.35
11 13 14 24 21 20 20 11 17
66 63 59 45 45 52 58 63 60
0.55 0.54 0.53 0.40 0.38 0.30 0.37 0.24 0.46
10 5 2 10 5 2 10 5 2
4.55 4.55 4.37 8.49 7.96 4.56 3.19 5.97 3.37
32 34 33 55 52 35 25 42 28
2.27 2.28 2.50 3.30 3.85 2.58 1.00 2.78 1.45
16 17 19 22 25 20 8 20 12
59 65 69 50 59 66 58 59 65
0.46 0.35 0.37 0.24 0.36 0.44 0.40 0.44 0.36
10 5 2 10 5 2 10 5 2
1.92 1.66 2.79 1.73 1.47 1.60 1.79 1.63 3.21
9 8 13 8 7 8 9 8 15
0.92 1.01 1.24 1.07 1.01 0.93 1.04 0.97 1.38
5 5 6 5 5 4 5 5 6
91 90 86 87 90 87 89 88 90
0.98 0.99 0.95 0.99 1.00 0.99 0.99 0.99 0.95
10 5 2 10 5 2 10 5 2
2.10 2.01 1.24 1.87 1.86 1.66 2.05 1.75 1.18
10 9 6 9 9 7 10 8 5
0.87 0.96 0.76 0.98 0.83 0.69 1.00 0.83 0.68
4 4 3 5 4 3 5 4 3
83 78 79 86 85 82 82 84 80
0.98 0.97 0.99 0.99 0.98 0.99 0.98 0.99 0.99
10 5 2 10 5 2 10 5 2
1.82 1.52 1.49 1.18 1.17 1.06 1.40 1.36 1.27
8 7 7 5 5 5 6 6 6
0.90 0.89 0.84 0.80 0.71 0.69 0.83 0.79 0.74
4 4 4 4 3 3 4 4 3
78 79 73 83 83 78 83 84 77
0.98 0.99 0.99 1.00 1.00 1.00 0.99 0.99 0.99
10 5 2 10 5 2 10 5 2
2.28 2.34 2.89 1.87 1.47 1.63 1.98 1.71 1.41
10 10 13 8 6 6 8 7 6
1.30 1.04 1.37 1.38 1.10 1.00 1.31 1.22 1.00
6 5 6 6 4 4 6 5 4
88 85 82 89 79 78 86 78 74
0.98 0.98 0.97 0.99 0.99 0.99 0.99 0.99 0.99
Relative RMSE (%)
Relative Bias (%)
(continued on next page)
173
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Table 2 (continued) Denoising
Stem modelling
Cloud Density (×1000)
RMSE (cm)
*
Bias (cm)
*
Segments detected (%)
r
D2
M1
10 5 2 10 5 2 10 5 2
3.06 3.58 2.82 1.88 2.47 2.51 2.96 1.63 2.02
12 14 11 8 10 9 12 7 8
1.50 1.66 1.46 1.28 1.50 1.54 1.50 1.10 1.34
6 7 6 5 6 6 6 4 5
72 68 64 81 71 66 82 71 73
0.95 0.93 0.96 0.99 0.97 0.97 0.96 0.99 0.98
10 5 2 10 5 2 10 5 2
2.17 2.16 3.05 1.70 1.84 1.26 2.22 1.67 1.19
9 8 12 7 7 5 9 7 5
1.30 1.32 1.10 1.31 1.19 0.98 1.48 0.96 0.73
5 5 4 5 5 4 6 4 3
65 64 58 71 66 58 77 64 58
0.98 0.98 0.95 0.99 0.99 1.00 0.98 0.99 0.99
M2
M3
D3
M1
M2
M3
Relative RMSE (%)
Relative Bias (%)
* Relative values were divided by the average diameter of all field measurements – for all segments. Fig. 10. Cluster analysis of diagnostic variables for all different denoising and modelling combinations, point cloud densities and tree species (c = cophenetic correlation).
denoising and modelling methods, close attention must be paid during the scanning to retrieve good quality data. In extreme cases the tree may appear deformed in the point cloud, yielding highly overestimated diameters on some of the automated methods. This effect can be seen in Fig. 8 (D2-M2, D2-M3 and D3-M2) as a cluster of points pertaining to a single Eucalypt tree. Liang et al. (2016) made a compilation of studies and results available for dealing with TLS data from forest inventories – mostly from temperate or boreal regions. Studies of TLS applications in the tropics are scarce, but with data becoming more easily available, the amount of studies is expected to rise in these areas shortly. The most robust denoising method was the Hough transformation (D1), as it preserved stem features further up in the canopies, also prompting better accuracy on the subsequent stem modelling phase, due to more well preserved stem point clouds. The best combination was denoising by Hough transformation (D1), followed by the IRTLS cylinder fit algorithm (M2) for stem modelling – although combinations involving M2 and M3 performed similarly, as seen in Table 2. The cylinder modelling approaches have better performance approximating stem features, as they can estimate the angles at which stem segments lean from a straight, vertical bole, while the circle based approach assumes perfect verticality of the stem, estimating only the segments’ radii on a horizontal plane. Combinations involving M2 and M3 performed similarly. Those slight differences may disappear with finer parameter tuning. For every
than cylinder based methods (RMSE: 3.7–9.17 cm). Olofsson et al. (2014) found diameter RMSE values in the range of 2.4–7.5 cm. They tested their adapted Hough transformation denoising and circle fit RANSAC algorithm on boreal tree species point clouds, obtaining best results for Scots pine (RMSE = 2.38 cm, bias = -0.6 cm), followed by Norway spruce (RMSE = 6.18 cm, bias = 1.54 cm) and deciduous species (RMSE = 7.48 cm, bias = 2.52 cm). Liang et al. (2014a, 2014b) also tested their algorithms (eigen decomposition and iterated reweighted least squares cylinder fit) on boreal species, finding a bias of 0.15 cm and RMSE values of 1.17 cm for Pine and 1.08 cm for Spruce. Furthermore, they also reported that automated estimations were as precise as visual assessment of the point cloud, and around 3 times less biased. The reported diagnostic measures were also obtained by contrasting their methods with manual field measurements of diameters along the stems. Most methods implemented on this study performed better on pine trees, followed by spruce and then eucalypt. Pine trees have less dense lower canopy, while spruce stems are often occluded by branches at low heights. Performance on eucalypt trees was poor in comparison, mainly due to wind effect during scanning given the sideways movement of the trees, which provoked the effect of “ghost stems” in the point cloud. In addition, the eucalypt stand was only seven years old, and the trees thin and slender, which explains the high relative values of RMSE and bias observed in Table 2. Our implemented algorithms were not able to properly model those trees, so in order to get good accuracy from stem
174
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
clusters. Testing nine different algorithm combinations we observed that estimations made on stem segments below 10 cm of diameter were generally much less accurate or not even detected. Without density reduction at higher portions of the canopy in the point cloud, high accuracy is likely achievable for even smaller diameters, decreasing this threshold. Since the most valuable, commercial wood is of larger diameters, concentrated on the first few meters above the ground, the modelling techniques have proven to be effective for measuring commercial wood on standing trees. From a biomass and carbon stock estimation perspective, more robust and refined stem extraction methods and reliable taper models are required to attain precise and accurate estimations of upper diameters from TLS point clouds, as well as increased level of detail on the tree modelling, including branches and possibly foliage. Some methods for full tree modelling are described by Dassot et al. (2012), Raumonen et al. (2013) and Hackenberg et al. (2015). Even though the data used here was obtained through terrestrial laser scanning, the algorithms tested can be applied on point clouds obtained using other technologies, such as mobile scanning and photogrammetric methods (Liang et al., 2014a). As the algorithms rely solely on geometrical features, requiring only 3D coordinates as input, any data source that provides this type of data can be targeted. The constant increase in popularity of LiDAR devices, which is associated to decreasing market prices, will make 3D point cloud data easily available and cheaper for a wider range of customers, therefore availability of robust, time efficient and flexible software for data processing is already a demand.
method implemented, both denoising and modelling, a set of general parameters was chosen based on empirical tests, on which the algorithms were capable to properly detect and subsequently model stems from varying tree architectures. There is no best global parameter combination for a given method and default parameters should work well on a wide variety of tree species and ages. Nevertheless, the parameter choice is likely optimizable for distinct target tree species/ architectures for all combinations, in order to improve robustness and accuracy for specific conditions, such as dense canopy or scattered branching. In this study, we opted for a default general parameter setting per method, regardless of tree species, in order to investigate the performance of fully automated procedures. Parameter tuning assessment under different forest conditions should be addressed in a simulation study. The RANSAC circle fit (M1) was the most stable stem modelling method, as its RMSE and bias values fluctuated the least. By taking a closer look at Table 2 it becomes clear that this happens due to the eucalypt trees. This method could process more robustly point clouds of worse quality, especially when applied in sequence from the Hough transformation denoising algorithm (D1). Diameter estimations above half of tree total height were less precise and accurate due to the longer distance from the laser sensor and occlusion by branches and the canopy layer. Larger laser footprint, as well as a sharper angle between laser beams and stem surfaces closer to the tree tops are also a possible source of error in the estimations. Some possible solutions to get further estimations of missing stem segments are regression, interpolation or species-specific taper equations based on the estimated values from the TLS data. Regarding some specifics of the denoising techniques, D1 is more robust but its performance is hindered in case of tilted trees and forking stems. The Hough transformation works well on slightly tilted trees, but in order to adapt the method to detect forked stems, extra parameters and conditions must be added to the algorithm. D2 and D3 overcome those issues by detecting plane-like features whose points can deviate from a strictly vertical distribution, regardless of stem forking. Those methods need denser point clouds at the canopy region though, as they are more prone to omit stem segments with scarce inliers on the stem surface. Regarding the modelling algorithms, M1 requires least amount of points, being less affected by absence of points at the stem surface. Its main issue relates to bent or tilted trees, as it will often overestimate the diameter for any stem segment whose section is not exactly horizontal. For generally straight trees this method is precise and robust, but for other cases, methods that extract a cylinder axis’ direction are better suited, such as M2 or M3. Those methods require connectivity among stem surface patches and a larger amount of points to yield good estimations of cylinder parameters, thus being more sensitive to point cloud density at the canopy layer. Fitting rings on a two-dimensional plane is more robust in this regard, since it is less demanding in terms of data amounts and patch connectivity, thus possibly a better choice for processing tree point clouds acquired by low cost or low resolution scanners. All methods were applied on reduced point density clouds, with no more than 3x105 points per tree, having little to no loss in terms of robustness and accuracy estimating diameters along the stems. Kelbe et al. (2015) estimated dbh from single scan, low density TLS point clouds, achieving a RMSE of 6.4 cm. According to Liang et al. (2014b), further improvement on TLS-based estimation of tree attributes could be achieved by using higher point density. Our study shows that such algorithms can provide accurate estimates on lower point densities, and rather than increasing point cloud density globally, adopting different point densities per height strata may suffice the goal of improved global precision/accuracy – e.g. by maintaining the original point density only at the canopy layer, avoiding missing stem segments. Fitting circles or cylinders on narrow features can be tricky, since the algorithms tend to falsely identify those shapes in small point
5. Conclusions All tested methods are accurate and precise on single tree point clouds of good quality. They perform generally better at lower portions of the stem. The Hough transformation (D1) followed by iteratively reweighted total least squares cylinder fit method (M2) presented best combined performance of denoising and stem modelling, respectively. D1 consists in a more robust denoising method as it is less sensitive to outliers, being able to capture higher percentages of stem features at upper tree heights. The cylinder based modelling methods (M2 and M3) have better modelling accuracy as they consider leaning angles in addition to radii estimations, thus being less prone to diameter overestimation. For roughly straight trees, a combination of D1 and M1 is recommended for robustness and stability, respectively. For varying tree architectures a cylinder modelling approach is recommended (M2, M3), and in case of forking stems, denoising methods D2 and D3 are better suited. All functions developed during this study are available as an open source R package on GitHub: github.com/tiagodc/TreeLS. Currently the package provides a collection of functions for manipulation of TLS point clouds, stem mapping at plot level, automatic stem isolation from single trees and taper extraction. Acknowledgments We would like to thank the Sven and Hildur Wingquist Foundation For Forest Research, for financing the Swedish data acquisition, managed by the forest remote sensing section at SLU Umeå. The Erasmus Mundus programme in Sustainable Forest and Nature Management, SUFONAMA. The Southern Swedish Forest Research Centre, SLU Alnarp. References Adler, D., Murdoch, D., 2016. rgl: 3D Visualization Using OpenGL. R Package Version 0. 96.0. < https://CRAN.R-project.org/package=rgl > . Aschoff, T., Spiecker, H., 2004. Algorithms for the automatic detection of trees in laser scanner data. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 36, 71–74.
175
Computers and Electronics in Agriculture 143 (2017) 165–176
T. de Conto et al.
Lukács, G., Martin, R., Marshall, D., 1998. Faithful least-squares fitting of spheres, cylinders, cones and tori forreliable segmentation. Proc. Eur. Conf. Comput. Vis. 1, 671–686. Maas, H., Bienert, A., Scheller, S., Keane, E., 2008. Automatic forest inventory parameter determination from terrestrial laser scanner data. Int. J. Remote Sens. 29, 1579–1593. Mengesha, T., Hawkins, M., Nieuwenhuis, M., 2015. Validation of terrestrial laser scanning data using conventional forest inventory methods. Eur. J. For. Res. 134, 211–222. Nelder, J.A., Mead, R., 1965. A simplex method for function minimization. Comput. J. 7 (4), 308–313. Olofsson, K., Holmgren, J., Olsson, H., 2014. Tree stem and height measurements using terrestrial laser scanning and the RANSAC algorithm. Remote Sens. 6 (5), 4323–4344. Pfeifer, N., Winterhalder, D., 2004. Modelling of tree cross sections from terrestrial laser scanning data with free-form curves. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 36, 76–81. Pueschel, P., Newnham, G., Rock, G., Udelhoven, T., Werner, W., Hill, J., 2013. The influence of scan mode and circle fitting on tree stem detection, stem diameter and volume extraction from terrestrial laser scans. ISPRS J. Photogramm. Remote Sens. 77, 44–56. R Core Team, 2016. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
Brolly, G., Király, G., 2009. Algorithms for stem mapping by means of terrestrial laser scanning. Acta Silv. Lignaria Hung. 5, 119–130. Choi, S., Kim, T., Yu, W., 2009. Performance evaluation of RANSAC family. Br. Mach. Vis. Assoc. 81. Chum, O., 2005. Two-view geometry estimation by random sample and consensus. Ph.D. Thesis. Czech Technical University, Prague, Czech. Dassot, M., Colin, A., Santenoise, P., Fournier, M., Constant, T., 2012. Terrestrial laser scanning for measuring the solid wood volume, including branches, of adult standing trees in the forest environment. Comput. Electron. Agric. 89, 86–93. Forsman, P., Halme, A., 2005. 3-D mapping of natural environments with trees by means of mobile perception. IEEE Trans. Robot. 21, 482–490. Gorte, B., Pfeifer, N., 2004. Structuring laser-scanned trees using 3D mathematical morphology. Int. Arch. Photogramm. Remote Sens. 35, 929–933. Gorte, B., Winterhalder, D., 2004. Reconstruction of laser-scanned trees using filter operations in the 3D raster domain. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 36, 39–44. Hackenberg, J., Morhart, C., Sheppard, J., Spiecker, H., Disney, M., 2014. Highly accurate tree models derived from terrestrial laser scan data: a method description. Forests 5 (5), 1069–1105. Hackenberg, J., Morhart, C., Sheppard, J., Spiecker, H., Disney, M., 2015. SimpleTree – an efficient open source tool to build tree models from TLS clouds. Forests 6 (12), 4245–4294. Henning, J.G., Radtke, P.J., 2006. Detailed stem measurements of standing trees from ground-based scanning LiDAR. For. Sci. 52, 67–80. Hilker, T., Coops, N.C., Culvenor, D.S., Newnham, G., Wulders, M.A., Bater, C.W., Siggins, A., 2012. A simple technique for co-registration of terrestrial LiDAR observations for forestry applications. Remote Sens. Lett. 3 (3), 239–247. Hopkinson, C., Chasmer, L., Young-Pow, C., Treitz, P., 2004. Assessing forest metrics with a ground-based scanning LiDAR. Can. J. For. Res. 34, 573–583. Illingworth, J., Kittler, J., 1987. The adaptive Hough transform. IEEE Trans. Pattern Anal. Mach. Intell. PAMI-9 (5), 690–698. Kaufman, L., Rousseeuw, P., 1990. Finding Groups in Data: An Introduction to Cluster Analysis. John Wiley & Sons Inc, New York. Kelbe, D., van Aardt, J., Romanczyk, P., van Leeuwen, M., Cawse-Nicholson, K., 2015. Single-scan stem reconstruction using low-resolution terrestrial laser scanner data. IEEE J. Sel. Top. Appl. Earth Observ. Remote Sens. 8 (7), 3414–3427. Lalonde, J.F., Vandapel, N., Huber, D., Hebert, M., 2006. Natural terrain classification using three-dimensional ladar data for ground robot mobility. J. Field Robot. 23, 839–861. Liang, X., Litkey, P., Hyyppä, J., Kaartinen, H., Vastaranta, M., Holopainen, M., 2012. Automatic stem mapping using single-scan terrestrial laser scanning. IEEE Trans. Geosci. Remote Sens. 50 (2), 661–670. Liang, X., Jaakkola, A., Wang, Y., Hyyppä, J., Honkavaara, E., Liu, J., Kaartinen, H., 2014a. The use of a hand-held camera for individual tree 3D mapping in forest sample plots. Remote Sens. 6, 6587–6603. Liang, X., Kankare, V., Yu, X., Hyyppä, J., Holopainen, M., 2014b. Automated stem curve measurement using terrestrial laser scanning. IEEE Trans. Geosci. Remote Sens. 52 (3), 1739–1748. Liang, X., Kankare, V., Hyyppä, J., Wang, Y., Kukko, A., Holopainen, M., Vastaranta, M., 2016. Terrestrial laser scanning in forest inventories. ISPRS J. Photogramm. Remote Sens. 115, 63–77. Lindberg, E., Holmgren, J., Olofsson, K., Olsson, H., 2012. Estimation of stem attributes using a combination of terrestrial and airborne laser scanning. Eur. J. For. Res. 131, 1917–1931.
176