FILTERING AIRBORNE LIDAR DATA BY AN IMPROVED MORPHOLOGICAL METHOD BASED ON MULTI-GRADIENT ANALYSIS

The technology of airborne Light Detection And Ranging (LIDAR) is capable of acquiring dense and accurate 3D geospatial data. Although many related efforts have been made by a lot of researchers in the last few years, LIDAR data filtering is still a challenging task, especially for area with high relief or hybrid geographic features. In order to address the bare-ground extraction from LIDAR point clouds of complex landscapes, a novel morphological filtering algorithm is proposed based on multi-gradient analysis in terms of the characteristic of LIDAR data distribution in this paper. Firstly, point clouds are organized by an index mesh. Then, the multigradient of each point is calculated using the morphological method. And, objects are removed gradually by choosing some points to carry on an improved opening operation constrained by multi-gradient iteratively. 15 sample data provided by ISPRS Working Group III/3 are employed to test the filtering algorithm proposed. These sample data include those environments that may lead to filtering difficulty. Experimental results show that filtering algorithm proposed by this paper is of high adaptability to various scenes including urban and rural areas. Omission error, commission error and total error can be simultaneously controlled in a relatively small interval. This algorithm can efficiently remove object points while preserves ground points to a great degree.


INTRODUCTION
The technology of airborne Light Detection And Ranging (LIDAR) can acquire three-dimensional coordinates of laser footprints on reflective objects scanned beneath the flight path.The points may drop on buildings, vehicles, vegetation, ground and so on.Separating ground and non-ground points is generally the first step of any application of LIDAR data, which is called filtering.
LIDAR data filtering is a challenging task, especially for area with high relief or hybrid geographic features.Many related efforts have been made by a lot of researchers in the last few years.The filter based on iterative TIN densification chooses some lowest points of neighborhood to form an initial sparse TIN, and then added some points to the TIN if the points meet some criteria (Axelsson, 1999;2000;Sohn and Dowman, 2002).This algorithm does well in handling surfaces with discontinuities, but it is more likely influenced by low outliers.And what is more, it cannot remove man-made objects such as bridges.Building TIN is the process that largely increases computation load for the magnitude of LIDAR data.Reference (Pfeifer et al., 2001) proposed a filtering algorithm by hierarchical robust interpolation using data pyramids similar to the process of terrain surface approximation above.The kind of filters based on mathematical morphology attracts much attention considering simplicity, obvious physical meaning.Reference (Vosselman, 2000;Vosselman and Maas, 2001) organizes raw points in a Delaunay triangulation, and conduct filtering by the erosion operation used for mathematical grey scale morphology.Filter functions are derived by a representative training dataset.A point is classified as a ground point based on distances and the filter functions.The shortcoming of this design is that it is difficult to find the filter function and the representative training dataset suitable for all kinds of scenes.The filters proposed in (Roggero, 2001;Sithole, 2001) are both variants on the morphological filter developed in (Vosselman, 2000;Vosselman and Maas, 2001).Their works are creating other complex filter functions that work better.Reference (Zhang et al., 2003) develops a progressive morphological filter.The laser points are firstly interpolated to generate a regular grid.Then morphological opening operation is performed iteratively to remove object points by gradually increasing the filter window size and the elevation difference thresholds.The iteration is repeated until the size of the filter window is greater than the maximum building size.If the LIDAR data is of huge size, a very heavy workload of computation is a must.
In this paper, a novel filtering algorithm is proposed based on multi-gradient in terms of the characteristic of LIDAR data distribution.This paper is structured as 5 sections.Section 2 presents multi-gradient suitable for LIDAR filtering.The new filter is proposed in detail in section 3. The experiments are described and discussed in section 4. Section 5 gives some brief conclusions.

