Int J Appl Earth Obs Geoinformation Detection of dead standing Eucalyptus camaldulensis without tree delineation for managing biodiversity in native Australian forest

In Australia, many birds and arboreal animals use hollows for shelters, but studies predict shortage of hollows in near future. Aged dead trees are more likely to contain hollows and therefore automated detection of them plays a substantial role in preserving biodiversity and consequently maintaining a resilient ecosystem. For this purpose full-waveform LiDAR data were acquired from a native Eucalypt forest in Southern Australia. The structure of the forest signi ﬁ cantly varies in terms of tree density, age and height. Additionally, Eucalyptus camaldulensis have multiple trunk splits making tree delineation very challenging. For that reason, this paper investigates automated detection of dead standing Eucalyptus camaldulensis without tree delineation. It also presents the new feature of the open source software DASOS, which extracts features for 3D object detection in voxelised FW LiDAR. A random forest classi ﬁ er, a weighted-distance KNN algorithm and a seed growth algorithm are used to create a 2D probabilistic ﬁ eld and to then predict potential positions of dead trees. It is shown that tree health assessment is possible without tree delineation but since it is a new research directions there are many im- provements to be made.


The importance of dead wood
The value of dead trees from a biodiversity management perspective is large. Once a tree dies, its woody structure remains for centuries and it contributes to forest regeneration while providing resources for numerous surrounding organisms (Franklin et al., 1987). More than 4000 species inhabit dead wood in Finland (Siitonen, 2001), where an estimate of 1000 species are threatened (Hanski, 2000). These species include animals, birds and other organisms, like fungi. Fungi contributes to wood decaying, formation of hollows and biodiversity, which supports the resilience of our ecosystem (Peterson et al., 1998).
In Australia, tree hollows play a significant role in managing biodiversity (Lindenmayer et al., 1997;Bennett et al., 1994). Nearly all arboreal mammals rely on hollows with the exception of the Koala (Phascolarctos cinereus) and perhaps Ringtail Possums (Pseudocheirus peregrinus) that preferentially make a stick nest. Additionally, numerous Australian bird species use hollows for shelters (Gibbons and Lindenmayer, 2002). Nevertheless, Australia has no real hollow creators unlike the northern hemisphere (e.g. Woodpeckers), and therefore it relies predominantly on natural processes of limb breakage, insect and fungal attack when access points are provided through damage caused by wind, storms and fire. This kind of hollows takes hundreds of years to form (Wormington and Lamb, 1999).
According to Gibbons et al. (2000), hollows are more likely to exist on dead trees or trees in poor physiological condition. In Australia, studies predict shortage of hollows for colonisation in the near future (Lindenmayer and Wood, 2010;Goldingay, 2009). A sample list of species that rely on hollows, provided by Forestry Corporation of NSW, is depicted at Fig. 1. Three of them are threatened (New South Wales Government, 2016). Consequently, automated detection of dead trees plays a substantial role in managing biodiversity.
As explained above, monitoring dead trees is essential for preserving a resilient ecosystem. Nevertheless, their distribution significantly varies making detection of them difficult (Kim et al., 2009). Remote sensing automates the process of monitoring forest and increases the spatial resolution of the monitored area.

