PROGRESSIVE DENSIFICATION AND REGION GROWING METHODS FOR LIDAR DATA CLASSIFICATION

At present, airborne laser scanner systems are one of the most frequent methods used to obtain digital terrain elevation models. While having the advantage of direct measurement on the object, the point cloud obtained has the need for classification of their points according to its belonging to the ground. This need for classification of raw data has led to appearance of multiple filters focused LiDAR classification information. According this approach, this paper presents a classification method that combines LiDAR data segmentation techniques and progressive densification to carry out the location of the points belonging to the ground. The proposed methodology is tested on several datasets with different terrain characteristics and data availability. In all case, we analyze the advantages and disadvantages that have been obtained compared with the individual techniques application and, in a special way, the benefits derived from the integration of both classification techniques. In order to provide a more comprehensive quality control of the classification process, the obtained results have been compared with the derived from a manual procedure, which is used as reference classification. The results are also compared with other automatic classification methodologies included in some commercial software packages, highly contrasted by users for LiDAR data treatment.


INTRODUCTION
In recent years, there has been a significant increase of the widespread use of airborne laser scanning systems for obtaining digital elevation models with high-resolution spatial accuracy.
The key advantage of these systems is that the point measurement is made in a direct manner using a laser distanciometry system.However, these systems have the problem of not carrying out any interpretation of the scene, so the points are distributed on the ground according to certain pattern depending of the scanning instrument of the LiDAR, the flight and data capture configurations and the terrain characteristics.This implies, in the case of the applications where the desired product is the ground surface, the need for a subsequent filtering process for classifying data as belonging to the ground or belonging to the different objects located on it.

LiDAR data filtering methodologies for DTM generation.
The interest on LiDAR data filtering methodologies has supposed the apparition of different algorithms and methods in order to differentiate between those the points belonging to the ground, from those not belonging to it.These methods have different approaches and usually they are classified into five different groups.

Morphological filtering.
Using mathematical morphology operators, such as erosion and dilation filtering procedures implemented in LiDAR data, are characterized by the use of a structural element.These methods can be used for removing non-ground objects from the DSM data, such as, building and trees.Many authors have proposed different filtering methods based on this type of operators.Linderberger (1993) proposed the first filter based in the use of mathematical morphology applied to the filtering of profiles captured by a airborne laser profiler, using a opening process and considering a structural element with constant size and a certain threshold.Vosselman (2000) introduces an distance-dependent threshold (slope) with a constant size structural element.This modification implies a larger tolerance for more distant points.
According the same approach, Sithole (2001) and Roggero (2001) incorporated the local terrain slope at each point to the structural element definition instead of working with an average slope for the entire area (Vosselman filter).Thus, these new filtering methods have structural elements with constant size but variable form.Kobler (2007) proposed the structural element rotation in order to adapt it to the local slope.Considering that the use of a fixed-size structural element can provide not satisfactory results, several authors use variable sizes for the structural element, instead the fixed sized of the abovementioned methods.For example, Kilian et al. (1996) perform various opening processes with different window sizes.They use the window size as weight in the calculation of the final terrain height.Lohmann et al. (2000) proposed the Dual-Rank filter, which first uses a morphological filter in order to remove the non-ground points before carrying out the calculation of the final model.Zhang et al. (2007) proposed a progressive morphological filter based in the progressive enlargement of the window size.Using similar approach, Arefi and Hahn (2005) propose a progressive morphological filter using a geodetic distance operator.

Progressive densification filtering.
The filter methods included in this group work follow a progressive approach, that starts from a small number of points that generates a first surface approximation.In successive iterations, new points are adding to those that have previously been classified as belonging to the ground.This method was proposed by Axelsson (2000).The procedure begins with a triangulation process obtained from the lowest points presented in the area, using a grid with large dimensions spacing.The rest of the ground points are progressive included through a iterative process.This iterative process is based in the analysis of each point according to the triangle where the point is located, considering the distance from the point to the triangle and the angles formed between this point and the triangle vertices.
Within this filter group, other authors such as Hansen and Vogtle (1999) use the point height respect to its position in the corresponding triangle instead the distance used in other methods.Sohn and Dowman (2002) add a initial descending densification before the final ascending densification.

