DETECTION OF TREE CROWNS BASED ON RECLASSIFICATION USING AERIAL IMAGES AND LIDAR DATA

Tree detection using aerial sensors in early decades was focused by many researchers in different fields including Remote Sensing and Photogrammetry. This paper is intended to detect trees in complex city areas using aerial imagery and laser scanning data. Our methodology is a hierarchal unsupervised method consists of some primitive operations. This method could be divided into three sections, in which, first section uses aerial imagery and both second and third sections use laser scanners data. In the first section a vegetation cover mask is created in both sunny and shadowed areas. In the second section Rate of Slope Change (RSC) is used to eliminate grasses. In the third section a Digital Terrain Model (DTM) is obtained from LiDAR data. By using DTM and Digital Surface Model (DSM) we would get to Normalized Digital Surface Model (nDSM). Then objects which are lower than a specific height are eliminated. Now there are three result layers from three sections. At the end multiplication operation is used to get final result layer. This layer will be smoothed by morphological operations. The result layer is sent to WG III/4 to evaluate. The evaluation result shows that our method has a good rank in comparing to other participants’ methods in ISPRS WG III/4, when assessed in terms of 5 indices including area base completeness, area base correctness, object base completeness, object base correctness and boundary RMS. With regarding of being unsupervised and automatic, this method is improvable and could be integrate with other methods to get best results. * Email: eng.siamaktalebi@gmail.com Mobile Phone: +98 9149718713


INTRODUCTION
Detection and classification of objects on earth were and still are important fields for researchers in different majors including Remote Sensing and Photogrammetry (Rottensteiner et el., 2011).As emerging new sensors like laser scanners, developing Photogrammetry field, and utilizing digital cameras, methods of detection and classification are got into new era.High resolution and high spectral digital cameras have lead researchers to develop and introduce new indices to detect a variety of objects on earth.LiDAR data give 3D coordinates of points directly, that this ability makes it a simple task to differentiate between smooth and rough planes.Smooth planes usually designate man-made objects and rough planes mostly designate natural grounds.Because LiDAR is an active sensor it has no problem dealing with shadow areas, while shadow areas are challenging in aerial images.At high resolution aerial images boundary of objects like Buildings is clearly notable, but in LiDAR data there are some problems in detecting such boundaries.Considering advantages and disadvantages both LiDAR and aerial imagery it seems integrating these two data sources is the best option (Rottensteiner et el., 2008).Tree detection in complex city scenes because of existing high buildings is a more difficult task than tree detection in plains and cities with low buildings.There are lots of methods and algorithms in detection and classification field but it is not possible to compare those together.This is because of lack of bench mark data sets.In other hand most of algorithms and methods have tested in different data sets.To overcome this problem and making it easier to compare methods together a working group established in ISPRS, named WG III/4.This WG grants a bench mark data set to participants and encourages them to test their methods on this data and send the results to WG for evaluation (F.Rottensteiner et el., 2011).In continue some of related works about tree detection are mentioned.We separated WG III/4 participants' related works from the others.

Related works
Tucker introduced Normalized Difference Vegetation Index (NDVI) on the basis of plants reflection properties in Red (R) and Infrared (IR) bands.This index is used to detect vegetation cover (Tucker, 1979).Ono introduced a new Shadow Index (SI) in 2007 (Ono, 2007), that Grigillio in his article used this index to detect buildings (Grigillio et el., 2011).Cheuk-Yan in his article enhanced bands which were in shadow areas using Gomma Correction (GC) and Linear Correlation Correction (LCC).According on his article detection of trees after enhancement was much easier, because enhanced areas were spectrally more near to non-shadow areas (Cheuk-Yan Wan et el., 2012).Cai introduced Shadow Index based on Hue Saturation Intensity (HSI) colour space and used difference, ratio and normalized difference to detect shadow areas, simultaneously he used NDVI to improve the results (Dong Cai et el., 2010).Sarabandi introduced a new transformation for shadow detection based on Color Invariant Indices (CII) (Sarabandi et el., 2006).

