Individual Tree Crown Delineation Using Multi-Wavelength Titan LiDAR Data

The inability to detect the Emerald Ash Borer (EAB) at an early stage has led to the enumerable loss of different species of ash trees. Due to the increasing risk being posed by the EAB, a robust and accurate method is needed for identifying Individual Tree Crowns (ITCs) that are at a risk of being infected or are already diseased. This paper attempts to outline an ITC delineation method that employs airborne multi-spectral Light Detection and Ranging (LiDAR) to accurately delineate tree crowns. The raw LiDAR data were initially pre-processed to generate the Digital Surface Models (DSM) and Digital Elevation Models (DEM) using an iterative progressive TIN (Triangulated Irregular Network) densification method. The DSM and DEM were consequently used for Canopy Height Model (CHM) generation, from which the structural information pertaining to the size and shape of the tree crowns was obtained. The structural information along with the spectral information was used to segment ITCs using a region growing algorithm. The availability of the multi-spectral LiDAR data allows for delineation of crowns that have otherwise homogenous structural characteristics and hence cannot be isolated from the CHM alone. This study exploits the spectral data to derive initial approximations of individual tree tops and consequently grow those regions based on the spectral constraints of the individual trees.


INTRODUCTION
The inability to detect the Emerald Ash Borer (EAB) at an early stage has led to a colossal loss of different species of ash trees across Canada.The rapid spread of the EAB is attributed to the inability to detect the infected ash trees at an early stage and mitigate the risk of exposure to the neighbouring trees in the vicinity.The current environmental risks being posed by the EAB has led to a need for an accurate Individual Tree Crown (ITC) delineation method to identify tree crowns that are at a risk of being infected or are already diseased.In this study, an ITC delineation method was proposed that employed airborne multispectral Light Detection and Ranging (LiDAR) to accurately delineate tree crowns.
Multispectral LiDAR is being increasingly used for analysing large scale forest scenes and improving ITC delineation.Traditional methods of ITC delineation have primarily resorted to exploiting structural information about the crown to develop segmentation algorithms such as edge detection (Brandtberg et al., 1998), (Culvenor, 2002), (Koch et al., 2006), (Popescu et al., 2003), (Pouliot et al., 2005), region growing (Erikson, 2004) and watershed segmentation (Chen et al., 2006), (Schardt et al., 2002).Though successful results have been obtained, erroneous segmentation still occurs in areas of dense forest scenes due to the overlap of different tree crowns.With the availability of multispectral LiDAR data, research is now oriented towards exploiting spectral information along with structural information for improving the accuracy of ITC delineation.
In a dense forest scene, tree crowns have varied sizes and overlap among different crowns can form a tree cluster.Due to presence of tree clusters, ITC delineation from structural information alone, becomes a relatively difficult task.Conifers tree tops are * Corresponding author easy to detect whereas deciduous tree tops are relatively difficult to delineate.Deciduous tree tops have a flat crown structure which causes adjacent deciduous trees to have a significant overlap.Existing ITC delineation methods often employ a localized treetop to enhance the segmentation results.However, methods based on edge detection, region growing and watershed segmentation inevitably fail to recognize a boundary between two overlapping deciduous crowns and a tree cluster may hence be detected as a single crown.Additionally, the structural information does not entail a homogeneous signature for the entire crown as the boundaries of the crown tend to have different heights relative to the tree top.This consequently results in distortions and over segmentation of a single crown type, as well.
In recent years both LiDAR and passive multispectral optical imagery have been combined in data fusion techniques to enhance ITC delineation.However, majority of the studies have focused on deriving forest parameters from the multi-spectral optical imagery data rather than augmenting the ITC delineation itself (Zhen et al., 2016).The traditional methods of ITC delineation from integrated data have used the active imagery (i.e.LiDAR) for ITC delineation and incorporated the multispectral imagery at a later stage.By incorporating the spectral information with structural information in the form of multi-dimensional segmentation, the boundaries between different tree crown types can be identified to enhance existing ITC delineation techniques.Additionally, the spectral information captures a homogeneous region for a single crown type which prevents over segmentation of a single crown.The spectral information coupled with structural information can be employed to successfully resolve over segmentation of single crowns and under segmentation of overlapping crowns.The results can be cross validated and supplemented by classification techniques to identify erroneous segmentation.a In this study, the spectral information is used to improve the ITC delineation at three stages.During DEM generation phase, preliminary classification of spectral data is used to eliminate non-ground classes with high accuracy.The remaining point cloud is further refined using a TIN (Triangulated Irregular Network) densification technique.The improved accuracy of the DEM inevitably propagates to accurate tree height representations in the CHM.During pre-segmentation phase, the spectral data is used for accurate tree top identification for seed selection.The homogeneous tree structures in dense forest scenes make it relatively difficult to identify tree tops.Using the spectral data, the individual tree species within a tree cluster can be detected and multiple tree tops corresponding to multiple trees can be identified.In the segmentation phase the spectral data is used to augment the region growing algorithm in detecting boundaries within dense forest scenes with homogeneous tree structures.

