# Detection of ground in point-clouds generated from stereo-pair images.

This paper proposes a new approach for constructing digital terrain models (DTM) from the point-clouds generated from airborne stereo-pair images. The method uses data decomposition based on the differential attribute profiles and [THETA]-mapping for the extraction of the most-contrasted connected-components. Their filtering is achieved based on multicriterion threshold function. The method is evaluated by comparing the output DTM with the reference Light Detection and Ranging data (LiDAR).Povzetek: V clanku predstavljamo novo metodo za konstrukcijo digitalnega modela reliefa iz oblaka tock, ki so generirani iz stereo-parov zracnih fotografij. Metoda uporablja podatkovno dekompozicijo na osnovi diferencialnih atributnih profilov in [THETA]-mapiranja, s katerima doseze zaznavo najbolj kontrastnih povezanih komponent. Razpoznavo tock terena dosezemo z veckriterijskim pragovnim filtriranjem. Metodeje evaluirana s primerjavo z digitalnim modelom reliefa, ustvarjenim iz podatkov LiDAR.

Keywords: digital terrain model, mathematical morphology, [THETA]-mapping

1 Introduction

Digital terrain models (DTMs) are essential part of various spatial analysis, geographic applications, and virtual reality systems [19, 6, 14]. In recent years, a considerable effort has been directed towards developing efficient approaches for accurate DTM generation.

When considering DTM generation from point-clouds, the most often used approaches can, according to the literature, be classified as slope-based, linear prediction-based, and morphological methods [20, 9]. Slope-based methods [18, 21] achieve point-filtering by comparing the gradients between neighbouring points. Consequentially, they have difficulties filtering points on step slopes and tend to smooth terrain undulations [20, 9]. Linear prediction-based methods, on the other hand, have difficulties filtering small and low objects as they rely on rough surface approximation to establish a liner prediction of the terrain [8, 2]. Actual filtering is usually achieved by observing the points' residuals from the predicted surface. Preservation of sharp terrain details (e.g. ridges) can, therefore, be exposed as another weakness [20, 9]. By applying operations of mathematical morphology [5, 11,4, 16], morphological filters proved to be fairly resistant to previously exposed drawbacks. However, they are severely dependent on the definition of the structuring element, as large objects (e.g. buildings) cannot be removed using a small structuring element, whilst large structuring element tends to flatten terrain details (e.g. mountain peaks) [20, 9, 5]. Several attempts have been proposed for optimal definition of a structuring element, the most efficient of which arc based on a multi-scale filtering. A set of filters of different scales is used for this purpose and different threshold values are usually defined for each of them. A progressive filtering was proposed by Chen et al. [5], where thresholding is applied on height differences achieved by each filter. On the other hand, Mongus and Zalik [11] proposed data-filtering by iterating thin-plate splines towards the ground, where resolution is increased at each iteration by including points, filtered according to their residuals from the previously estimated surface. This, so-called, hierarchical multiresolution filtering has recently been improved by Chen et al. [4], Pingei et al. [16] have, on the other hand, based their approach on the slope estimation achieved by linearly increasing filtering scale. Since all of these methods are adopted for processing high-resolution point-clouds containing vast amounts of points (e.g. LiDAR data), iterative approaches may not always be appropriate. Mongus and Zalik have [12] proposed an efficient multiscale approach that avoids iterations by using attribute filters based on the max-tree data structure. Although the method proves efficient when filtering LiDAR data, its accuracy is not guaranteed when filtering low-resolution point-clouds (such as those generated from stereo-pair images) since it is based on the standard deviation of point heights.

This paper presents a new method for estimation of digital terrain model from point-clouds generated from airborne stereo image pairs. By considering [THETA]-mapping, the proposed method is an extension of [12], where a different set of attributes is used for the filtering. Section 2 explains theoretical foundation of connected operators from mathematical morphology that allows their efficient estimation. The method is explained in Section 3. Section 4 gives the results, whilst Section 5 concludes the paper.

2 Theoretical background

Let g : E [right arrow] R be a regular grid, where E [subset] [Z.sup.2] and p [member of] E is a grid point. Consider a level-set [E.sub.l] [subset] E given by the hight-level l as [E.sub.l] = {p | g[p] = l). A connected component from [E.sub.l] is named a flat-zone of g. A filter that acts on flat-zones rather than individual grid-points is named a connected operator [17]. A connected operator can eider remove a flat-zone (by merging it with some other flatzones) or leave it perfectly preserved, but it cannot brake it. If the decision about which flat-zones to merge is based on some of their attributes, this type of operator is named an attribute filter [1], Consider a set of all thresholded sets T = {[T.sub.l]} of g, each obtained by