A. Moussa, Canada (CAL):
This method has used the combination of ALS point cloud and orthophoto data.In this method objects are classified into building, tree and ground segments using a rule-based segmentation.Employing spectral data obtained from aerial images the segmentation result is improved.At last morphological operations are used to smooth labeled image (Moussa & El-Sheimy, 2012).
1.2.2J. Niemeyer, Germany (HAN): This method uses a supervised ALS point's classification on the basis of Conditional Random Fields (CRF), which employs a statistical model of context (Niemeyer et el., 2011)

INPUT DATA
Our input data in this article are IR, R and G bands from aerial images and ALS point clouds.IR, R and G bands have wave lengths as below:  IR covering from 0.79 to 89 µm,  R covering from 0.61 to 0.68 µm and  G covering from 0.50 to 0.59 µm.ALS point clouds have 4 kinds of data including: First Pulse Range (FPR), Last Pulse Range (LPR), First Pulse Intensity (FPI) and Last Pulse Intensity (LPI).In this article we would use FPR.

METHODOLOGY
Fig. 1 shows a flowchart for our methodology.This method can be divided into 3 sections.First section uses aerial images and both two other sections use ALS point clouds as input data.In the first section NDVI and SI indices are calculated using aerial image bands.Using a ratio between these indices we would obtain a new index, called ESI.This index shows vegetation in shadow areas with more values per pixel.This index and NDVI are reclassified into some classes using equal intervals.Then a linear equation is wrote between these indices, which this equation causes vegetation cover whether in sunny areas or in shadow areas get closer to top classes.At the end of this section a rule-based function is used to obtain the first section result layer.The problem with this section is the fact that the result layer includes grasses and bushes which are not desirable.In two other sections we would try to overcome these problems.In the second section a RSC is calculated using DSM.This index could eliminate entities which have a low value e.g.grasses.In the third section employing DSM and morphological operations we would obtain DTM, going ahead nDSM would be calculated using DSM and DTM.At the end of this section a rule-based function is used to reach the last result layer.The third section has ability to separate vegetation with low height, e.g.bushes, from higher vegetation, e.g.trees.
At the end we have three result layers from three sections.Using a point-wise multiplication of three last layers the total last result layer would obtained.Morphological operations could be used to smooth the result layer.In the following we would explain each section with more details. Where In this index brighter pixels show vegetation cover.For computing total slope Eq. 9 is used (Ritter, P. 1987 Where Z = threshold value

nDSM:
In this section our aim is to separate entities with a specific height from other entities.For this purpose DTM is calculated employing Geodesic Dilation, and then, using DTM and DSM we have extracted nDSM (H.Arefi, M. Hahn, 2005):

Third Section Result Layer:
At last using a threshold we would end up to the result layer: Where P = threshold value

Final Processing
Now we have three result layers from three sections.In this part using a point-wise multiplication the total last result layer is computed:

PREPROCESSING INPUT DATA
In this paper we have worked on all of three areas.ALS strip bands 9, 5 and 3 are respectively used for Area 1, Area 2 and Area 3. First pulse range data from ALS points are used to calculate DSM.With knowing the fact that point density in ALS point cloud is 4 points/m 2 , employed grid intervals for DSM is 25cm.These steps are done using ENVI.As there is no Ground Control Points (GCPs) provided, for geo-referencing aerial images the DSM of each area is employed.To reduce processing process and make aerial images pixel size equal to DSM pixel size, aerial images are resampled to 25cm pixel size using Nearest Neighborhood (NN) algorithm.Then areas' aerial images and DSMs are cropped.These steps are done using Georeferencing panel in ArcGIS.

EXPERIMENTAL RESULTS
This section is performed using MATLAB.To be ensuring parameters are chosen correctly they are obtained for Area 1 and then used for the other two areas.Firstly NDVI and SI are calculated and then NDVI is divided on SI.This gave us ESI.Both NDVI and ESI indices are scaled and transferred to [0,255].In continue after eliminating 2% of ESI and NDVI indices as outliers, they are reclassified into 25 classes by equal intervals.Now employing Eq. 6 and replacing X, A and B parameters respectively with 18, 0.3 and 0.7 LP layer is computed.Using 14 as a threshold in Eq. 7 the result layer in the first section is obtained.In the second section slope and RSC layers are computed using DSM and then employing 14 as a threshold in Eq. 12 the result layer for second section is computed.In the third section replacing h with 13m in Eq. 13 nDSM is calculated and then using 1m as a threshold in Eq. 14 the result layer in this section is obtained as well.Now we have three result layers from three sections.By multiplying these result layers the last result layer would be obtained.This layer is smoothed using morphological closing followed by morphological opening (Fig. 5).Structure element used in morphological operations is a disk with radius 1.

EVALUATION
Evaluation of the last results is performed by ISPRS WG III/4, so we are able to compare our results with other participants' results (Table 1).They used Fig. 6 as reference label images.This evaluation is performed on the basis of the method described in (Rutzinger et el., 2009).Our results details are reachable on the ISPRS website.In the website our method is introduced as TAF.Considering Table 1 it is obvious our method is not the best method, but it is at the acceptable ranking.If we want to compare this method on the basis of boundaries RMS, this method in the Area 1 is in the second place by 1.5m and in both Area 2 and Area 3 it is in the first place by 1.4m and 1.2m values respectively.This is because completeness and correctness indices for this method are close enough, which means this method has good precision comparing to other methods.On the basis of completeness on both per-pixel and per-object level this method is in the third place in Area 1, in the sixth place in Area 2 and in the third place in Area 3. In the basis of correctness on a per-pixel level this method is in the fourth place in Area 1 and Area 2 and in the fifth place in Area 3. Fig. 7 represents the results in a graphical way.In this figure we can see large trees are detected very well, but there are some problems in detecting small trees.
Because of using SI in this method trees which are in shadow areas are detected.Regarding to the results in Area 3 we can say this method has some problems in buildings boundaries.

CONCLUSION AND FUTURE WORK
Considering the results it is obvious indices used in this paper are good for the tree detection matter.This method is related to choosing precise values for some parameters, which are chosen manually.By changing these parameters the results could be improved, but it is a challenging job to select parameters manually.So Artificial Intelligence based methods, such as ANFIS, for determining parameters seem to be perfect.Replacing some parts of algorithm with alternatives could improve the results.For example, nDSM in this paper is calculated using Geodesic Dilation which it could calculate with other methods like Morphological Opening.Another example is method for calculating Shadow Index.We used just one method to calculate SI, which there are many methods to compute SI.
Figure .Representation of methodological workflow in brief 3.1 First section 3.1.1NDVI: Using IR and R bands NDVI is calculated: Figure .This figure shows the way of combining ReNDVI and ReESI 3.1.5First Section's Result Layer: At the end of this section a threshold used to get the last result: Our study area's data is gathered by ISPRS WG III/4.These data are captured over Vaihingen, Germany.ALS data has 4 points/m 2 point density and aerial images have 8cm ground pixel size and 11bits radiometric resolution.The field is divided into 3 areas.WG III/4 participants should test their methods for detection and extraction in these areas and submit their results to the committee to be evaluated.Each area could be seen fully or partially in some images.ALS point strips which cover three areas are shown inFig.4 (Franz Rottensteiner et el., 2011).

Figure .
Figure .This figure shows test areas aerial images (a) and ALS points' strips (b).

Figure 5 .
Figure 5.This figure shows the last result layer smoothed by morphological operations.White pixel show trees.a) Area1, b) Area2 and c) Area3

1.2.3 D. Grigillo and U. Kanjir, Slovenia (LJU): ALS
. At the end the result label image is smoothed by morphological operations.

1.2.4 Yao, Germany (TUM): In
this method both aerial images and range data are used in supervised classification.To obtain training data 10% of area is digitized manually.

1.2.5 Q. Zhan, China (WHU): Images
are ortho-rectified using DTM.ALS points, spectral and range data are used in a supervised classification. ): 3.2.2Rate of Slope Change:Zhilin Li, 2005 in his book wrote Eq. 17 to calculate RSC(Zhilin Li et el., 2005), but here we computed this layer in other way: