DETECTION OF PLANAR POINTS FOR BUILDING EXTRACTION FROM LIDAR DATA BASED ON DIFFERENTIAL MORPHOLOGICAL AND ATTRIBUTE PROFILES

This paper considers a new method for building-extraction from LiDAR data. This method uses multi-scale levelling schema or MSLS-segmentation based on differential morphological profiles for removing non-building points from LiDAR data during the data denoising step. A new morphological algorithm is proposed for the detection of flat regions and obtaining a set of building-candidates. This binarisation step is made by using differential attribute profiles based on the sum of the second-order morphological gradients. Any distinction between flat and rough surfaces is achieved by area-opening, as applied within each attribute-zone. Thus, the detection of the flat regions is essentially based on the average gradient contained within a region, whilst avoiding subtractive filtering rule. Finally, the shapes of the flat-regions are considered during the building-recognition step. A binary shape-compactness attribute opening is used for this purpose. The efficiency of the proposed method was demonstrated on three test LiDAR datasets containing buildings of different sizes, shapes, and structures. As shown by the experiments, the average quality of the buildings-extraction was more than 95%, with 96% correctness, and 98% completeness. In terms of quality, this method is comparable with TerraScan R ©, but both methods significantly differ when comparing correctness and completeness of the results.


INTRODUCTION
Recent advances in airborne Light Detection and Ranging (Li-DAR) technology have led to the development of accurate, reliable, and fast data-acquisition systems that go hand-in-hand with the increasing demand for generating accurate building representations (Lukač et al., 2013).Over the past few years, a considerable number of methods have been developed for this purpose that either rely on LiDAR data alone, or fuse it with aerial images.Although the latter contributes to the accuracy of buildingextraction by providing supplementary information beyond the scope of LiDAR systems (e.g.surface reflectance on multiple electromagnetic spectrum bands), they may not always be available and several drawbacks may be related to the noise that is present within these images, such as shadows, clouds, and highrise buildings (Meng et al., 2009).In any case, the processing of massive amounts of geometric data, with a lack of topology as obtained by LiDAR systems, is still a challenge.
Traditionally, a normalised digital surface model (nDSM) is computed for this purpose by subtracting the digital terrain model (DTM) from the LiDAR point-cloud (Lohmann et al., 2000).Over the past few years, several advances in building-extraction have been made that are based on mathematical morphology.Most of these methods apply morphological operations on a 2.5D grid generated from a LiDAR point-cloud in order to cope with the lack of topology.One early example was based on maximum filtering (Vosselman, 2000), and dual-rank morphological filtering (Lohmann et al., 2000).Tarsha-Kurdi et al. (2006) used directional gradient filters on interpolated DSM for detecting offterrain segment edges.The morphological opening was therefore applied in order to remove any pronounced segments remaining on the ground, and morphological closing for filling the gaps within non-ground segments.By using the least-square method, the authors detected non-ground planes that belonged to buildings.Zhang et al. (2006) used progressive morphological filtering (Zhang et al., 2003) to find non-ground points by gradually increasing the filtering scale, and thresholding the height differences.Afterwards, the authors used region growing and planefitting to detect buildings' regions.Vu et al. (2009) considered LiDAR data within a multi-scale morphological space.The authors analysed elevation clusters' features across the scale-space in order to detect buildings.Meng et al. (2009) proposed groundfiltering based on the elevation differences between neighbouring pixels, and then removing the majority of non-building points using morphological operators.The remaining non-building regions were removed by area and compactness analysis.Chen et al. (2012) extended progressive morphological filtering for detecting non-ground points, where the authors applied a regional growing and adaptive random sample consensus (RANSAC) algorithm for the detection of buildings' segments.They considered distance, standard deviation, and normal vectors' attributes.Recently, Cheng et al. (2013) have proposed a building extraction method using the reverse iterative mathematical morphological (RIMM) algorithm.This method is an extension of (Zhang et al., 2003), as it avoids relying on a constant slope by dynamically calculating those threshold values applied regarding height differences.Although this method overcomes the problem of dealing with buildings of different sizes, it still uses predefined structuring elements and is, therefore, heavily subjected to buildings' shapes.
This paper proposes a morphological approach for building extraction from LiDAR data.The method is independent of the sizes as well as the shapes of buildings by applying differential attribute profiles, whilst traditional morphological operators are used for removing non-building points during the data denosing step.Section 2 provides a brief introduction to those morphological operators used in this paper.The new method is proposed in Section 3. The results are given in Section 4, whilst Section 5 concludes the paper.

