VOXEL-BASED APPROACH FOR ESTIMATING URBAN TREE VOLUME FROM TERRESTRIAL LASER SCANNING DATA

The importance of single trees and the determination of related parameters has been recognized in recent years, e.g. for forest inventories or management. For urban areas an increasing interest in the data acquisition of trees can be observed concerning aspects like urban climate, CO2 balance, and environmental protection. Urban trees differ significantly from natural systems with regard to the site conditions (e.g. technogenic soils, contaminants, lower groundwater level, regular disturbance), climate (increased temperature, reduced humidity) and species composition and arrangement (habitus and health status) and therefore allometric relations cannot be transferred from natural sites to urban areas. To overcome this problem an extended approach was developed for a fast and non-destructive extraction of branch volume, DBH (diameter at breast height) and height of single trees from point clouds of terrestrial laser scanning (TLS). For data acquisition, the trees were scanned with highest scan resolution from several (up to five) positions located around the tree. The resulting point clouds (20 to 60 million points) are analysed with an algorithm based on voxel (volume elements) structure, leading to an appropriate data reduction. In a first step, two kinds of noise reduction are carried out: the elimination of isolated voxels as well as voxels with marginal point density. To obtain correct volume estimates, the voxels inside the stem and branches (interior voxels) where voxels contain no laser points must be regarded. For this filling process, an easy and robust approach was developed based on a layer-wise (horizontal layers of the voxel structure) intersection of four orthogonal viewing directions. However, this procedure also generates several erroneous “phantom” voxels, which have to be eliminated. For this purpose the previous approach was extended by a special region growing algorithm. In a final step the volume is determined layer-wise based on the extracted branch areas Ai of this horizontal cross-section multiplied by the thickness of the voxel layer. A significant improvement of this method could be obtained by a reasonable determination of the threshold for excluding sparsely filled voxels for noise reduction which can be defined based on the function change of filled voxels. Field measurements were used to validate this method. For a quality assessment nine deciduous trees were selected for control and were scanned before felling and weighing. The results are in good accordance to the control trees within a range of only -5.1% to +14.3%. The determined DBH values show only minor deviations, while the heights of trees are systematically underestimated, mainly due to field measurements. Possible error sources including gaps in surface voxels, influence of thin twigs and others are discussed in detail and several improvements of this approach are suggested. The advantages of the algorithm are the robustness and simple structure as well as the quality of the results obtained. The drawbacks are the high effort both in data acquisition and analysis, even if a remarkable data reduction can be obtained by the voxel structure.


INTRODUCTION
Until recently, forest stands have been the main focus concerning national carbon inventories and conservation, but recently the importance of urban trees is being increasingly investigated.Urban habitats differ significantly from natural systems (Norra, 2009) due to site conditions (e.g.contaminants, regular disturbance, climate, and species composition).In contrast to forest stands urban trees are planted for ecological and social services for the local population (Konijnendijk et al., 2005).Therefore, the urban tree species composition is dependent on site conditions, the arrangement of recreational areas and (gardeners) experience and taste.Consequentially, urban trees differ in structure from those grown in forest stands, and allometric relations developed for forest situations cannot be applied to the urban environment.New allometric relations must be developed for a rapid and easy measurement of urban tree biomass (Zianis, Mencuccini, 2004).For this purpose an extended approach was developed for a fast and non-destructive determination of the branch volume of a tree, DBH (diameter at breast height) and tree height in urban environment from point clouds of terrestrial laser scanning (TLS).In the next chapter, related literature is compiled and commented, followed by a section about the used test material, data acquisition and preprocessing.In chapter 4, the data analysis is described in detail; results of practical application including a quality assessment are shown in chapter 5, while conclusions and an outlook are given in chapter 6.

RELATED WORK
A common and well-known task is the estimation of forest carbon stocks.Good et al. (2001) present a compilation for estimating the biomass of individual trees, including double regression, ratio sampling, Randomized Branch Sampling (RBS) and Importance Sampling (IS).The most common method is the application of allometric functions (e.g.Zianis, Mencuccini, 2004).Due to the potential of terrestrial laser scanners to rapidly and accurately acquire suitable data for estimation of tree volume (Lefsky, McHale, 2008), several research groups are currently adapting and enhancing this approach.Based on techniques of Computer Graphics and Remote Sensing (an overview is given in Bucksch et al. (2010)) the aims of this research are volume determination (e.g.Lefsky, McHale, 2008;Raumonen et al., 2011;Dassot et al., 2010), structural and radiative consistency (Côté et al., 2009), estimation of branch parameters, including length, radii or straightness of stem and branch (Bucksch, Lindenbergh, 2008;Bucksch et al., 2010;Gorte, Pfeifer, 2004a,b;Pfeifer, Winterhalder, 2004;Pfeifer et al., 2004), branch size distribution (Raumonen et al., 2011), dimensionality of tree stems (Gatziolis et al., 2010), angles between branches, and the reconstruction of the outer hull (Pfeifer et al., 2004).In addition to measurements of diameters, three basic approaches can be found: (i) extraction of a skeleton of a tree (Bucksch, Lindenbergh, 2008;Bucksch et al., 2010;Gorte, Pfeifer, 2004a,b;Côté et al., 2009;Xu et al., 2007), (ii) fitting predefined primitives, e.g.cylinders (Pfeifer et al., 2004;Lefsky, McHale, 2008;Raumonen et al., 2011), free-form curves (B-Splines) (Pfeifer, Winterhalder, 2004;Yan et al., 2009) or mesh in general sense (Xu et al., 2007), (iii) structuring the point clouds in a voxel domain (Gorte, Pfeifer, 2004a,b;Gatziolis et al., 2010).The main challenges in this field of research are the large data requirements (varying from several thousand to commonly 2 million points), noise (rarely treated, e.g.Gorte, Pfeifer, 2004a,b, but discussed as error source, e.g.Bucksch et al., 2010), varying point densities (may lead to over-or underestimation of tree volume, e.g.Leksky, McHale, 2008), gaps or holes, the vast variety of branch diameters and lengths with a tree and a reliable estimation of complete stem and crown volumesrarely presented in these approaches.

DATA ACQUISITION AND PRE-PROCESSING
The test site was located at the urban area of Karlsruhe city in the Upper Rhine Valley, Germany.From about 560 trees scheduled to be cut and weighed by the Gardening Department of Karlsruhe (GBA) during winter period 2009/2010, nine deciduous trees were chosen as suitable control data for these investigations.Additional aspects were a low level of interpenetration of the crown with adjacent trees, the absence of leafs, fruits, ivy or mistletoe and a variety of different tree species (Acer spp.(denoted by A1, A2, A3), Tilia spp.(T1, T2), Robinia spp.(R1, R2), Betula spp.(B1) and Carpinus betulus (C1); cf.Table 1).Before felling, each tree was scanned from three to five positions with a Leica HDS 6000 using the finest scan incrementation of approx.2mm at an average distance of 10m and a point accuracy of about ±3mm.To obtain a highly accurate registration of these individual point clouds (using Cyclone 5.8 ® ), five sphere targets had been positioned around each tree.The standard deviation ranges from ±1mm to ±2mm for all investigated trees.To convert green weight at felling time to dry volume and carbon content several wood samples at different heights were collected and analysed after felling.
A point cloud reduction using an average point spacing of 2mm was also carried out (Cyclone "unify" procedure) as well as a manual elimination of disturbing objects like traffic signs, parked cars or branches of neighbouring trees.

DATA ANALYSIS
Data analysis starts with the generation of the voxel structure and the assignment of the laser points to this voxel domain.Voxels without points are not considered in the subsequent analysis.The necessary definition of a suitable voxel size depends on the applied scan resolution of the terrestrial laser scanner and the dimension of the finest branches to be measured.A certain number of points per voxel should be achieved, e.g. with a voxel edge length of 1cm like in this application and diameters of finest branches of 3-4mm about 20 points per voxel will be obtained, which can be subsequently distinguished from voxels containing noise or so-called "phantom" points, well-known in the domain of terrestrial laser scanning.Based on this voxel structure a specific approach was developed which analyses the data in horizontal layers of one voxel thickness.

Noise Reduction
Two different kinds of noise reduction are carried out: (i) elimination of isolated voxels and (ii) elimination of voxels with a small number of laser points.In the first case those voxels should be excluded which are not connected to the tree, i.e. which are isolated (cf.Gorte, Pfeifer, 2004a).Based on a N26neighbourhood (comparable to a N8-neighbourhood in 2D) or a N124-neighbourhood, the central voxel is eliminated if all adjacent voxels are empty, because it implies that there is a gap between the surface of the tree and the particular points.
A second method of noise reduction was applied to eliminate voxels which are sparsely filled by laser points, i.e. voxels containing fewer points than a pre-defined threshold.For rough volume estimations a threshold may be defined by experience, but our investigations have shown that this threshold varies dependent on the tree type and therefore must be determined by a specific method explained in section 4.4.