STUDY AREA AND DATA
In this study the LiDAR data were obtained by Optech Titan multispectral LiDAR system which offers three spectral bands at wavelengths of 532 nm, 1064 nm and 1550 nm.The scan was performed over West Rogue, Scarborough area located southerneast of Toronto, Ontario, Canada.The raw LiDAR point data were provided in the form of three scans corresponding to the three wavelengths.The three scan strips were captured from different viewing angles.The raw LiDAR scan is illustrated in figure 1.

METHODOLOGY
As presented in figure 3, the LiDAR point cloud was initially processed to generate DEM, DSM and subsequently CHM.The intensity data was normalized and the spectral information was extracted.The combined spectral and structural information was used to identify tree tops and perform region growing segmentation.The steps are described in detail in the following subsections.

Intensity Normalization
The intensity data was normalized with respect to range using the following equations:

Digital Surface Model Generation
In this study, each of the three zones was processed to initially generate a gridded DSM.The DSM was generated by only accounting for the first returns and discarding subsequent returns.The resulting point cloud was then gridded, using Inverse Distance Weighted (IDW) interpolation, to later filter out the low frequency elevation trend (DEM) in the data.The grid size was determined by computing the average distance between pulses in the irregular point cloud using the following equation (Zhang, 2009): A grid size smaller or equal to the average pulse spacing was then selected to prevent any unintended interpolation.For the three zones the optimal grid size was determined to be 0.5 meters.The gridded DSM for the three datasets are illustrated in figure 4. and Zone 3 (Bottom -Forest scenes with high relief)

Digital Elevation Model Generation
The DEM was computed using preliminary classification results and an improved progressive TIN densification algorithm (Zhao et al., 2016 andAxelsson, 2000).Each zone was initially classified into three classes (Trees, Buildings and Terrain) using K-means clustering algorithm.The non-terrain features were subsequently removed and a TIN based model was generated from the remaining terrain points (Zhao et al., 2016).
Although the classification itself can be sufficient in removing majority of the non-ground points, it can inevitably lead to over erosion of the ground data.Due to over erosion of the ground data the interpolation can be inaccurate and hence densification of the ground points is required.The progressive TIN densification algorithm triangulates the classified ground points and iteratively introduces new ground points to the dataset.The overview of the process is illustrated in figure 5.
Figure 5. Flowchart for DEM generation The criteria for the addition of the point to the triangulation were measured with respect to the normal distance from the point to the façade of the TIN and the angles to each of three vertices of the TIN (Axelsson, 2000).The initial and final triangulation stages are captured in the figures 6 and 7, respectively.

Individual Tree Crown Delineation
The CHM was generated by eliminating the low frequency elevation trend (DEM) from the DSM.The structural information pertaining to the height and size of the tree crowns was used to determine seeds for the region growing algorithm.

Masking Non-Tree Pixels:
The structural information from the CHM and the spectral information was communally used to generate a non-tree mask.A Normalized Differenced Vegetation Index (NDVI) was generated, using the normalized intensity data in the visible and near infrared part of the spectrum, to isolate vegetation pixels.The tree height information from the CHM was then used to further filter out the non-tree pixels.A rule based approach was used to identify non-tree pixels: a pixel was labelled as non-tree if the NDVI and CHM values were below a threshold 0.35 and 4 meters respectively.The threshold values were determined empirically.

