RANGE IMAGE SEGMENTATION FOR TREE DETECTION IN FOREST SCANS

To make a tree-wise analysis inside a forest stand, the trees have to be identified. An interactive segmentation is often labourintensive and time-consuming. Therefore, an automatic detection process will aspired using a range image. This paper presents a method for the segmentation of range images extracted from terrestrial laser scanner point clouds of forest stands. After range image generation the segmentation is carried out with a connectivity analysis using the differences of the range values as homogeneity criterion. Subsequently, the tree detection is performed interactively by analysing one horizontal image line. When passing objects with a specific width, the object indicates a potential tree. By using the edge points of a segmented pixel group the tree position and diameter is calculated. Results from one test site are presented to show the performance of the method.


INTRODUCTION
Terrestrial laser scanning became more important for forest applications in the last few years.Where a scan is done within few minutes, the post processing of the scan data is a lengthy task.The tree position and even segmented trees are of interest to seek a tree wise evaluation as for example tree height determination, extraction of the diameter at breast height (DBH) (Aschoff et al., 2004, Henning & Radtke, 2006;Watt et al., 2005;Maas et al., 2008, Wezyk et al., 2007;Hopkinson et al., 2004) and skeletonisation of the branch structure (Bucksch & van Wageningen, 2006, Schilling et al., 2012).
Separation of point clouds is a time consuming process when an interactive separation is done.The 3D point clouds have to be moved to different views and points that do not belong to the tree have to be removed (Hopkinson et al., 2004).An easy and automatic method is the separation using vertical cylinders inside the 3D point cloud (Maas et al., 2008).Assuming that the tree position is given, a cylinder with a constant radius is placed at this position.However, a big disadvantage is the constant radius which cuts the branches in the case they are oversized.Other objects which are included inside the cylinder are also separated.
The aim of the work presented here is the automatic separation and detection of trees inside a forest stand using a range image created from one single scan.A method will be shown to segment a range image and perform stem detection.Advantageous is the identification of occluded regions inside a 2D view rather than using 3D point clouds.Neighbouring pixels with a big difference in range value indicate two separate objects, where the farthest one is partly hidden with a high probability.
The outline of the paper is the following: Section 2 gives an overview of the data set used.In Section 3 the algorithm is shown, starting with a separation of the ground points, the range image creation and the segmentation.Subsequently the tree detection follows in Section 4. Results from one data set will be presented in Section 5. Finally a discussion will be given and a summary and an outlook close the article.

DATA SET
One data set is chosen to evaluate the presented methods for tree segmentation.The site are recorded with a Riegl LMS-Z 420i laser scanner on a research site located in the Tharandter Wald, about 25 km south-west of the city of Dresden in Germany (50°57′49′′N, 13°34′01′′E and 380 m a.g.l.).The laser scanner has a range s of 2 ≤ s ≤ 800 m and average measurement accuracy between 5 mm and 7 mm (Riegl, 2009).
The data set is one of four ground based scans of a 147 year old beech stand (Fagus sylvatica) with a sloping terrain (6°) and was recorded in March 2007 (leaf-off).The test site (Figure 1) has not much ground vegetation.Further stand features are given in Spank (2010) and Bienert (2013). (2)

ALGORITHM
The algorithm is divided into 4 steps: 1) separation of the ground points; 2) creation of the range image; 3) connectivity analysis with the adjacent pixels; 4) generation of the segmented objects as 3D point clouds.

Separation of the ground points
A first and crucial step is the separation of the ground points, due to the fact that the segmentation will be carried out with a connectivity analysis of adjacent pixels (Bienert et al., 2010;Douillard et al., 2011).This method prevents that the trees are interconnected via the ground pixels.It is possible to delete the points interactively or in automatic process.
The automation process operates as follows: First of all, a digital terrain model (DTM) with a predefined grid size s DTM is generated with a vertical point density analysis.The Z coordinates of all points inside the grid are listed concerning their frequency.In the height of the ground, the highest point density will be assumed.Therefore the first peak inside the histogram indicates the height of the ground and thus the height of the grid.
The next step is the transformation of the whole laserscanner point cloud into a voxel space with a voxel size s voxel .The bounding box of the point cloud and the voxel size define the dimension of the voxel space.Voxels which contain at least one laserscanner point have a voxel value of 1 (filled voxel) and voxels without points (empty voxel) a value of 0.
Afterwards, the DTM raster points are also converted into the same voxel space.Each voxel which contains a DTM raster point ('DTM voxel') belongs to a group of voxels ('ground segment') which will be later deleted including all contained points.Starting with an arbitrary 'DTM voxel', a 3D connectivity analysis is done only based on all 'DTM voxels'.All neighboured voxels, which are filled with laserscanner points, belong to the 'ground segment'.This connectivity analysis is only performed with all 'DTM voxels' for the central element of the structure element, to avoid a region growing to the lower parts of the tree stems.The DTM raster size s DTM and voxel size s voxel defines the order of neighbourhood of the structure element.Clearly, the first order neighbourhood is sufficient if the grid size of the DTM is smaller than the edge length of a voxel.However, once the DTM grid s DTM is larger than the voxel size s voxel , decides the ratio o N of the DTM grid size to the voxel size (Eq.( 1)), with which the neighbourhood connectivity analysis is performed to add the voxels to the 'ground segment'.
Finally, the containing laserscanner points of all voxels of the 'ground segment' will be separated automatically.The separated point cloud without ground vegetation is used in the further process.

Creation of range image
Producing a range image, the algorithm can be divided in three different steps according to (Vosselman & Klein, 2010): • definition of the image raster, • determination of the range value, • determination of the colour value of the pixel representing the range.

Image grid
The resolution (c,r) of the range image grid is given by the difference of the maximum and minimum scanning angles (φ, θ) of the scan scene divided by the scan resolution α.

Range and colour value
Avoiding gaps between the filled pixels the used scan resolution is lower than the original scan resolution.The range value is given by the 3D distance measurement to the point.In case of several points inside a pixel, the point with the smallest distance are stored as range value.Also a 2D distance can be used, built by the coordinates X and Y, assuming the scanner was levelled before scanning.Especially for forest applications, the type of distance is important, while producing a range image.A range image on the basis of a 2D distance shows more homogenous steps between two neighbouring pixels, especially in the upper stem sections (Figure 3a).This is advantageous in the segmentation of stems.
As pictured in Figure 3b) the range values (indicated through colour values) show larger steps in the upper regions of the stems which are caused by an increasing vertical angle and thus having a longer 3D distance.
For processing the image, a range matrix R with the dimension r x c of the image is build and filled with the transformed range value j i q , (Eq. ( 4)).j i q , with (i = 0,...r; j = 0,...c) is an integer variable calculated with the minimum distance inside the plot, the distance s i,j of the pixel on position i,j and the desired resolution s res .Parameter s res is a preset value given by the user and represents the distance resolution (default value 0.001 m).
Pixel with no points inside belong to the background.For visualisation purposes the colour is assigned by linear interpolations of the values between minimum and maximum distances.