Related work
Remote sensing was introduced for automating detection of dead trees since fieldwork is a time consuming task, considering the variance spread of trees and the spatial resolution of the area of interest. From a classification perceptive, the task of identifying dead standing and dead fallen trees is different. Fallen trees are identified by detecting segments or line-like features on the terrain surface using LiDAR (Polewski et al., 2015;Mcke et al., 2013). Regarding standing dead trees, their shape (reduced number of leaves or broken branches) (Yao et al., 2012) and light reflectance (less green light illuminated) (Pasher and King, 2009) are important factors for identifying them.
Previous work on dead standing trees detection performs single tree crown delineation before health assessment (Yao et al., 2012;Shendryk et al., 2016b). Tree crown delineation is usually done by detecting local maxima from the canopy height model (CHM) and then segmenting trees using the watershed algorithm (Popescu et al., 2003). Improvements has been achieved by introducing markers controlled watershed (Jing et al., 2012) and structural elements of tree crowns with different sizes (Hu et al., 2014). Additionally, Popescu and Zhao (2008) analyse the vertical distribution of the LiDAR points in conjunction with the local maximum filtering of CHM.
In the case of Eucalyptus camaldulensis, tree delineation is a challenge due to their irregular structure and multiple trunk splits. Local maxima filtering, used for tree detection, leads to over-segmentation because each tree trunk split forms a local maxima. Shendryk et al. (2016a) published a Eucalyptus delineation algorithm that performs segmentation from bottom to top; the trunks point cloud is separated from the leaves and individual trunks are identified before the segmentation. Nevertheless, the density resolution starts from 12 points/ m 2 and goes up to 36 points/m 2 around forested areas. For small research projects capturing this high resolution is acceptable, but for larger areas, the density of the emitted pulses is above the optimal resolution for a cost effective versus quality acquisition (Lovell et al., 2005). The investigated area of this paper is 95,196 ha and the entire area was scanned. The resolution of our acquired data has a minimum specification of four pulses per square meter. This is considered an optimal resolution in respect to the cost. Nevertheless, due to tree heights (up to 43 m) and acquired pulse density, there are not enough peak returns from trunks to enable bottom to up delineation. An example of the missing information about the trunks is depicted in Fig. 2.

Study area
The study area (Fig. 3) is a native River Red Gum (Eucalypt camaldulensis) forest of size 95,196 ha 2 in south-eastern Australia. The regeneration of the Eucalyptus camaldulensis is extremely dependant on floods and therefore, their distribution in respect to density, health and age is highly variance (Kerle, 2005). Additionally, the height of Eucalyptus camaldulensis reaches up to 40 m and their structural complexity is high with multiple trunk splits (Wilson and of N.S.W, 1995). Fig. 4 shows the structure of the forest and the shape variations of dead trees.

Acquired full-waveform (FW) LiDAR
FW LiDAR are acquired by RPS Australia East Pty Ltd from 900 m above ground level, using a Trimble AX60 Airborne LiDAR sensor. The wavelength of the emitted laser was 1062 nm, the maximum scan angle was 60°, and the pulse rate was 400 kHz. The acquisition was held from March 6th till March 31st, 2015. The collected LiDAR were delivered into 206 flightlines, of which 13 are cross runs for geometric correction. There is also a 30% of swath overlap. The footprint spacing of the pulses is 4.3 footprints per square meter.
Traditional ways of interpreting FW LiDAR, suggests extraction of a denser points cloud using Gaussian decomposition (Neuenschwander et al., 2009;Reitberger et al., 2008). This project implements on the open source software DASOS (Miltiadou et al., 2015), developed by the authors of this paper. Influenced by Persson et al. (2005), who used voxelisation to visualise the waveforms, DASOS uses voxelisation for both visualisations and classification. It further normalises the intensities so that equal pulse length exists inside each voxel, making intensities more meaningful. The literature is also moving towards voxelisation with promising results obtained by recent publications (Cao et al., 2016;Hancock et al., 2017). Generally, DASOS aims to ease manipulation of FW LiDAR in a volumetric representation. Its first two features (polygal mesh extraction and aligned metrics with hyperspectral imagery) was presented at the ISRSE conference in 2015. This paper presents a new feature implemented and how it is useful for assessing forest health.

Field data
The field data were collected in July 2015 during the winter season of Australia by Interpine Group Ltd and Forestry Corporation of NSW. There are 33 plots with radius 35.68 m and area 0.4 ha allocated randomly inside the study area. A total of 2386 trees were individually measured. Tree measurements include the geo-location, the trunk diameter at the breast height (1.3 m), species and health conditions (i.e. dead or alive). The geo-location of each tree was calculated at postprocessing using the recorded magnetic bearing from the centroid of the plot and the distance from the centroid. Here, it is worth mentioning that a single tree may be recorded as multiple trees if there is a trunk split bellow the breast height of 1.3 m. Furthermore, 91.59% of the • There is a high variance in the density of the forest.
The testing/training samples of small trees may contain information from nearby alive trees or ground.
• A tree may have dead branches but still be alive.
• Eucalyptus camaldulensis have irregular shapes and multiple trunk splits making tree delineation to require very dense acquired data.
• The pulse density of the acquired data does not allow bottom to top tree delineation. Additionally, crown detection from DEM leads to over-segmentation due to the multiple trunk splits.
• Dead trees are identifiable by their light reflectance, but the acquired data do not contain this kind of information. Therefore, the classifier depends predominantly on tree shape, which is not an independent factor of identifying dead trees; a tree may not have leaves but still be alive.
• If a tree has a trunk split below the 1.3 m height, then it is recorded as multiple trees; inconsistent meaning of the "one tree" term.
• Field data contain small trees, that are non detectable within the FW LiDAR.
• Unknown accuracy of the geo-locations of the trees; Some trees appear to be on the ground, once visualised on top of DEM: M. Miltiadou et al. Int J Appl Earth Obs Geoinformation 67 (2018) 135-147 trees are River Red Gum.
The field data contain information for 260 dead trees. Nevertheless, not all of them are considered useful for biodiversity. Dead trees with big Diameter at Breast Height (DBH) are more likely to contain hollows. Additionally, trees with DBH smaller than the footprint spacing of the LiDAR are not identifiable from the FW LiDAR. Table 1 shows the number of dead and alive trees in respect to their DBH. Please note that for training the classifier equal numbers of dead and alive trees have been used to reduce statistical bias.
Within the field data, some plots exist on two flightlines due to flights overlap. Overlaps exist at the edges of the flightlines and their scan angle significantly varies. For that reason, each unique set of field plot and flightline is considered as a test/training plot. This results into 50 plots used for training and cross-validation. Table 2 underlines the challenges faced while working on the detection of dead standing Eucalyptus camaldulensis. The challenges are categorised into three groups (the nature of the study area, the acquired data, the field data) and they all influence the quality of the classifier.

Methods and algorithms
This sections explains the methodology applied. In a few words, Random Forest is applied for identifying the most significant features. Then, a weighted KNN algorithm is used for generating a probabilistic field. Once ground pixels are removed, a seed growth algorithm segments potential dead trees and positions are finally assigned. More details, including intermediate steps for noise reduction, are explained below.
But before proceeding to the details, it worth highlighting the reasons of choosing those algorithms. As shown in Fig. 4, the shapes of the dead trees significantly vary from one to another. Therefore, during statistical analysis they may form multiple clusters with variant shapes and tangle between the clusters of the alive trees. For that reason, it is preferable to use a classification algorithm that supports non-linear boundaries and makes no assumptions of them. All parametric classifiers (e.g. Support Vector Machine and Bayesian) create models that summarise the properties of the class of interests, while the k-nearest neighbour (KNN) algorithm preserves all the information in some way. Nevertheless, KNN is prone to noisy parameters and for that reason the Random Forest Algorithm is used first for identifying the most significant features and reducing dimensionality.

Subtract DTM from FW LiDAR
For this paper, a new feature on DASOS was added for subtracting pre-calculated DTMs generated using the Quick Terrain Modeller. The Table 3 Explanation of a sample, the most relevant processed features, exported by DASOS.