Surface based filtering.
Same to the algorithms based on the progressive densification, all methods classified into this group, use an initial surface reconstruction from a point cloud for a further filtering of the whole dataset.In progressive densification methods, the point assigned to the ground class are increased step by step, but the surface based methods typically start with a previous hypothesis that all point belong to the ground.Iteratively the influence of the non ground points will be reduced.
One of the most popular methods of this group is the method proposed by Kraus and Pfeifer (1998), known as robust interpolation method.This filter integrates data filtering and DTM interpolation in one single process.The purpose of this algorithm is to determinate which is the individual weight of each item in the modeled surface that represents the ground.Finally, the points are classified as ground points or not ground points, depending on whether or not exceeds a threshold value of different in height with respect to the final DTM surface.This method is improved in Pfeifer et al. (2001) and Briese at el. (2002).Moreover, there are proposals such as Elmqvist et al. (2001) which introduce the inner and outside strengths concepts to accomplish the ground surface location.Brovelli et al. (2004) proposed a method based on the splines surface calculation and edge extraction techniques for the ground points classification.

Clustering and segmentation based filtering.
This group deals with the localization of the homogeneous classes (clustering;Nardinocci et al. (2003), Jacobsen and Lohoman 2003), and Vosselman and Sithole (2005) or Filin and Pfeifer (2006).).This clustering can be done directly in the object space, using region growing techniques.Usually the homogeneity criterion is the normal vector or its variation, resulting flat surfaces in the first case, or variation surfaces in the second one (

Others.
Within this group, it is considered all those methodologies that have not fit in the above groups, for example, the repetitive interpolation filter proposed by Klober et al. (2007).
Last generation aerial fullwaveform laser systems present the capability of storage the full wave of each received pulse.The received waves represent the sum of the reflections of all the intercepted surfaces in the laser terrain footprint.The form of the received waves has been objective of different studies oriented to the improve of the object detection capability (Chauve et al., 2007;Lin et al., 2008).Also, it has been established the use of wave parameters (width and pulse amplitude, etc.).Such scan systems provide additional information of interest, although Lin and Mills (2009) show that these systems still require a significant research effort to demonstrate their potential in different applications.Doneus et al. (2008) and Wagner et al. (2008) use a modified robust interpolation filter method in order to include the information about the pulse width..But tests on different terrain types, particularly those giving raise problems with traditional systems (non full-waveform) are still under research.
Certainly, the LiDAR data filtering is a challenge issue.This interest has led to the emergence of a large number of different method, which show (and prove) the great effort made in this matter by the most important worldwide geomatics research centers.However, most authors conclude that with current methods is not possible to establish a fully automated filtering procedure for any data point configuration and scene.For this reason, it is interesting the integration of different methodologies in a unique filtering processing.

PROPOSED METHODOLOGY
We present a filtering method for LiDAR data classification according to an approach which combines the use of different filtering methods in one process.The main objective of this classification is to separate the points that are located on the terrain (ground points) and the points that are located on another objects (buildings, trees, etc.).The proposed methodology is based on the combined use of a segmentation process and a progressive triangulation densification.The proposed methodology is implemented on four stages: a preparation phase, a segmentation process and retrieval of existing segments using region growing techniques, a progressive triangulation densification oriented to the low height areas extraction, and, finally, a fusion process of the above results and the ground point classification and DTM generation.

Data preparation
The proposed methodology does not work directly using the original data.First, a process oriented to the detectionelimination of outliers (low-height data) presented in the area must be accomplished; next, a regular grid is generated.This grid must be a correct representation of the original data in cloud point structure.
For the regular grid generation, two considerations must be taken into account: the grid spacing and the assigned value to each grid position.The cell size (spacing) must be according to the average data density, in order to minimize the loss of spatial resolution and the presence of no-data cells.Once the grid spacing is defined, there will be cells where there are available several points and other cells that do not contain any point.Since the objective of this classification is to obtain the ground data, the lower point of each cell is selected, using only simple pulses and the last echoes of the multiple pulses.Since there are positions without any original data inside, it is necessary to perform an interpolation process (distance interpolation from neighboring positions).The whole process to obtain the regular grid is shown in Figure 1, and it is designed to obtain digital surface model of the scene in raster format (DSM).Middle: Initial regular grid generation -red pixels representing no data grids-; Lower: Final regular grid.

Segmentation process
The segmentation process aims with grouping cells in regions according to a similarity criterion.For this process, the different segments or regions presented in the scene must be generated, using a region growing process considering a binary space, where all the pixels will belong to a one of the two possible classes (ground or not-ground).
This methodological propose introduces the "planar point" concept, where a given point defines a local plane.For this consideration, it is necessary to obtain the plane define for each point and its neighbourhood.The accuracy of the adjustment is used to establish whether the point have a local behaviour as flat or not, generating for this purpose an image of the scene assigning to each position the precision adjustment centered at that position.For the generation of the binarized space from the previous image, it is necessary to introduce also the "barrier point" concept or point of geometric discontinuity.This concept refers to a point where there is an abrupt change in the geometric continuity defined at the scene of points corresponding to the outer edges and interiors of buildings, vegetation zones, bridges or other structures, etc., where precision adjustment is low and it should behave as a barrier in the later region growing process.
Once the threshold indicated to differentiate the "planar points" from the "barrier points", it is possible to generate the binarized image which will be the segmentation by region growing procedure.This region growing process ends with the image segmentation and it obtains the different segments present in the scene.The results of the process of segmentation can be seen in Figure 2.

Progressive densification process
Parallel with the segmentation process, another process oriented to the labeling of low points that belongs to the ground will be developed though a progressive triangulation densification according to a methodology similar that the proposed by Axelsson (2000).This process starts from an initial triangulation created from the low-points of the DSM.
The process will be a densification process of the triangulation based in the use a double threshold algorithm (angle and distance) that will incorporate this triangulation ground points with the required density in each case.In this way, it is possible to obtain a significant representation of the positions corresponding to the ground that could be a approximation of DTM in TIN structure.

Fusion of the results of both processes and final data classification.
Once obtained, the first segments in the scene (segmentation process) and on the other side, the first approximation of the DTM represented by the DSM positions labeled as ground points (triangulation process), it can be carried out a merging process of these results oriented to the initial DTM densification.
This process will consist of two distinct sub-processes.First, using triangulation points as seed points, being classified as ground segments all those segments that include any point of each triangulation.Next, a triangulation densification of the initial triangulation must be made.In this sub-process, we incorporate all those points of these segments that they are not previously incorporated.
With this process, we can obtain the final triangulation that represents the final DTM.Using this surface, the original data can be classified into two classes: ground data and non-ground data, analyzing the distance between the points and the generated surface.Figure 3 shows the densification process and the data classification.

RESULTS
To analyzing the results from the proposed methodology, the approach has been applied in an area of 550m x 550m size with a steep relief with important variations in height and complex buildings interspersed with cultivated area, as well as a road with several bridges over a river flowing along a ravine quite pronounced.As can be seen in table 1, which summarizes the data characteristics in this area, there are more than 700.000 points, with a 2.3 points/m 2 resolution and fitted with a 0.65m spacing regular grid.This is a highly complex urban area for the ground classification process, which is an ideal data set to test the efficiency of the proposed methodology.The analysis has been undertaken at punctual level, always taking as reference the classification made manually by an operator (TS-R), and comparing the obtained results in each case.From this analysis the coincidences are established (skill points), the "real ground points" classified as "non ground points" (omission error) and "real non-ground points" misclassified as "ground points" (commission error).
It can be seen clearly that the results obtained with the proposed methodology, although it presents more commission errors, the sum of both errors (omission and commission) is lower than those obtained using the commercial software.Additionally this method (CL-O) obtains more than 78% of success.

CONCLUSIONS
The use of progressive densification method for ground data classification of LiDAR cloud 3D points has a great efficiency.The final obtained results show an important dependence on the preset thresholds and the terrain characteristics of the scene.The automatic selection of the thresholds is complex considering the large variety of terrain types and characteristics of the input data.
Techniques of classification based on the use of growth regions are very fast and they have an easy implementation, although the results will depend on the availability of an efficient seeds search procedure.We propose a procedure based on the progressive triangulation densification for the terrain regions growing.
The combined use of both methods, progressive densification at a preliminary stage for the seed establishment and subsequent region growing is proposed.The approach provides goods results in an efficient manner in scenes with different characteristics, specially in urban and suburban areas.

Figure 1 .
Figure 1.Regular grid generation process.Upper: original data; Middle: Initial regular grid generation -red pixels representing no data grids-; Lower: Final regular grid.

Figure 2 .
Figure 2. Segmentation process.Upper left: MDS; Upper right: accuracy of local planes adjustment; Lower left: binary image; Lower right: segmentation.

Figure 3 .
Figure 3. Upper: Terrain segments generation process; Lower: left) Definitive triangulation densification; right) Final classification of ground data.
The obtained results are compared to a reference classification and several automatic classifications performed with commercial software.The reference classification has been made manually using TerraSolid (TerraScan module)classification CL-R-.Besides, some automatic classifications for comparison have been made with the same software (TerraScan) using different parameters configurations, considering the recommended values from the user's manual (classification TS-1, TS-2 and TS-3) (see Table2 Several configurations of the automatic classification process of TerraScan software.The obtained results of these classifications, and the results obtained from the proposed methodology application (classification CL-0) is shown in the figure 4.

Figure 4 .
Figure 4. Obtained results in the data classification for the different approaches.