Filling of Hollow Surfaces (Interior Voxel)
One challenge of voxel-based approaches for determination of tree volumes is the determination of hollow surfaces, i.e. voxels lying inside the stem and branches which contain no laser points but have to be taken into account for volume calculations.Mathematical morphology (dilation, closing) may be a possible method (e.g.Gorte, Pfeifer, 2004a), but the problem in this application is that the dimension of the hollow surface exceeds the space between narrow branches leading to erroneously filled areas which cannot be reversed (e.g. by erosion).Therefore, another procedure was chosen based on the intersection of the four areas of obstructed vision from four orthogonal directions.It is robust, fast, does not need prior information, and fits to all forms of branches and stems.Additionally, repeated iterations (cf.mathematical morphology) are not required.Within one horizontal voxel layer, the algorithm searches the main four planimetric directions (coordinate directions) for surface voxels.If such voxels are detected all subsequent voxels in this row or column are set to value v=1, except the original surface voxels (Figure 1b).After executing this procedure for all four directions (Figure 1b-e) and intersecting the areas of obstructed vision (v=1), i.e. adding up the voxel values, all interior voxels contain the value v=4 (Figure 1f) and can be marked as "filled" voxels whereas all voxel values below 4 are set back to v=0.Therefore, it can be distinguished between original surface voxels, artificially "filled" voxels (v=4) and background voxels (v=0).Erroneously filled voxels are most likely to occur in tree crowns with narrow branches and twigs (Figure 2a).To solve these problems an extended approach was created (based on a region growing algorithm) for segmentation of regions of filled voxels.Starting at a filled voxel the N4 neighbourhood (due to the four viewing directions) is analysed.
If a neighbouring voxel is also a filled voxel it is accepted and the growing process is continued until the border of a filled region is reached.Then the border voxels are checked if at least one unfilled neighbouring voxel (background voxel) can be found, i.e. not exclusively filled or surface voxels.In this case the whole filled region is eliminated, i.e. its values are set to v=0 (Figure 2b).Small gaps often remain at the acquired surfaces of stem or thicker branches caused by occlusion effects of surrounding branches (even if scanning from four or five different directions) leading to the erroneous omission of real tree volume.Therefore, this rigorous algorithm could be modified by accepting a filled region, even if not all but most border voxels have an adjacent surface voxel (e.g.95% or 90%).

Layer-wise Volume Estimation
The estimation of the layer-wise tree volume is based on the segments generated in the previous filling process described above.Assuming that the diameter of a branch does not change significantly in a horizontal layer of one voxel thickness, the volume can be approximated by an oblique cylinder (Figure 3), i.e. by the cross-sectional area A i (=segment area) of a branch multiplied by height h i (=layer thickness).Due to the fact that surface voxels are normally not completely filled by laser points, the total area A i of a cross-section, i.e. the sum of the area of surface voxels plus filled voxels, overestimate the real volume.On the other hand regarding only the filled (interior) voxels the real volume will be underestimated.Therefore, both areas are averaged for the volume estimation.The final tree volume results from the sum of each branch volume of all layers.The diameter at breast height, DBH, is derived from averaging the cross-sectional areas, A i , of the stem between a height of 1.1m and 1.5m.Tree height is determined as the difference between the lowest and the uppermost voxel layer of the point cloud.

Determination of Threshold for Noise Reduction
A detailed inspection of histograms of voxel point densities show that a high number of voxels contain only a few points.Due to the scan resolution and the mean scanning distance of about 10m even small twigs are acquired by approximately 10-15 points per voxel (edge length 1cm).Taking occlusion effects into account this number may be reduced to 8-12 points per voxel.It can therefore be assumed that voxels containing fewer points can be stated as noise.To define a suitable threshold the function "change in filled voxels" (Figure 4) can be analysed.For this purpose, the method described in section 4.2 is iteratively applied starting from 1 point per voxel for elimination.With increasing number of points per voxel, noise is eliminated and, in consequence, a lot of erroneously filled voxels disappear and a bigger change in filling is obtained.Due to the fact that real twigs and branches have a certain minimum number of points per voxel, there is no larger change in filling by increasing this threshold until reaching this number.Exceeding this number of points per voxel, real surface voxels are eliminated and the surfaces are increasingly perforated, which leads again to a larger change of fillings.Therefore, a suitable threshold can be defined as minimum of the function "change in filled voxels" (Figure 4).This method delivers appropriate results in most of the investigated cases.In only one exception however, this function shows a continuous decrease of filled voxels without a minimum of change.In this case the threshold for elimination of noise voxels is set to a feasible minimum value, e.g.t=1 (all voxels with only one point are eliminated).