Tree Top Identification:
In this study, multi-scale method was adopted to account for different crown sizes in the forest scene (Schardt et al., 2002), (Gong et al., 2002).The multiscale approach was used to detect horizontal cross sections of the tree, as compared to tree tops, to initialize a large seed for seeded region growing segmentation.This approach also minimized the false tree top detection originating from complex deciduous tree structure.In addition to the structural information from the CHM, the spectral information was also exploited in the initial seed selection.The three spectral bands and the CHM were morphologically decomposed by an opening operation using a disk structuring element (SE).To account for different crown sizes, the decomposition was performed iteratively with different sized disk SEs.The regional maximum was then selected as the seed.Different scale spaces were used as estimates of the different horizontal cross sections of the tree crowns.The larger scale spaces were consequently omitted to prevent initialization of a large seed, as that could falsely identify a region larger than the actual crown size.

Region Growing:
With the seeds identified, region growing with neutrosophic logic was performed (Shan et al., 2008).The traditional methods of region growing rely on some measure of distance to merge pixels to their respective seeds, whereas neutrosophic logic introduces a degree of indeterminacy when evaluating the cost function to grow the seeds.The criterion for addition of a pixel to its respective seed depends on two quantities: degree of truth and level of indeterminacy.The degree of truth measures the normalized difference of the pixel value, of the pixel to be added to the seed, from the mean of the seed.The degree of indeterminacy introduces a measure of variance of the small circular neighbourhood surrounding the pixel.In this study, a small degree of indeterminacy implied that the region around the pixel was homogeneous hence representing a single tree.On the contrary a large indeterminacy value was indicative of a more heterogeneous region hence indicating presence of another tree crown type.For merging the pixels to their respective seeds, a rule based approach was used.If the degree of indeterminacy was smaller than a threshold value than the degree of truth of the individual pixel was calculated.A level of indeterminacy higher than the threshold value was indicative of the presence of another crown type and hence the degree of truth of a small circular region surrounding the pixel was calculated instead.The mathematical notation of the criteria is illustrated in expression 4 below (Shan et al., 2008).

RESULTS AND DISCUSSION
The initial seed selection, for the three zones, is illustrated in figure 9.The grown regions, with non-tree pixels masked out, are illustrated in figure 10.For zone 1, majority of the canopy was located in an urban area with predominantly isolated tree cluster types.Few cases of dense overlapping tree clusters were observed.Zones 2 and 3 had presence of dense forest scenes with homogenous morphology.
Based on visual observation tree top identification was relatively accurate in isolated and dense tree clusters, owing to the exploitation of the spectral bands.Existing ITC delineation techniques, based on active imagery, only examine structural information of the tree canopy to initialize a tree top.In areas of isolated tree clusters this approach is relatively accurate as the tree tops can be identified based on morphological filtering of the height models.In areas of dense forest scenes, the broad structure of deciduous trees can cause significant overlap between trees and tree top identification based on tree height models is relatively difficult.In the method proposed, the spectral information was used to identify different tree species within a dense forest scene with broad deciduous tree structure.The identification of individual tree species was then used to morphologically decompose the CHM and identify tree tops of trees with significant overlap.Such decomposition, based on CHM itself, only identified a single tree top for trees with significant overlap.Figure 11 illustrates the tree top identification based on the spectral constraints and compares it to tree top identification based tree height information alone.The identified tree tops were grown based on the spectral and structural constraints.The growth of the seed towards the boundary of the tree can introduce distortions as the height model varies from the tree top to the boundary of the tree.Based on visual observation of segmented trees in figure 10, it was noted that trees in dense tree clusters were delineated based on spectral characteristics instead of morphological constraints; hence the segmentation was uniform from tree top to the outer boundaries of the tree.Since spectral information preserves a homogeneous structure for the entirety of the tree, the distortions were minimized during the segmentation.Within a forest scene, the contrasting spectral signatures among trees were used to separate and identify individual trees within a dense structurally homogeneous tree cluster.In zones 2 and 3 the criteria for seed growth was predominantly based on the spectral constraints as the CHM represented a homogeneous region due to the overlap among trees.In each spectral band different tree species were identifiable and hence the regions were grown if the merge criteria were satisfied in multiple spectral bands.The growth of the seed was contingent upon the spectral and structural distance to the tree pixel (degree of truth) and a measure of spectral variance around the tree pixel (level of indeterminacy).The degree of indeterminacy introduces a measure of homogeneity around the tree pixel.This is indicative of presence of other tree species which are not identifiable based on a measure of distance alone.In the study the indeterminacy degree was used to identify whether the tree pixel belonged to the growing seed based on the neighbouring pixels around it.The level of indeterminacy was also a measure of location of the tree pixel within the tree.A more indeterminate (higher variance) region indicated a higher probability of the tree pixel being on the boundary of the tree.Criteria for merging the pixel was then based on a circular neighbourhood around the current tree pixel rather than the tree pixel itself.For more determinant (lower variance) region the tree pixel was relatively definite to be part of the tree represented by the growing seed and hence the merge criteria were based on that individual tree pixel alone.
Cases of over segmentation were also observed in zones 2 and 3. Initial seed selection based on the spectral constraints was imperative in identifying individual trees in a cluster, however there were cases where multiple trees were detected in a single tree crown.To mitigate this, different sized structuring elements were used to morphologically decompose the zones but a large sized seed would inevitably identify multiple tree crowns.Seed selection was therefore limited to only a range of structuring element sizes.This effect was prominent in large homogeneous tree clusters, in zone 2, that were not separable even in the spectral domain.An overall improvement in the accuracy of seed selection and the final segmentation was observed with the incorporation of spectral data.