DIFFERENTIAL PROFILES
Capital letters are used in this paper to denote binary sets and the lowercase letters denote grey-scale functions.Capital and lowercase Greek letters denote the operators.Accordingly, G = {pn} is a binary set containing N foreground points pn ∈ E, where n ∈ [1, N ], and E ⊂ Z 2 , whilst g : E → R is a regular (greyscale) grid given by a mapping function that maps a grid-space E to height-values from R. Traditional morphological operators (Shih, 2009) are defined by a structuring element w (or wS, where w is square-shaped and S is its scale), whilst attribute filters (Wilkinson, 2007) are defined by a generic attribute function Λ and an attribute threshold λ.The following notations denote the fundamental morphological operators: • Γ(G) is a binary morphological opening acting on G, • γ(g) is a grey-scale morphological opening acting on g, a binary closing that is by duality property (Ouzounis and Wilkinson, 2007) equal to a complement of the opened complement of G, and • φ(g) = −γ(−g) is a grey-scale closing, where complement is substituted by negation.
The multi-scale grid segmentation that is the focus of this paper is based on a morphological concept know as a granulometry (Maragos, 1989).A granulometry is an ordered set of morphological filters that reduces the content of a grid by filtering it on an increasing scale.When considering traditional morphological operators, a granulometry is defined as an ordered set of structuring elements w = {wS i |i ∈ [0, I]}, where S0 = 0 and Si−1 < Si.As shown by Pesaresi and Benediktsson (2001), when registering the differences between successive members of a granulometry at a point-level, a band-pass grid decomposition is achieved known as differential morphological profiles (DMPs).Denoted as δ w , DMP essentially assigns a response-vector with I elements to each point by A particular member of DMP, denoted as δ w[i] (g), is referred to as i th scale-zone of g, whilst γw S I (g) is a grid multi-scale residual.Using similar principle, Ouzounis et al. ( 2012) recently extended the concept of DMPs to differential attribute profiles (DAPs).Denoted as δ Λ λ , where λ = {λi}, λ0 = 0, and λi−1 < λi, DAP with I members is obtained by Those members of DAP, given by δ Λ λ[i] (g), are considered as attribute-zones of g and γ Λ λ I (g) is a grid's attribute residual.However, this generalisation is complicated when considering all the possible attribute functions based on which filtering can be achieved.Namely, when considering e.g.shape-attributes (for example shape-compactness), where the relationship between two connected regions C ⊆ C from a grid-space E (i.e.C , C ⊆ E) does not necessarily mean Λ(C ) < Λ(C), the subtractive filtering rule needs to be considered in order to achieve a stable decomposition (Urbach et al., 2007).For simplicity, a further condition for complementing Equation 2 is, therefore, increasing the property of Λ given by (3) A straightforward example of an increasing attribute is the area of the connected region C from a grid-space E.
Decomposition achieved by DMPs and DAPs allows for automated grid segmentation by registering two characteristic properties from each response vector: • r(g) containing the maximal response obtained at each particular point, and • q(g) containing the scale at which the maximal response has been induced.
Several approaches have been proposed based on this notion (Pesaresi and Benediktsson, 2001;Beucher, 2007;Hernández and Marcotegui, 2011;Ouzounis et al., 2012).In this paper, we define a mapping θ w (g) : g → (r(g), q(g)) at each particular point p as where is the maximum.Although the given definition is based on δ w , the same definition is used to define θ Λ λ (g) based on δ Λ λ (g).

THE EXTRACTION OF BUILDINGS
The proposed method operates in θ w (g) and θ Λ λ (g) scale-spaces in order to achieve building extraction from LiDAR data.The following four steps are used: • Initialisation arranges LiDAR points into a grid, • removal of outliers is a data denoising step, where levelling based on θ w is used to remove those features that are too small to be considered as buildings, • binarisation is a grid decomposition step, where candidates for building-regions are obtained according to their sizes and surfaces based on θ Λ λ , and • recognition of buildings, where regions' boundary-shapes are taken into account.