MULTI-GRADIENT OF POINT CLOUDS
The LIDAR data is a discrete three-dimensional point cloud that is stored according to the acquiring time.The distribution of raw points is irregular, and the data is huge.Therefore one highly effective method of data organization is in demand.In this paper, point clouds are divided by an index mesh which can support effective neighboring search as well as maintain the high resolution potential of raw data.Because the laser pulse of LIDAR is emitted according to a certain frequency, the average point spacing can be calculated, which can be the mesh interval.During any neighboring computation or morphologic operation, International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-1/W1, ISPRS Hannover Workshop 2013, 21 -24 May 2013, Hannover, Germany the corresponding mesh cell is firstly found out.Then all points in the cell are taken out to carry out computation or judging respectively.
In this paper, morphological gradient with symmetrical structuring element is adopted, which can overcome the influence of irregular distribution of LIDAR point clouds.
f denotes the data to be processed.D f is the definition domain of f. b is the structuring element whose size is (2w+1)×(2w+1).
Then the morphological operations applied to LIDAR are defined as follows: For the laser footprint p whose coordinates are (x, y, z), the elevation after dilation operation f⊕b is: (1) The elevation after erosion operation fΘb is: The multi-gradient of every LIDAR point can be divided into two kinds.One is G1 which is calculated by the formula (3).The other is G2 which is calculated by the formula (4).The size of dilation and erosion structuring element is 3×3 for gradient computation.

MORPHOLOGICAL FILTERING CONSTRAINED BY MULTI-GRADIENT
Low outliers are usually produced by system errors or multipath errors.These points always exist alone, and their elevation is far smaller than points in neighbourhood.Possible low outliers can be located by the gradient G1, and then are judged by analyzing whether the points are isolated.
Afterwards, non-ground points that are higher than ground are separated from ground points by an improved iterative opening operation constrained by multi-gradient.
Morphological opening operation is composed of erosion and dilation operations.The morphological opening operation of dataset f with structuring element b is defined as follows: fοb=(fΘb)⊕b (5) As the nature of opening operation, the objects that are smaller than structuring element are removed by opening operation.It makes filtering more efficiently to improve the traditional opening operation by multi-gradient, reducing the probability of rejecting ground points as object.
Firstly, the jumping points lying on ground like point b in figure 1(a) are chosen by multi-gradient to carry out erosion operation.
For an arbitrary LIDAR point p, if its G1 is higher than a certain threshold, and G2 is smaller than another threshold.Then, those points in the neighboring area around p are eroded with a structuring element of size (2w+1)×(2w+1).For instance, if the distance of a point n and point p is smaller than w, and point n is higher than point p, then the elevation of n after erosion is equal to the elevation of p.
Secondly, the points whose elevations become small after erosion are taken to carry on dilation according to formula (1) with the same size structuring element.
If the difference between the elevation of a point after dilation and its original is larger than a predefined threshold dh, the point is judged as object point.
The above opening operation is executed iteratively with both the size of structuring element and threshold increasing gradually until the structuring element is larger than the biggest building.The structuring element size w and the threshold dh are increased by an empirical formula (6).pointSpacing is the average point spacing of point clouds.w=4×dh/pointSpacing (6)

EXPERIMENTAL RESULT AND DISCUSSION
In order to test the filtering algorithm proposed, 15 sample data provided by ISPRS Working Group III/3 (Sithole and Vosselman, 2003;2004) are employed for experiments.These sample data include those environments that may lead to filtering difficulty as described in table 1.
The output of filtering algorithms can be numerically assessed by three types of error.S1 is the total number of ground points.S2 is the total number of object points.E1 is the count of ground points rejected as object.E2 is the count of object points accepted as ground.International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-1/W1, ISPRS Hannover Workshop 2013, 21 -24 May 2013, Hannover, Germany

CONCLUSION
In order to deal with the difficulties in LIDAR data filtering, a novel morphological filter is proposed in this paper.The jumping points lying ground are judging by multi-gradient, which are used to constrain the iterative morphological opening operation.The method is tested using the sample data which ISPRS provide for testing.The experimental results show that this algorithm is highly robust in all kinds of complex scenes, which efficiently removes object points while preserves ground points to a great degree.All types of error are controlled simultaneously in a relatively small interval.

Figure 1 .
Figure 1.Distribution sketch of LIDAR points.As shown in the figure 1(a), point a and point b are ground points, while point c are object point.Point a, b and c in the figure 1(b) are all ground points.G1 of point b in figure 1(a) and G1 of point b in figure 1(b) are maybe both higher than predefined threshold.But G2 of point b in figure 1(a) is small, while G2 of point b in figure 1(b) is large.So the jumping points lying on ground like point b in figure 1(a) can be identified according to multi-gradient, which are used to constrained the following morphological filtering.

Figure 2 -
4 are DSM (Digital Surface Model) generated by raw point clouds and DEM (Digital Elevation Model) generated after filtering of some sample data.

Table 1 .
Sample data for testing algorithm