RESULTS AND QUALITY ASSESSMENT
The scans of nine deciduous trees resulted in point clouds of 20 to 60 million points.After applying both methods of noise suppression a strong data reduction could be achieved.The presented approach of volume calculation yielded good results ranging between -5.1% and +14.3% compared to the control volumes derived from the weights of the felled trees (Table 1) with a slight trend to overestimation for trees with dense structure of twigs.Nevertheless, certain differences could be detected.One problem of tall trees in this study is the metrological reduction of scan density with increasing heightand thus distance.Small twigs at the top of a tree would require a particularly higher scan density.Another problem is a significant amount of noise caused by different determining factors.On the one hand, a tree is not a totally static object, i.e.
due to the sequential scanning mode inevitable movements of branches and twigs lead to noise of the acquired surface point.
On the other hand, noise (and so-called "phantom" points) is introduced by the large numbers of edges where the laser beams are split, introducing errors in distance measurements and reducing the accuracy from millimetres to centimetres.This effect occurs mainly in trees with small ending twigs and a high degree of arborisation, and is further exacerbated by disturbing objects like thorns, dried leaves or fruits.Additionally, registration accuracy decreases with increasing tree height due to the limited spatial distribution of the targets (Pfeifer et al., 2004).Another major influence concerns occlusion and shadowing of interior branches dependent on the degree of arborisation (even if scanning from 4 or 5 directions), which,especially in this approach,leads to missing filled voxels.In addition to tree volume, DBH and the height of treesoften used as variables for biomass functionsare calculated (Table 2).Estimates of DBH correspond to the measured control values with only marginal deviations (with one exception).For all trees, height determination is systematically lower than related manual measurements.Approximately, the differences lie between 0.5m and 1.5m caused by the uncertain acquisition of the highest part of crowns from terrestrial positionby TLS as well as by manual measurement.

CONCLUSIONS
The presented method delivers good results in accordance with the control data.It was found to be more accurate than other methods such as detailed analysis of aerial images or measurement of main branch perimeters by specialists climbing the trees and predicting the connected branch volumes.The usage of a voxel structure, the noise suppression and the processing of the data in horizontal layers leads to a significant data reduction and enables the processing of large data sets consisting of 20 to 60 million points.The applied extended algorithm for filling hollow surfaces is generally applicable for stems as well as for complex crowns.The advantages of this method are its simplicity and robustness, i.e. no complex skeleton extraction or fitting of cylinders is necessary.
Moreover, volume extraction is independent of branch length, diameter or orientation.It needs no advanced input information and can be used for a wide range of applications.The drawbacks of this approach are the necessity of high resolution scans and the extensive computing time for determination of a suitable threshold for noise reduction.In future studies, it is suggested that the number of analysed trees be increased to obtain a better statistical basis, to allow extrapolation to larger areas.Additionally, further investigations have to be carried out concerning the influence of the number and position of necessary scans, possibly in dependence on the tree structure.

Figure 1 .
Figure 1.Principle of filling process, a) Original surface voxels of a horizontal layer, b)e) Marking of voxels of obstructed vision in four orthogonal directions, f) Intersection result of the four areas of obstructed vision, only voxels containing value v=4 are filled voxels for volume estimation

Figure 2
Figure 2. a) Correctly filled interior voxels and erroneously filled voxels between branches, b) Result after elimination of erroneously filled voxels / voxel segments by a region growing approach analysing the borders of these segments

Table 1 .
Estimated and control volume of scanned treesWhile the stems and thicker branches are usually completely filled, filling often failed at higher parts of the crown, independent of the orientation of the branches.The degree of correct fillings also depends on the appropriate suppression of noise.The presented method complies with theoretical expectations; only one tree shows a different behaviour (continuous decrease of the number of filled voxels).

Table 2 .
Estimated and control DBH and tree heights