No
Label Description 1 Height_Middle_Column The height of the middle column of the shape Height_Mean The Mean height of all the columns Height_Median The Median height of all the columns 1 Height_Std The Standard Deviation of the heights 2 Top_Patch_Len_Std The Standard Deviation of all the top patches 3 Dis_Std The Standard Deviation of the distances between the central voxel and every non-empty voxel 4 Per_Int_Above_Iso Percentage of voxels that contain an intensity above the isolevel 5 Top_Patch_Len_Mean The Mean length of all the top patches Top_Patch_Len_Median The Median length of all the top patches 7 Dis_Mean Mean distance from the central voxel to every nonempty voxel 8 Dis_Median Median distance from the central voxel to every non-empty voxel 9 Sum_Int_Diff_Z The Mirror Summed Difference of the intensities using the middle column in the z-axis as the axis of symmetry 10 Sum_Int_Diff_X The Mirror Summed Difference of the intensities using the middle column in the x-axis as the axis of symmetry Fig. 6. Information about the feature vectors created.
M. Miltiadou et al. Int J Appl Earth Obs Geoinformation 67 (2018) 135-147 subtraction of the DTM is done during the voxelisation; each terrain value is subtracted from the position of each waveform sample and not from the origin of the pulse since the terrain at the origin and the terrain at the position may vary. Fig. 5 shows an example of a DEM generated before and after the subtraction using DASOS.

The new feature of DASOS; feature vectors
This paper presents the new feature of DASOS, which is useful for characterising object inside the 3D space (e.g. trees). For each column of interest, inside the voxelised FW LiDAR, information from around its local area are exported as a feature vector. Multiple feature vectors are listed within.csv files for easy interpretation within statistical software packages. There are two types of exported information: processed and raw. The processed option lists information like distribution of nonempty voxels and standard deviation of heights (Table 3). The raw option lists the voxel intensity values of the local area. Additionally, there are two shapes of the local area from where the features are obtained (cuboid and cylinder). The size of each shape is user defined. Here, this new feature of DASOS is used for generating feature vectors used as a likelihood in the classifier. Fig. 6 depicts the divisions of the datasets and related feature vectors generated. As aforementioned, there are 50 plots (including overlaps). For cross-validation, these plots were randomly divided into 5 equal training datasets. Another one was created by merging three of them to check whether the increased training samples improves classification accuracy. The feature vectors generated for each field plot are divided into two categories (processed and raw intensities) and two sub-categories (cylinder and cuboid shape), resulting into four types of feature vectors per plot. For each type, three.csv files are generated: dead trees, alive trees and columns from unknown population. The first two are used for training and the last one for testing. Each feature vector represents the area within either a cuboid or cylinder. The dimensions of those shapes are slightly smaller than the average size of the dead trees to reduce noise. Additionally, the values of the vectors are normalised to belong in the range [0, 100].

Random Forest
Random Forest generates multiple regression trees by randomly sampling the data at its nodes and chooses the best predicting variables for each sampled data. The variable importance is defined by the influence it has to the classification once this variable is modified and the rest remain unchanged (Liaw and Wiener, 2002). This paper uses the R package for finding the most relevant features exported by DASOS (Section 4.2 in identifying dead trees).
It worth highlighting that Random Forest failed to find any relation between the feature vectors with the "Raw Intensities" due to the irregular shapes of Eucalyptus camaldulensis and the variant scan angle of each field plot. Nevertheless, "Raw Intensities" may be useful in classifying trees with smaller shape variance (e.g. pine trees).
Regarding the "Processed Intensities", Fig. 7 shows a list with the feature importance according to Random Forest and Table 3 gives the explanation of each important feature identified. The most important one is the standard deviation of heights. This is reasonable since the  canopy of dead trees has bigger height variance in comparison to alive trees whose canopy is leafy. Please note that the order of the significant features slightly vary according to the sub-dataset used.

Probabilistic field derived from weighted K-Nearest Neighbours Algorithm
Once the ten most significant features are identified, a weighted KNN is applied to generate a probabilistic field. As mentioned in Section 4.2, positive training feature vectors from dead trees and negative from alive trees has been created. To reduce bias, the same number of dead and alive trees is used at each training dataset.
Let's assume that T is a training dataset with n feature vectors: = … T x f x n N : ( , ( )), 1 .
The function f(x n ) ∈ {0, 1} return 0 for alive trees and 1 for dead trees. Therefore, the dataset T has this form: (2) Every feature vector t q ∈ T contains the 10 most important features, as identified by Random Forest (t = {t 1 , t 2 , …, t 10 }). Every feature is associated with a weight value according to its importance 1 0 ). Additionally: Let's define a data vector x = (x 1 , …, x 10 ) of an unknown population. How do we calculate the probability of vector x to belong to the dead trees population? At first, the weighted Euclidean distance from x to every t q ∈ T is calculated as follow: Then the k−nearest training samples are selected. Here, k = 7 was considered sensible, but testing different values of k could lead to the optimal k value. The nearest 7 indices of the training samples are selected as follow: 7 is a subset of the training samples T and contains the k-nearest indices to x, which could be positive, negative or both.
For each ∈ v V i a distance-weight u i is calculated: Finally, the probability of being a dead tree is calculated as follow: where the function δ(a, b) returns 1 if a is equal to b and 0 otherwise. For each column of the voxelised FW LiDAR, a testing data vector x is created and its probability of being dead is calculated. Fig. 8 shows the probability field of each column to belong to the dead trees population. The big circle indicates the limits of the field plot and the small circles the actual dead trees locations.

Filtering
As shown in Fig. 8, there is 'Salt and Pepper' noise. This occurs when no pulses pass through the corresponding column of the voxelised space or when the height of the shape used to calculate the corresponding feature vector is bigger than the canopy height (common at ground). It is removed using a median 3 × 3 filter which assigns to every empty pixel the median value of its non-empty neighbouring pixels (Fig. 9a). A smoothing filter is afterwards applied for further noise reduction (Fig. 9b).

Removing ground pixels
Removing the ground pixels is a trivial task because the DTM has already subtracted. A histogram of the heights, defined by the top nonempty voxels times voxel size, was generated. There are three welldefined classes: ground, trees and noise (Fig. 10b). Ground and noise are removed and this processed is illustrated in Fig. 10. 4.7. Dead and alive thresholding, filtering, segmentation and position assignment Fig. 11a show how the the dead trees are thresholded from the alive M. Miltiadou et al. Int J Appl Earth Obs Geoinformation 67 (2018) 135-147 trees. To reduce over-detection of dead trees, the outliers are filtered out using a 3 × 3 median kernel (Fig. 11b). Afterwards, the pixels are grouped using a seed growth algorithm (Algorithm 1, Fig. 12a). Each group corresponds to a predicted dead tree. Its predicted location is the average geo-spatial location of the pixels that belong to that tree (Fig. 12b).
Algorithm 1. Seed growth algorithm for grouping nearby pixels of dead trees 1: P⟵ all pixels classified as dead 2: s⟵ 0 3: while not reached the end of set P do 4: get next pixel p ∈ P that is not assigned to a segment 5: assign pixel p to segment s 6: find K(p 1 , …, p n ) such that K ⊆ P and every p i ∈ K is a neighbour of p 7: ∀ p i ∈ K, p ⟵ p i and repeat from line 5 8: all pixels of segment s has been labelled 9: s ⟵ s + 1