Connectivity Analysis
After range image generation, as described in section 3.2, a second connectivity analysis is done on the basis of the image.Within a preset order of neighbourhood adjacent pixels are analyzed.To bridge occlusions, caused by twigs and branches in the upper stem part in front of the analysed object, a large structure element is used (e.g.27 px x 27px).
Whether a pixel belongs to the same object as its neighbour pixel is defined when a homogeneity criterion is fulfilled.The difference q  of the range value j i q , of the central pixel and its actual neighbour pixel decides whether the pixel affiliates to the object given by the central pixel.If q  is smaller than a threshold max q  then the neighbour pixel belongs to the object of the central pixel.Each pixel is checked for connectivity (Figure 4), while placing the kernel on every filled pixel.b) segmented range image (Bienert, 2013).
Due to the fact that the scanning procedure measures equiangular points, the spacing between two points increases with the distance to the scanner.Consequently, the difference q  of the range values of neighbouring pixels increases also with the distance to the scanner.To overcome this effect a floating threshold max , j i q  is used to assign the pixels to the same objects.This threshold is varied from pixel to pixel and depends on the distance of the pixel to the scanner.Equation ( 5) can be used for range images based on 3D distances as well as range images with 2D distances.If the last pixel is allocated to an object, the segmented pixels are transformed into a 3D point cloud (Eq.( 6)).Each segmented object is saved in a separate point cloud file.The transforming is done using following formulas: where X, Y, Z... coordinates in the scanner`s own coordinate system θ... vertical angle φ... horizontal angle α... scan resolution j... column i... row , ... 3D distance from scanner to pixel (i,j) , ... 2D distance from scanner to pixel (i,j) By using the calculated angles φ and θ, the coordinates of the 3D point could be falsified.The angles are not constantly caused by the scanning accuracy.Therefore the effect is noticeable in far distant regions.An alternative is the application of the original scan angles, which will be stored during the range image procedure in a matrix.

TREE DETECTION
The tree detection is performed on the basis of the segmented image.For that purpose a horizontal image line (interactively chosen) is analysed (Figure 5b).The changeover of coloured pixels indicates potential tree objects.Two different types of pixel crossings are distinguished:  from background (black) to an object (coloured) and the other way around,  from object to another object (colour changing).
Two neighbouring objects are an indicator for occlusion and so the following tree diameter determination of the hidden tree could be estimated to small.In that case, the hidden objects are marked for the further processing.Hidden objects are distinguished, comparing the range values j i q , of the edge pixels in case of a colour changing.A greater range value stands for an edge pixel of an assumed hidden object.Figure 5 shows an example of the steps of data processing after separation of the ground points (Figure 5a).The segmented range image contains 21 segmented objects (Figure 5b).By analysing one image line, 15 pixel crossings are recognized.(Bienert, 2013).
If a colour change is recognized, the edge points P 1 (X 1 ,Y 1 ) and P 2 (X 2 ,Y 2 ) are calculated for each edge pixel using Eq. ( 6).Using the edge pixels of an object (Faro, 2006) the diameter d ~and an approximate position tree P ~( tree s ~, φ tree ) can be determined.The 2D-distance given by both edge pixels defines the diameter d (Eq.( 8)).Obviously, d ~ is an approximation and it will always be determined too small.
The approximate position can be calculated by the triangle spanned by the scanner position, one edge point and the position tree P ~(Figure 6).The approximate tree position tree P ~is given by the distance tree s ~ and the horizontal angle φ tree : To decide if the object is a tree, two aspects have to be investigated:  visible objects (not hidden) must have a diameter d ~of > 3.5 cm,  partly hidden objects are stored always, even if the diameter d ~is smaller than 3.5 cm.

Separation of the ground points
For removal of all ground points, a DTM grid with a spacing s DTM of 0.5 m and a voxel space with a voxel size s voxel of 0.2 m were used.Applying Eq. ( 1) a neighbourhood o N of 2 was set to delete the ground points.Figure 7 pictures a detailed view of the point cloud with the DTM points (blue dots) and the filled voxels of the 'ground segment' (green boxes).The ground points inside the voxels are not shown, due to the fact of clarity of the figure.Some parts of the ground vegetation are preserved during the removal.This is caused, when the ground vegetation is higher than the order of neighbourhood of the structure element inside the voxel space.

Segmentation
The data set is a detailed scan from the ground inside a beech stand and represents an area of 154° ≤ φ ≤ 228° and 50° ≤ θ ≤ 111°.As a base for the image segmentation a 2D range image was used with a range between 5.14 m ≤ s ≤ 30.00 m.The image resolution is 2474 px x 2028 px with an angular pixel spacing of 0.03° (equivalent to the scan resolution).For generation the range matrix R a range resolution s res of 0.001 m was used (see Eq. ( 4)). Figure 8 shows the final 2D range image without ground points.
Due to the fact of heavy branching in the upper stem sections, a huge kernel (27 px x 27 px) was chosen to perform the segmentation.A fixed range r fixed of 3 x s res was used.Thus the homogeneity criterion max , j i q  (Eq.( 5)) reaches from 3 x s res up to 18 x s res from the minimum up to maximum range inside the picture.Figure 9 presents the segmented image.
Figure 9. Segmented range image with image line (white horizontal line at y'=1655 px) and 3 partly hidden trees (detected during the image line analysis, red arrows).The blue arrow marks the stem which is missed in the tree detection with the 'validation algorithm' (see section 5.3).

Tree detection validation
To validate the presented tree detection method, the point cloud was processed with another tree detection method ('validation algorithm').The algorithm searches point clusters inside a horizontal slice, which was filtered from the point cloud in a constant height above the digital terrain model.Finally, the detected clusters are analysed in a circle fitting procedure to ensure the candidate is a tree.The algorithm is shown in detail in (Bienert et al., 2007).The data set includes at breast height (1.3 m above terrain) 18 visible trees (complete and partly occluded).The 'validation algorithm' classified 18 objects as trees.17 of 18 were trees and one tree was not found.
The image line (y'=1655 px) contains 23 segmented objects (line position is shown in 9).From the outset 5 of them are classified as 'no trees', because of a diameter smaller than 3.5 cm.Three of the 18 remaining ones are partly hidden (see arrows in Figure 9).Finally 18 trees were found.Table 1 shows the detection results of both algorithms in detail.Besides the detection results, the diameter d, the distance s, the azimuth a and the differences between the diameter stem heights (ΔZ) are also shown.It is obvious in Table 1 that the distance s shows different minima.The reason is given by the missed tree (located at 6.29 m) in the 'validation algorithm'.However, both methods can be used in combination to increase the reliability of the tree detection.Table 1.Results of the tree detection compared with a 3D cluster search algorithm ('validation algorithm').

DISCUSSION
The success of the segmentation is coupled on the correct selection of the kernel size and the homogeneity criterion, which differs in dependence of the tree species and the scan resolution.A huge kernel bridges over areas which are partly hidden by small vegetation elements.As homogeneity criterion the smallest space between two neighbouring objects do not have to be undercut, to secure delimitation.A small kernel and a small homogeneity criterion lead to an oversegmentation.On the contrary a homogeneity criterion which allows a huge range difference produces undersegmented areas.Problems with the right segmentation occur especially in the crown regions, where the occlusion is too heavy.
Certain scanning conditions have to be fulfilled for a successful tree detection.First of all, the scanner has to be levelled to achieve vertical tree stems inside the range image.Also sloping terrain causes problems by using a fixed horizontal image line for tree detection.The effect is visible for far distant trees.If the slope is too heavy, the image line crosses the upper stem parts for downslope scan directions and runs below the stems (ground was removed) for upslope scan directions.
The data set shows a tree detection rate of 100%.Advantageous in comparison to the validation algorithm (working on a point cloud cut with a cluster searching) is the identification of partly hidden objects.Consequently, a statement is given if the diameter is estimated too small.Comparing the diameters, differences up to 15 cm are given.This effect is explainable by the different stem heights ΔZ passing the image line caused by the perspective.Three diameters were estimated below the reference diameters and 15 diameters above (up to 1.26 m).Finally all trees could be found with a good approximately position and diameter.

SUMMARY AND OUTLOOK
The paper presented a fast and effective method to detect trees inside a range image.The algorithm can handle different data sets, concerning the viewing directions and tree species.For an optimization of the segmentation results the homogeneity criterion and the pixel neighbourhood has to be adapted to the scan resolution and the width of the twigs, to overcome hidden parts.Furthermore the range image allows the utilization of image processing methods for the 3D scan data.It is shown, that the tree detection on base of the segmented range image, even for partly occluded stems, is working reliably.The estimated tree positions and diameters (of completely visible stems) might be comparable with the established tree detection methods, if the same stem height is used.
Future work deals with the estimation of diameter at breast height.Involving the DTM height, the breast height of stem can be determined, to make the diameter comparable with manual measured diameters gained from inventories.

Figure 1 .
Figure 1.Range image of the beech stand (range between 5.14 m ≤ s ≤ 30.00 m).

Figure 2 .
Figure 2. Principle of generation the 'ground segment':horizontal view of voxel space and DTM grid points with moving structure element with second order neighbourhood.

Figure 3 .
Figure 3. Detailed coloured range images of beech crowns recorded from the ground.a) 2D range image; b) 3D range image.

Figure 4 .
Figure 4. Principle of connectivity analysis with a 3 x 3px kernel: a) range image with range values j i q , ; the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey 4 (Figure5c) show a colour changing, whereas the other objects are surrounded from background.

Figure 5
Figure 5. a) Grey value coded range image of a detailed scan; b) segmented image with horizontal line; c) analysed image line with numbered tree objects (Bienert, 2013).

Figure 6 .
Figure 6.Principle of the determination of the tree position using edge points.

Figure 7 .
Figure 7. Detailed view of the separated point cloud with DTM points (blue dots, spacing 0.5 m) and voxels of the ground segment (green voxel, voxel size 0.2 m).

Figure 8 .
Figure 8. 2D range image after automatic removal of the ground points (range up to 30 m).