CONCLUSION
Based on the multispectral LiDAR a new approach to ITC delineation was developed.The spectral data were used to improve the quality of the segmentation in three phases.Spectral data was exploited initially in the pre-processing phase to improve tree height representations in the CHM.Using constraints in the spectral and structural domain tree pixels were isolated and tree tops were identified.Multi-dimensional segmentation was performed to delineate tree crowns based on the unique spectral signatures of tree species.The method was effectively able to incorporate spectral data in isolating different tree species and hence augment the segmentation.The under segmentation of trees due to the homogeneity of the morphological information in the CHM was resolved by introducing spectral data in the segmentation procedure.The proposed method was seen as an improvement over single channel segmentation.
Treetop identification remains an important step to accurately delineate ITCs.Dense forest scenes with broad structured trees can then be identified with larger tree tops proportional to their structure to improve segmentation.Classification of different tree species based on spectral signatures can be used to supplement the existing segmentation techniques and validate the shape of the extracted tree crowns.In dense forest scenes classification of different tree species can further supplement and constrain the growth of the seed beyond the structural bounds of the tree.Work is being continued in improving segmentation and tree top identification from classification of different tree species. 6.

Figure. 1
Figure. 1 Titan Multispectral LiDAR Scanner DataThree zones were extracted from the raw LiDAR point cloud based on the variation of cover types and the elevation profiles.Zone 1 captures predominantly urban cover types with presence of isolated tree clusters.The terrain exhibits low relief and the elevation profile is primarily homogenous.Zone 2 captures a mixture of urban and forest scenes with high relief.The elevation profile is heterogeneous with higher elevations observed near the forest scenes.Zone 3 predominantly represents forest scenes with significantly high relief.The trees in zone 3 have significant overlap.Each zone measures around 170 by 145 meters.The location and proximity of the zones is illustrated on the satellite image shown in figure2.

Figure. 2
Figure. 2 The three study zones overlaid on the true colour composite of worldview image

Figure 3 .
Figure 3. Proposed Methodology for ITC Delineation from Multi-Spectral LiDAR Data 1) where   = The coordinates of the laser scanner  ,    ,   = The coordinates of a LiDAR point  = Flying height from GNSS ℎ = Height of the LiDAR point 2    =  * ( 2 ) (2)   where   = Normalized intensity  = Original intensity   = The maximal range where d = average distance between pulses  = laser pulse density in returns/m 2

Figure 4 .
Figure 4. DSM -Zone 1 (Top Left -Urban with low relief), Zone 2 (Top Right -Urban and Forest scenes with high relief)and Zone 3 (Bottom -Forest scenes with high relief)

Figure 8 .
Figure 8. DEM with 0.5-meter grid size -Zone 1 (Top Left), Zone 2 (Top Right) and Zone 3 (Bottom) 4) where   = Individual pixel  = Image domain   = The current seed (  ) = The indeterminacy value around pixel   (  ) = The degree of truth of pixel   ′(  ) = The degree of truth of a circular region around pixel

Fig. 11
Fig. 11 Boundaries of identified tree tops superimposed on Band 1 (Top Row -With spectral constraints, Bottom Row -Without spectral constraints)