Initialisation
Since connectivity between points is required when using morphological operators, the input LiDAR data point-cloud is arranged into a regular grid.The grid-space E is defined by the bounding-box of the input dataset, whilst the resolution of a grid Rg is defined according to the LiDAR point-density DL as Rg = 1.0/DL, making the accuracy of the method in geometrical terms proportional to the data density.When there is more than one point within a particular grid-cell, the lowest point is selected as the representative point pn ∈ E in order to minimise the amount of vegetation and the number of other non-building points above the ground, whilst an interpolation is used to define the value of empty grid-cells, i.e. g[p * n ] = U N DEF , where p * n ∈ E and U N DEF denote the undefined values.Inverse distance weighting (IDW) was selected in our case, amongst several spatial interpolation techniques (Chaplot et al., 2006).As explained by Lloyd (2010), IDW produces a smooth oscillation-free surface without introducing any additional outliers to the data by where pn is a point from the data-sampling neighbourhood W p * n of p * n , dp n is the Euclidean distance between p * n and pn, and r is the power parameter that defines the smoothness of the interpolation.As shown by Chaplot et al. (2006), accurate results are obtained by letting r = 2, where W p * n contains no less than the three closest points.

Removal of outliers
This denoising step is primarily focused on filtering low and high outliers (Sithole and Vosselman, 2004;Mongus and Žalik, 2012).Thus, sharp peaks and steep valleys need to be considered.A levelling schema, know as MSLS-segmentation (Pesaresi and Benediktsson, 2001) is applied for their detection.Consider two DMPs obtained from g: θ w (g) containing positive response values and θ w (−g) containing negative response values (note that this is equivalent to replacing γ with φ in Equation 1 and multiplying the obtained responses by −1).At a particular point pn, the following labelling schema is used to obtain a characteristic scale q(g) of g and corresponding response r(g): Accordingly, when q(g)[pn] = q(−g)[pn], a point pn belongs to a steep valley and is a suitable candidate for a low outlier, whilst q(g)[pn] = q(g) [pn] indicates a point belonging to a sharp peak, i.e. a suitable candidate for a height outlier.When q(g)[pn] = 0, a point is not considered as a noise-point and should, therefore, not be filtered.Thus, a threshold value t[pn] at a point pn is, in our case, defined as Note that a higher threshold value is used when considering convex peaks in order to minimise the distortions of buildings' geometries caused by the filter, whilst controlling the maximal filtering scale wS I contained in w allows manipulation over the removed features.Specifically, major portions of vegetation-points can be removed by recognising them as high outliers, whilst points lying within the buildings (usually recorded when targeting glass buildings) are removed as they are recognised as low outliers.In our case, w = {w0, w1, ...wS I } is used, where SI = 3.0m was heuristically defined.A set of outliers To is then given by The outliers from To are then interpolated using eq.6 and the obtained result is shown in Fig. 1.

