Informatica 39 (2015)271-275 271
Detection of Ground in Point-clouds Generated from Stereo-pair Images
Domen Mongus and Borut Žalik
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
Keywords: digital terrain model, mathematical morphology, ©-mapping Received: June 24, 2014
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 ©-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 zračnih fotografij. Metoda uporablja podatkovno dekompozicijo na osnovi diferencialnih atributnih profilov in ©-mapiranja, s katerima doseže zaznavo najbolj kontrastnih povezanih komponent. Razpoznavo tock terena dosežemo z veckriterijskim pragovnim filtriranjem. Metode je evaluirana s primerjavo z digitalnim modelom reliefa, ustvarjenim iz podatkov LiDAR.
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 are 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]. Pingel 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 ©-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.
272 Informatica 39 (2015) 271-275
D. Mongus et al.
2 Theoretical background
Let g : E ^ R be a regular grid, where E c Z2 and p e E is a grid point. Consider a level-set E; c E given by the hight-level l as E; = {p | g[p] = l}. A connected component from E; 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 flat-zones) 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;} of g, each obtained by
T = {p | g[p] > l}.
(1)
YA(s)[p] = V{1 I P G Ck, A(Ck) > «},
(2)
AaA(g) = {yAL (g) - yA(g)},
(3)
g'[p] = v Aa(g)w,
g°[p] = A i I yA- (g)[p] - yA(g)[p] = g'[p],
(5)
where is infinum (i.e. the lower-bound). Consider a set of peak connected-components Cp = {C;p} containing a
ck
point p, i.e. Cp connected-component Cm, Af(g) is identified by
p e Cp. The most-contrasted with the respect to the given
max = \J l | ag°[p]-i < A(CP)
(6)
A peak connected component C;k e T; is defined by its height level l and its component-at-level index k . Let an attribute function A(C;k) that estimates a particular attribute of Ck, e.g. its area, diameter, or bounding-box. For simplicity, let A be increasing, thus, satisfying the condition Ck1 c Ck2 ^ A(Ck2) < A(Ck2). An attribute filter yA acting on g is at a particular point p defined by
where \J is supremum (i.e. the upper-bound). In other words, an attribute filter y At removes all the peak connected components not satisfying an attribute threshold condition a by assigning to each point p the maximal height-level at which it still belongs to a peak connected component Ck with A(Ck) > a. Since Vg, yA(g) < g, yA is an anti-extensive morphological filter named attribute opening. Its dual, an attribute closing ^f, is defined as
Yt(g) = -^t(-g).
A decomposition named DAP or differential attribute profile A has recently been proposed by Ouzounis et. al. [15]. A is based on progressive content reduction by filtering g at an increasing scale. Consider an ordered set of attribute thresholds a = {a4}, where i G [0,1] and ai-1 < aj, A is obtained by
where i e [1,1]. Thus, AA(g) is an I-long response vector registering the differences introduced by each particular yA, whilst yA (g) is a grid residual.
Recently, Mongus and Zalik [12] proposed ©-mapping that registers the most-contrasted connected-components and estimates their arbitrary attributes by observing characteristic values contained in AA. Namely, ©-mapping estimates the most-contrasted connected-components from g by registering the maximal responses from AA(s) and filtering scale at which they are obtained. Formally, ©(g,A, a) : g ^ (g',g°), is at p given by
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 Cm„x 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 ©-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 Rg is defined by the point-density Dl as Rg = 1.0/Dl. 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]
g[pn] =
E

g[pn] d
—r Pn
E.
PnEWpn aPn
(7)
(4)
where pn is a grid-point from the neighbourhood Wpn of p;, and dPn is the Euclidean distance between p*n and
Detection of Ground in Point-clouds.
Informatica 39 (2015)271-275 273
pn. 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 Wpn 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)
A(C P
\ max
gc[p] = 1-m
9n * DT(CPax)
(9)
where DT(Cmax) is a function that estimates the average distance of a grid-point contained within Cmax to the closest background point.
After g', g°, and gc are estimated, a set of ground grid-points G is recognized with a multicriterion threshold function given by
a structuring element. In our case, final DTM is obtained by
DTM [p]
g[p] ; g[p] - Yw (g) < R/2.0
Yw (g)[p] ;	otherwise
(11)
where w is box-shaped structuring element of size 5 x 5.
where i G [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 ©-mapping:
-	g' describes the height difference or residual of the most-contrasted connected-component from its background and is estimated by eq. 4,
-	g° describes the area of the most-contrasted connected-components according to eq. 5,
-	gc is a function describing shape-compactness of the most-contrasted connected-components and is estimated based on a well-known distance transformation as [13]
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. 1a). 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. 1b 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, ©-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.
g = {p | g'[p] < tr, g°[p] < ts, gc[p] < tc}, (10) 6 Acknowledgments
where tR, tS, and tC 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 yw, where w is
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".
274 Informatica 39 (2015) 271-275
D. Mongus et al.
Acknowledgement
Acknowledgement text.
_....... ..
k/y . :t* . .
A J
(a)
(c)
(d)
Figure 1: Estimation of DTM from (a and b) stereo-pair images and (c and d) the reference LiDAR data.
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. Legue-dois, 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 Photogram-metry 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. Bamler, N. Yague-Martinez, and H. Breit. Spaceborne spotlight SAR interfer-ometry 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.
Detection of Ground in Point-clouds.
Informatica 39 (2015)271-275 275
[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:1-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 andMachine 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.
[18]	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 ofSelected 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 dualdirectional slope filter. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 38(7B):628-632, 2010.