[T.sub.l] = {p | g[p] [greater than or equal to] l}. (1)

A peak connected component [C.sup.k.sub.l] [member of] [T.sub.l] is defined by its height level l and its component-at-level index k. Let an attribute function A([C.sup.k.sub.l]) that estimates a particular attribute of [C.sup.k.sub.l], e.g. its area, diameter, or bounding-box. For simplicity, let A be increasing, thus, satisfying the condition [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. An attribute filter [[gamma].sup.A.sub.a] acting on g is at a particular point p defined by

[[gamma].sup.A.sub.a](s)[p] = [disjunction] {l|p [member of] [C.sup.k.sub.l] A([C.sup.k.sub.l]) [greater than or equal to] a}, (2)

where [disjunction] is supremum (i.e. the upper-bound). In other words, an attribute filter [[gamma].sup.A.sub.a] removes all the peak connected components not satisfying an attribute threshold condition a by assigning to each point p the maximal heightlevel at which it still belongs to a peak connected component [C.sup.k.sub.l] with A([C.sup.k.sub.l]) [greater than or equal to] a. Since [for all]g, [[gamma].sup.A.sub.a] (g) [less than or equal to] g, [[gamma].sup.A.sub.a] is an anti-extensive morphological filter named attribute opening. Its dual, an attribute closing [[phi].sup.A.sub.[alpha]], is defined as [[gamma].sup.A.sub.a] (9) = -[[gamma].sup.A.sub.a] (-g).

A decomposition named DAP or differential attribute profile [DELTA] has recently been proposed by Ouzounis et. al. [15]. [DELTA] is based on progressive content reduction by filtering g at an increasing scale. Consider an ordered set of attribute thresholds a = {[a.sub.i]}, where i [member of] [0,I] and [a.sub.i-1] < [a.sub.i], [DELTA] is obtained by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (3)

where i [member of] [1,7]. Thus, [[DELTA].sup.A.sub.a] (g) is an I--long response vector registering the differences introduced by each particular [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a grid residual.

Recently, Mongus and Zalik [12] proposed [THETA]-mapping that registers the most-contrasted connected-components and estimates their arbitrary attributes by observing characteristic values contained in [[DELTA].sup.A.sub.a]. Namely, [THETA]-mapping estimates the most-contrasted connected-components from g by registering the maximal responses from [[DELTA].sup.A.sub.a](s) and filtering scale at which they are obtained. Formally, [THETA](g, A, a) : g [right arrow] (g', [g.sup.o]), is at p given by

g'[p] = [disjunction] [[DELTA].sup.A.sub.a] (g)[p], (4)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (5)

where [disjunction] is infinum (i.e. the lower-bound). Consider a set of peak connected-components [C.sup.p] = {[C.sup.p.sub.l]} containing a point p, i.e. [C.sup.p.sub.l] = [C.sup.k.sub.l] | p [member of] [C.sup.k.sub.l]. The most-contrasted connected-component [C.sup.k.sub.max] with the respect to the given [[DELTA].sup.A.sub.a] (g) is identified by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (6)

where max defines the height-level of the most-contrasted connected-component. Note that possibly no response was obtained at a given p, meaning that the corresponding peak connected-components are not in contrast against their surroundings and, therefore, belong to the grid residual, i.e. background. In any case, an arbitrary attribute of [C.sup.p.sub.max] can then be measured and used as an attribute in multicriterion threshold definition.

3 Ground extraction from point-clouds

The proposed method generates a digital terrain models from point-clouds obtained by stereo-pair images in the following tree steps:

--Initialization is the first step of the method were input point-cloud is sampled into a grid,

--Point filtering is performed in the space of the most-contrasted connected-components obtained by [THETA]-mapping, and

--Construction of DTM is the final step of the method, where removed points are interpolated.

Each of these steps is discussed in continuation.

3.1 Initialization

In order to apply morphological operators on point-clouds, points are firstly sub-sampled into a regular grid g. The resolution of the grid [R.sub.g] is defined by the point-density [D.sub.L] as [R.sub.g] = 1.0/[D.sub.L]. When a particular grid-cell contains more than one point, the hight level of the grid-point is defined by the lowest one since it has the highest probability of being a ground point. On the other hand, interpolation is used in order to estimate the hight levels of the undefined grid-points g[p*] = UNDEF, obtained when there are no points contained within the corresponding grid-cells. In our case, the height level at p* is estimated using inverse distance weighting (IDW) as [10]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

where [p.sub.n] is a grid-point from the neighbourhood [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] of [p.sup.*.sub.n], and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the Euclidean distance between [p.sup.*.sub.n] and [p.sub.n]. Parameter r defines the smoothness of the interpolation. According to the evaluation of the spatial interpolation methods described in [3], accurate results are obtained when [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] contains at least three closest points and r = 2.

3.2 Ground filtering

In order to achieve extraction of the most-contrasted connected-components, the underlying definition of DAPs is given first. In compliance with demanded increasing property of the attribute used for grid decomposition, the proposed method constructs DAPs according to the area of the contained peak connected components A. An area threshold vector a is given as

a={20.0 * i}, (8)

where i [member of] [0,1], Note that a is given in square-metres, thus, its definition should be adjusted when the input point-cloud is not georeferenced. In any case, the following attributes of the most-contrasted connected-components are estimated by [THETA]-mapping:

--[g.sup.I] describes the height difference or residual of the most-contrasted connected-component from its background and is estimated by eq. 4,

--[g.sup.o] describes the area of the most-contrasted connected-components according to eq. 5,

--[g.sup.c] is a function describing shape-compactness of the most-contrasted connected-components and is estimated based on a well-known distance transformation as [13]

[g.sup.c][p] = A([C.sup.p.sub.max]) / 9[pi] * [bar.DT] ([C.sup.p.sub.max]), (9)

where [bar.DT]([C.sup.p.sub.max]) is a function that estimates the average distance of a grid-point contained within ([C.sup.p.sub.max]) to the closest background point.

After g', [g.sup.o], and [g.sup.c] are estimated, a set of ground gridpoints G is recognized with a multicriterion threshold function given by

G = {p | g' [p] [less than or equal to] [t.sup.R], [g.sup.o] [p] [less than or equal to] [t.sup.s], [g.sup.c][p] [less than or equal to] [t.sup.c]}, (10)

where [t.sup.R], [t.sup.S], and [t.sup.C] are residual, size, and compactness thresholds, respectively.

3.3 DTM construction

In the final step of the method, DTM is contracted by interpolating the heights of the non-ground points NG = E/G using IDW, as given by eq. 7. However, using r = 2 may not always be appropriate as it may produce some sharp unnatural terrain features. Additional smoothing is, therefore, performed based on morphological opening [[gamma].sub.w], where w is a structuring element. In our case, final DTM is obtained by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)

where w is box-shaped structuring element of size 5x5.

4 Results

In order to evaluate the method, a point-cloud has been generated from georeferenced stereo-pair image as proposed in [7] with approximately 17.000 points. Average point spacing was below 3.1m and average absolute height error was 5.3m in comparison to the reference data (see Fig. 1 a). The reference data was acquired with LiDAR technology. The referenced point-cloud contained more than 1.6 millions of points with average point-spacing below than 0.25m and average absolute height error below 0.1m (see Fig. 1c).

The reference DTM was obtained with [12] and was used for the evaluation of the proposed method (see Figs, lb and d). The results show that the proposed method is capable of removing important portion of noise as the average absolute difference of DTMs was lower than the average error of the point-clouds. Namely, the error is reduced to 4.8m. However, significant portion of DTM's details is missing due to the lower point-cloud resolution.

5 Conclusion

The paper proposes a new method for estimation of DTMs from point-clouds, generated by stereo-pair areal images. The method determines non-ground regions by estimating their geometrical characteristics, namely their sizes, shape compactness, and height differences from the background. As confirmed by the results, [THETA]-mapping provides sufficient solution for this purpose as great majority of errors were introduced by interpolation and lower data accuracy in comparison to LiDAR data.

6 Acknowledgments

This work was supported by the Slovenian Research Agency under grants L2 - 3650 and P2 - 0041. This paper was produced within the framework of the operation entitled "Centre of Open innovation and Research UM". The operation is co-funded by the European Regional Development Fund and conducted within the framework of the Operational Programme for Strengthening Regional Development Potentials for the period 2007-2013, development priority 1 : "Competitiveness of companies and research excellence", priority axis 1.1: "Encouraging competitive potential of enterprises and research excellence".

Domen Mongus and Borut Zalik

University of Maribor, Faculty of Electrical Engineering and Computer Science

E-mail: domen.mongus@um.si and borut.zalik@um.si, http://gemma.uni-mb.si

Received: June 24, 2014

Acknowledgement

Acknowledgement text.

References

[1] E. Breen and R. Jones. Attribute openings, thinnings and granulometries. Computer Vision and Image Understanding, 64(3):377-389, 1996.

[2] M. A. Brovelli, M. Cannata, and U. M. Longoni. LiDAR data filtering and DTM interpolation within GRASS. Transactions in GIS, 8(2): 155-174, 2004.

[3] V. Chaplot, F. Darboux, H. Bourennane, S. Leguedois, N. Silvera, and K. Phachomphon. Accuracy of interpolation techniques for the derivation of digital elevation models in relation to landform types and data density. Geomorphology, 77 (1-2): 126-141, 2006.

[4] C. Chen, Y. Li, W. Li, and H. Dai. A multiresolution hierarchical classification algorithm for filtering airborne LiDAR data. ISPRS Journal of Photogrammetry and Remote Sensing, 82:1-9, 2013.

[5] Q. Chen, P. Gong, D. Baldocchi, and G. Xie. Filtering airborne laser scanning data with morphological methods. Photogrammetric Engineering & Remote Sensing, 73(2): 175-185, 2007.

[6] R. Dinuls, G. Erins, A. Lorencs, I. Mednieks, and J. Sinica-Sinavskis. Tree species identification in mixed baltic forest using LiDAR and multispectral data. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 5(2):594-603, 2012.

[7] M. Eineder, N. Adam, R. Barnier, N. Yague-Martinez, and H. Breit. Spaceborne spotlight SAR interferometry with TerraSAR-X. IEEE Transactions on Geoscience and Remote Sensing, 47 (5): 1524-1535, 2009.

[8] H. S. Lee and N. Younan. DTM extraction of LiDAR returns via adaptive processing. IEEE Transactions on Geoscience and Remote Sensing, 41(9):2063-2069, 2003.

[9] X. Liu. Airborne LiDAR for DEM generation: some critical issues. Progress in Physical Geography, 32(1):31-49, 2008.

[10] C. D. Lloyd. Local Models for Spatial Analysis (2nd ed.). CRC Press, 2010.

[11] D. Mongus and B. Zalik. Parameter-free ground filtering of LiDAR data for automatic DTM generation. ISPRS Journal of Photogrammetry and Remote Sensing, 66 (1): 1-12, 2012.

[12] D. Mongus and B. Zalik. Computationally efficient method for the generation of a digital terrain model from airborne LiDAR data using connected operators. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, In press:l-12, 2013.

[13] R. S. Montero and E. Bribiesca. State of the art of compactness and circularity measures. International Mathematical Forum, 4 (25-28): 1305-1335, 2009.

[14] A. O. Onojeghuo and G. A. Blackburn. Characterising reedbeds using LiDAR data: Potential and limitations. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, (In press): 1-7, 2012.

[15] G. K. Ouzounis, M. Pesaresi, and P. Soille. Differential area profiles: Decomposition properties and efficient computation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(8): 1533-1548, 2012.

[16] T. J. Pingel, K. C. Clarke, and W. A. McBride. An improved simple morphological filter for the terrain classification of airborne LIDAR data. ISPRS Journal of Photogrammetry and Remote Sensing, 77:21-30, 2013.

[17] P. Salembier and M. H. Wilkinson. Connected operators: A review of region-based morphological image processing techniques. IEEE Signal Processing Magazine, 136 (6): 136-157, 2009.

118] J. Shan and A. Sampath. Urban DEM generation from raw LiDAR data: A labeling algorithm and its performance. Photogrammetric Engineering & Remote Sensing, 71(2):217-222, 2005.

[19] B. Sirmacek, H. Taubenbock, P. Reinartz, and M. Ehlers. Performance evaluation for 3-D city model generation of six different DSMs from air- and spaceborne sensors. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 5(1):59-70, 2012.

[20] G. Sithole and G. Vosselman. Experimental comparison of filter algorithms for bare earth extraction from airborne laser scanning point clouds. ISPRS Journal of Photogrammetry and Remote Sensing, 59(1-2):85-101, 2004.

[21] C. Wang and Y. Tseng. DEM generation from airborne LiDAR data by an adaptive dual-directional slope filter. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 38(7B):628-632, 2010.

Printer friendly Cite/link Email Feedback | |

Author: | Mongus, Domen; Zalik, Borut |
---|---|

Publication: | Informatica |

Article Type: | Report |

Date: | Sep 1, 2015 |

Words: | 3046 |

Previous Article: | Cervix cancer spatial modelling for brachytherapy applicator analysis. |

Next Article: | IJCAI 2015--the worst best ever. |

Topics: |