Binarisation
The binarisation step searches for flat regions in g with areas from a predefined range, in order to obtain a set of suitable building- candidates.A new morphological algorithm is proposed that uses differential attribute profiles δ Λ λ (g) based on attribute function Λ, defined as the sum of the second-order morphological gradients.Let functions ex and in estimate the external and internal morphological gradients, respectively (Shih, 2009).The secondorder morphological gradient that essentially captures the rate of change in the gradient of g, is then estimated as in ( ex (g)).For a particular connected region C ⊆ E, a sum of gradients Λ(C) is obtained as The points on the lower sides of the edges are emphasised, as the external gradient is computed first.Thus, flat surfaces above a neighbourhood (e.g.rooftops) have significantly lower Λ-value than those with rough surfaces (e.g.trees).Consequently, when applying attribute opening γ Λ λ based on Λ with an attribute threshold λ, regions with larger areas are removed, if their surfaces are flat, rather than when they are rough.The ratio R a λ between the area A(C) of the removed regions C and the value of λ[i] can, therefore, be used for an accurate characterisation of surfaceflatness.Based on this notion, building-extraction from LiDAR data can be achieved over the following three steps: • δ Λ λ (g) is computed first, where the following definition of λ is used λ = {0, λmin, λmin + λ∆, λmin + 2λ∆, ..., λmax}, (12) • each member δ Λ λ[i] (g) is thresholded in order to obtain a set of filtered points T i B as • finally, binary area opening ) is performed on each T i B in order to remove those regions that dissatisfy the ratio criterion R a λ .
Thus, the detection of flat-regions for building-extraction is essentially achieved based on an average second-order morphological gradient contained within a region.Although this type of average is a non-increasing attribute, the proposed framework provides an elegant solution to avoid dealing with the subtractive filtering rule that is otherwise required in these cases (Urbach et al., 2007).The results obtained in the case of the test-set from Two typical errors introduced by the proposed algorithm are highlighted in Fig. 2b.In the first (left) case, the building-region extends over the neighbouring trees, whilst a brick wall with flat top-surface was detected in the second (right) case.Both types of errors can be successfully removed during the building recognition step, applied over a set of building-candidates G, obtained by

Recognition of buildings
In the final step of the proposed method, those flat regions that do not describe buildings are removed from the binarised grid G. Firstly, any thin portions of regions that may extend from the buildings over the neighbouring objects (e.g.trees, as shown in Fig. 2b) are removed using binary morphological opening, i.e Γw 3 (G) is used in our case.Binary morphological closing Φw 3 (Γw 3 (G)) is then applied in order to remove any holes that may be present within the remaining regions.Finally, the shapes of the regions are considered and shape-compactness attribute opening is applied in order to remove long and thin regions.Let Cn ∈ G be a connected set of points from G, and n a connectedset-index.The shape-compactness Ψ(Cn) of Cn is defined as (Nixon and Aguado, 2012): where A(Cn) and P (Cn) are the area and perimetre of Cn, respectively.Shape-compactness attribute opening is then defined as Thus, the set of building regions GB is given by

RESULTS
The accuracy of the proposed method was examined on three testcases within various urban features containing buildings of different sizes and shapes.The test datasets DS1, DS2, and DS3 were located at Ljubljana -Črnuče, Medvode, and Maribor, in Slovenia, respectively (DS1 is shown in Fig. 2, whilst DS2, and DS3 can be seen in Fig. 3).The data-densities of the test-sets were as follows: 12.35 points/m 2 in the case of DS1, 9.97 points/m 2 in the case of DS2, and 6.31 points/m 2 in the case of DS3.
A different set of parameters was used for each test-case, as shown in Table 1.λmin corresponds to the area of the smallest building contained within the dataset and is slightly larger in the case of DS1 than in the cases of DS2 and DS3.The minimal compactness of the building region defined by φ does not significantly influence the results.Consequentially, φ-values are relatively close, where slight differences in the definition do not contribute more than 2% to the quality of the building-extraction, as shown in Table 2.However, the definition of λ∆ has a more significant influence on the results as it defines the area by which a building can be attached to the ground and still successfully detected.The proposed method is relatively sensitive to this particular parameter as a larger λ∆-value increases the probability of false-positives (e.g.detection of plateaus), whilst lower λ∆-values increases the probability of false-negatives (e.g.attached buildings).The quantitative evaluation of the method was achieved by comparing the extracted building-regions with the reference data obtained by the commercial software TerraScan R , where any errors were corrected manually.The following quality metrics were considered, according to the ISPRS guidelines (Rutzinger et al., 2009): where T P , F P , and F N denote true positives, false positives, and false negatives, respectively.A comparison of the results obtained by the proposed method and the initial results obtained by TerraScan R is shown in Table 2.
The evaluation of the proposed method showed an over 95% average quality of building-extraction, achieving 99.66% average completeness and 96.33% correctness.Although the average quality is comparable with TerraScan R , comparison between the completenesses and correctnesses indicates some significant differences in the performances of both methods.Namely, the proposed method achieves lower average correctness in comparison to TerraScan R , with the worst quality of building-extraction achieved in the case of DS1 (as shown in Figs. 1 and 2a).In comparison with DS2 and DS3, DS1 contained a considerable number of plateaus.Since some of them were recognised as buildings, the correctness of the building extraction decreased.However, this might be prevented by supplementing the method with an efficient ground-extraction algorithm (e.g.Mongus and Žalik (2012)).Another major contributor to errors in correctness is the fact that the proposed method does not consider a sufficient mechanism for preventing the spreading of the building-regions throughout the neighbouring objects, e.g.trees where up-to 1.5m overspreading was noticed.Note that binary opening only mitigates this problem but does not solve it.On the other hand, TerraScan R showed difficulties when detecting small buildings and often rejected the boundaries of larger ones.The average completeness of TerraScan R is, therefore, lower than the one achieved by the proposed method.However, since the proposed method is based on connected operators that are unable to breakup flat-zones, it is incapable of detecting buildings below the attaching point.Consequentially, the attached buildings remained undetected (i.e.garages placed below the surrounding terrain in the case of DS2) or only partially detected (in several cases, this effect reached up to 3.5m).This is the reason for the significant majority of errors in completeness introduced by the proposed method and can be observed in the case of DS2 (as shown in Fig. 2b).On the other hand, DS3 was included within the test dataset with the purpose of examining the behaviour of the method in regard to the shapes and sizes of the contained buildings.As shown in Fig. 3c, a great majority of errors were again related to false-positive plateaus and partially detected attachedbuildings.The results, therefore, indicate that this method is capable of detecting buildings of vastly different sizes and shapes.
As the method achieved the worst results in the case of a testset with the highest data-density (i.e.DS1), whilst the results are comparable in the cases of DS2 and DS3, there is no indication that the method is subjected to data density.