Evaluation
The results are cross-validated according to the predicted locations of the dead trees and their distance from the actual dead trees. Three different test were undertaken during evaluation: does the increased training samples improve dead tree detection? What shape (cylinder or cuboid) performs better? Are the predictions better than a random prediction? The random prediction was generated by uniformly distributing locations with density equal to the average density of the    Tables 4 and 5 shows the precision and recall percentage achieved using a cylindrical shape to extract features for the likelihood. The D1, D2, …, D5 corresponds to the five divided datasets in the cross-validation. The D_1_2_3 uses all the training samples from D1, D2 and D3 to train the classifier. As shown in the corresponding charts of precision ( Fig. 13 and recall Figure 14), the increased amount of training samples do not improve the prediction. The bigger the training set is, the shorter the distances to the k nearest neighbours should be. But the field data contain noise (Section 3), which compensates the value of the increased samples. A less noisy training dataset with trees of similar heights should improve the results.
The second comparison between using a cuboid and cylindrical shape to extract features. Tables 4 and 5 and Figs. 13 and 14 show the results of the prediction using a cylindrical shape to extract features, while Tables 6 and 7 and Figs. 15 and 16 show the results obtained using a cuboid shape. The average results are similar but the cuboid shape has a wider range of predicted results. Trees do not have corners and therefore the information retrieved with the cuboid shape are less meaningful. Nevertheless, the cuboid shape is slightly bigger and therefore collects better information from big trees but more noise from small trees. This justifies the wider range of good/bad results.
Finally, there is a comparison between the average results obtained and the evaluation of the random prediction are shown in Tables 8 and  9 and Figs. 17 and 18. During the evaluation a predicted tree is classified as TP if there is a dead tree within the given distance. As the minimum distance to a dead tree increases, the number of FP reduces. The precision of the random classifier significantly increases after the 6m distance from an actual dead tree, because its distribution is uniform. Additionally, the dimensions of the shapes used are 7.2 m diameter and width for cylinder and cuboid respectively. Therefore, evaluation above this distance is meaningless.
From the aforementioned figures, it is shown that the methodology proposed performs better than random. This indicates that forest health assessment is possible without tree delineation. Additionally, field surveying could be improved by better planning using the results of this systems. Of course, this is a new research direction and there are many improvements to be made (e.g. extracting features resistant to tree height variations).

Conclusions and future work
The importance of dead wood in our ecosystem is large and it is monitored for managing biodiversity. The study area of these projects is a native Australian forest with Eucalyptus camaldulensis, where shortage of hollows available for colonisation is predicted. Dead trees are more likely to be aged and contain hollows and therefore detecting them is essential for protecting their inhabitants.
This paper proposed a new direction for detecting dead standing Eucalyptus camaldulensis and presents the new feature of DASOS (extraction of feature vectors). Previous work on forest health assessment uses tree delineation but this leads to over-segmentation when applied to Eucalyptus camaldulensis due to their irregular shapes and multiple trunk splits. Additionally the density of the acquired LiDAR makes bottom to top delineation impossible since information about the trunks are missing. Therefore, this paper investigates the possibility of detecting dead trees from voxelised FW LiDAR without tree delineation.
Field data were provided by Forestry Corporation of NSW, Australia and Interpine Group Ltd, New Zealand. The GPS positions of the dead and alive trees within the plots are given but these data contain noise. Additionally, some small trees listed within the field data are not  M. Miltiadou et al. Int J Appl Earth Obs Geoinformation 67 (2018) 135-147 detectable from the FW LiDAR.
Regarding the methodology, feature vectors characterising dead and alive trees are generated. Then those feature vectors are run over the volume to generate a 2D image for each testing plot with the probability of a pixel to be from a dead tree or not. Salt and pepper noise removal and smoothing filters are applied. The ground is afterwards removed and a threshold is defined to separate pixels containing dead and alive trees. Then extra filtering is applied and a seed growth algorithm is used to label each segment. Each segment corresponds to a dead tree prediction and its estimated location is the average position of the pixels that belong to the segment.
The results have been cross-validated. Overall, there are three outcomes. The increase amount of training samples does not improve the results of the classification because of the noise within the field data. The feature vectors derived from cylindrical shape are more reliable because the range of the recall and precision is smaller. By the end, the most important outcome is that the results was clearly better than the random prediction, justifying that it is possible to identify dead trees without tree delineation.
Nevertheless, this is the first attempt to assess health forest without tree delineation and therefore many improvements could be made. Fore example: • Check and improve accuracy of the field data using visualisations of the FW LiDAR.
• Adjust features extracted so that they are not height dependant.
• Categorise the trees according to their size and derive a sets of training feature vectors according to their height.
• Once the seed growth algorithm is applied, check the size and shape of the segments and the possibility of dividing big segments into multiple dead trees.
• Add negative samples from the ground and the edges of alive trees, because the system is confused at the edges of the alive trees.      M. Miltiadou et al. Int J Appl Earth Obs Geoinformation 67 (2018) 135-147