CONCLUSIONS
The paper proposes new morphological approach for the extraction of buildings from LiDAR data.This method achieves an accurate removal of high and low outliers based on MSLS-segmentation schema, whilst a new morphological algorithm is proposed for the extraction of flat regions.Differential attribute profiles are estimated based on the sum of the second-order morphological gradients, whilst area opening performed on each attribute-zone essentially leads to the thresholding flat regions according to the average contained gradient (i.e. the ratio between the sum of gradients and the area of the region).Since the errors introduced by the extractions of flat regions are characteristic, they are successfully removed according to the boundary shape during the final step of the method.Although sufficient accuracy of the proposed method was shown in comparison to TerraScan R in all the testcases, the great majority of errors could be related to particular limitations of the method.Namely, this method does show difficulties when detecting attached buildings, leading to proportional errors in the completeness of building detection, whilst errors in correctness are related to the limited detection of discontinuous terrain features (e.g.plateaus) and the overspreading of building regions throughout neighbouring vegetation.The development of sufficient solutions to these tree issues will be considered during our future work.Moreover, the mathematical stability (i.e.scale-invariance, idempotence, and increasing property) of the proposed morphological filter still needs to be proven and the extension of the method for point-clouds obtained by multi-view stereo (MVS) and structure from motion (SfM) will be considered.

Figure 1 :
Figure 1: Denoising of (a) LiDAR data generated grid, where (b) a majority of the vegetation is removed.

FigFigure 2 :
Figure 2: Detection of flat regions in (a) denoised grid g and (b) the results obtained based on δ Λ λ (g), colour-coded according to i, where λ = {0, 50, 100, ..., 25000} was used.Two typical errors are highlighted.In the left case, points from the neighbouring vegetation are accepted as building points whilst, in the right case a brick wall with a planar top surface is accepted.

Figure 3 :
Figure 3: The results obtained by the proposed method on testsets (a) DS1, (b) DS2, and (c) DS3, where (left) the heightmaps are coloured according to the ground-truth data.The yellow colour represents true-positives, the red colour false-positives, and the blue colour is used to show false-negatives.The extracted building points are (right) triangulated and visualised together with the generated DTM, based on Mongus and Žalik (2012).

Table 1 :
Parameters used for buildings-extraction from LiDAR test datasets.

Table 2 :
Quantitative evaluation of the proposed building-extraction method in comparison with TerraScan R .