A Line-based Spectral Clustering Method for Efficient Planar Structure Extraction from Lidar Data

Planar structures are essential components of the urban landscape and automated extraction planar structure from LiDAR data is a fundamental step in solving complex mapping tasks such as building recognition and urban modelling. This paper presents a new and effective method for planar structure extraction from airborne LiDAR data based on spectral clustering of straight line segments. The straight line segments are derived from LiDAR scan lines using an Iterative-End-Point-Fit simplification algorithm. Adjacency matrix is then formed based on pair-wise similarity of the extracted line segments, and a symmetric affine matrix is derived which is then decomposed into eigenspace. The planar structures are then detected by mean-shift clustering algorithm in eigenspace. The use of straight line segments facilitates the processing and significantly reduces the computational load. Spectral analysis of straight line segments in eigenspace makes the planar structures more prominent, resulting in a robust extraction of planar surfaces. Experiments are performed on the ISPRS benchmark LiDAR data over three test sites containing a variety of buildings with complex roof structures and varying sizes. The experimental results, which are quantitatively evaluated independently by the ISPRS benchmark test group, are presented. The results show that the proposed method achieves on average 80% of completeness with over 98% of correctness. Better performance is observed over larger size of buildings (>10m²) with over 92% of completeness and nearly 100% of correctness in all test areas, indicating the robustness and high reliability of the proposed algorithm.


INTRODUCTION
Light Detection and Ranging (LiDAR) has been one of the major sources for spatial information collection.The airborne LiDAR measurements are made by a scanning mechanism to generate profile slices of landscape.The recorded point clouds directly provide 3D information of the underlying terrain and objects.LiDAR has the advantages of high accuracy, rapid acquisition and high resolution which are extraordinary appropriate for large-scale elevation collection and object extraction.
Planar structures are widespread in urban landscape and are required in a variety of mapping tasks.For building extraction and urban modelling, accurate planar structures are fundamental features in data-driven reconstruction approaches (Oude Elberink and Vosselman, 2009).Automated recognition of object, such as vegetation and buildings, can be made more efficiently by employing geometric and topologic analysis of planar segments (Xu et al., 2012).Moreover, planar structures prove to be useful in data quality control (Vosselman, 2012).Multi-source registration and robot navigation also benefit from recognised planar structures.
To address the problem of converting raw LiDAR data to embedded planar surfaces in an efficient and low-cost manner, automated extraction of planar structures from airborne LiDAR data has been an active research for decades.Existing algorithms can be generally categorized as point-based and linebased.Point-based algorithms mainly utilise point coordinate attribute.Planes can be detected via region growing approaches to group locally consistent points on a large dataset (Vosselman et al., 2004;He et al., 2012) or by iterative segmentation of point clouds using RANSAC algorithms (Brenner, 2000;Schnabel et al., 2007).Point clouds can also be transformed in different feature spaces to facilitate plane detection.3D Hough transform has been proved useful by expressing points as 3D curves in Hough Space for intersection inspection (Borrmann et al., 2011).Research also shows promising results have been achieved by clustering in Gaussian Sphere feature space (Dorninger and Pfeifer, 2008;Sampath and Shan, 2010).
Unlike the point-based approaches, the other group which is also adopted in this paper employs straight line segments as the basic element for planar surface extraction.The main idea is based on the observation that, in each scan line, the points belonging to a planar surface form a straight line segment.Linebased method was introduced in Jiang and Bunke (1994) using range image, and then adopted for urban mapping (Hebel and Stilla, 2008) and robot navigation (Georgiev et al., 2011) employing LiDAR point cloud.Generally those methods follow the typical split-and-merge paradigm.In step of splitting, raw data is interpreted as straight line segments which have higher level geometric information than points.In merging step, a seed plane is derived by similarity measurement among neighbouring line segments, along with region growing.The advantages of the line-based algorithm are the high data compression rate and the higher geometric level.However, line-based approaches, to date, inherit the deficiency of greedy procedure that only one model can be detected in each iteration (Huang and Brenner, 2011).Two problems embed in this procedure.First, the remaining lines are iteratively calculated until they have been assigned and removed, which result in high computational complexity for multiple surfaces extraction.Second, line segments which belong to two adjacent planes may be too early removed by the firstly found plane.In order to tackle these issues, we introduce a novel line-based spectral clustering method for detecting and estimating a set of plane structures simultaneously in global impression.

METHODOLOGY
The proposed algorithm for planar surface extraction in airborne LiDAR data essentially includes procedures of straight line segment extraction, line-based spectral clustering and planar structure estimation which will be detailed in the following sections.

Straight line segment extraction
We adopted the Iterative-End-Point-Fit simplification scheme to extract straight line segments from each scan line, as described in Nguyen et al. (2007).Given a successively recorded points from each scan line, the algorithm starts with the single edge by simply connecting the first and last points.Then the edge is iteratively split by the furthest perpendicular distance point until all point-to-edge distances are within a user specified tolerance ( ).The output is a patch of successive connected line segments.
To avoid irrelevant line segments, three filters are defined and applied.Firstly, we eliminate the line segments with small proportion of last pulse count to point size in each line segments.Intuitively, small proportion indicates more non-last pulse points which tend to vegetation.Furthermore, we reject the line segments which contain less than four last pulses.In this way, most vegetation and small object which belong to irrelevant features can be easily filtered out.Finally, we estimate coarse ground elevation level by splitting a scan line into slices and then linking the local lowest points.Those line segments closing to ground are eliminated as ground segments.

Line-based spectral clustering
In this section, we focus on line-based clustering.To be more specific, we formulate and apply spectral clustering on straight line segments to decompose spectral space into meaningful subparts denoted as subspaces and each subspace represents a planar surface.Spectral clustering algorithm is first suggested in Donath and Hoffman (1973).It has recently been employed for data analysis such as image motion (Yan and Pollefeys, 2006) and mesh smoothness segmentation (Liu and Zhang, 2004) in computer vision.Rather than directly projecting local features in space, spectral clustering views the data clustering as a graph partitioning problem without make any assumptions on the form of the data clusters.Adjacency matrix is formed from data and then projected to the low-dimensional spectral space, in which high similarity primitives are encouraged moving toward each other while others moving increasingly apart, as Polarization Theorem stated (Brand and Huang, 2003).As a result, spectral clustering often outweighs traditional clustering algorithms, such as -means clustering (von Luxburg, 2007).The comprehensive introduction of spectral clustering can be found in (Vidal, 2011).
Generally, spectral clustering algorithms first construct affinity matrix by measuring pair-wise similarities and then use the eigenvectors of the affinity matrix to enhance clustering of data points in the eigenspace.However, line segments by themselves fail to derive co-planar similarities directly.Therefore, we design the spectral clustering flow chart as illustrated in Figure 2 and the elaborative description is given in subsections.
Figure 2: Flow chart of spectral clustering.

2.2.1
Retrieve local best-fit plane subspace.One observation is that a line segment and its nearest neighbours often belong to the same subspace.Hence, the subspace for a line segment ∈ can be obtained by the best-fit plane of line segments in its local neighbourhood.Then, if two line segments and lie in the same subspace , their locally estimated subspaces and should be similar.We adopted the method proposed by Zhang and Faugeras (1992) to define line segment neighbourhood.Particularly, line segments fall in a cylindrical space with radius r and its axis coinciding with the direction of the query line segment are defined as neighbourhoods of the segment.It means distances between the line segment and its neighbour segments should be less than r.Furthermore, we reject line segments as neighbours of the query line segment when they are in the same scan line.It is due to all line segments in the same scan line always lay in the same vertical profile.
With the neighbourhood, we then search for the locally best-fit plane for each line segment.Commencing with a line segment ∈ and its neighbourhood ⊂ , we then iterative generates candidate planes { } formed by two end-points (source and target) from and the mid-point of , ∈ .The line segment is selected as an inlier of P if offset between the mid-point of x and P is smaller than a threshold and • is closed to zero, where is the normal of and is the direction of .
with the largest number of inliers is verified as the best-fit plane of the line segment .
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey 2.2.2 Graph Laplacian construction.The similarity of each pair of subspaces is measured using angular difference and geodesic distance between the two subspaces.The angular difference is formed as Where # is the angle between the normal vectors of the two subspaces.The geodesic distance is defined as Where dist(m , S 1 ) is the Euclidean distance from mid-point m of to subspace .Then the weighted ℝ 3×3 adjacency matrix is generated by an exponential kernel as Where : and : ! are the scaling parameters of the Gaussian Kernel.Clearly, when 5 is close to 1, it means both DiffAng and DistGeo are small, then we expect and are in the same cluster.
To make the adjacency matrix unidirectional, the ℝ 3×3 affinity matrix is expressed as = (5 + 5 )/2, such that 0 < = < 1 .To normalise affinity matrix, a diagonal matrix is defined as the degree of , where is the sum of 's i-th row.Afterwards, the normalised affinity matrix denoted as symmetric graph Laplacian is expressed as Where G = /L L .Straight line segments are now represented in a ℝ 3 feature space by each row vector.

Spectral subspace clustering.
Let G and G are two points in ℝ 3 .The correlation of them matches the cosine angle between two points in feature space as With # = arccosRG S G T is the angle between two unit vectors.
In ideal case, two vectors in one cluster should have the same coefficients (M(NN = 1), such that the grouping ordered adjacency matrix appears a block-diagonal structure.Based on block-diagonal structure, its eigenvalues and eigenvectors are the union of the eigenvalues and eigenvectors of each block.However, correlations from real-world data are not exact 0 or 1, which means off-diagonal blocks are non-zero.On the other hand, the largest eigenvectors will be stable to small changes of G HIJ if their eigenvalues has small eigengap, as matrix perturbation indicated in Stewart (1999).Therefore, the union of the largest eigenvalues and eigenvectors represents the most energy-preserving of eigenvalues and eigenvectors of G HIJ (Ng et al., 2001).Furthermore, as Polarization Theorem stated in Brand and Huang (2003), the spectral decomposition to successively lower rank results in the sum of squared anglecosines ∑(cos# ) ! strictly increasing.It means the correlations will migrate away from 0 towards two poles +1 or -1, also known as truncated eigenbasis amplifying.In a low-dimensional feature space, the projected points with high correlations step together while the low correlation points move more apart.As a result, the clustering in the new spectral space representation is more likely to success.The graph matrix defined in Equation ( 4) is projected to eigenspace for clustering.The decomposition in lower dimensional is derived by choosing the largest eigenvectors from eigenvalue decomposition of L WXY .The collected eigenvectors form the new matrix Z ∈ ℝ 3× .Let [ in the -th row of Z, the point with vector [ represents the -th straight line segment in eigenspace.To project data matrix into a unit sphere, [ is normalized to be a unit vector [ \ , such that [ \ = [ /‖[ ‖.
We subsequently cluster those points in unit sphere via meanshift (Comaniciu et al., 2002) to detect planes.The selection of mean-shift clustering is that mean-shift is a non-parametric algorithm which iteratively computes the mean shift vector by translating density estimation windows until convergence.Mean-shift clustering processes without the prior knowledge of cluster quantities.Additionally, it can handle arbitrarily shaped clusters which make it more reliable for extraction applications.
We demonstrate the clustering procedure by a simple roof structure in Figure 3.The rooftop consists of 1219 LiDAR points (Figure 3(a)).The points can be represented by 82 straight line segments obtained from LiDAR scan line segmentation as shown in Figure 3(b).The line segments are further processed using the proposed algorithm.The local bestfit planes are firstly retrieved, and the weighted adjacency matrix is formed containing 82×82 similarity measurements.We only accept the weights over 0.7 and assign others as zero.
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey In Figure 3(c), the accepted coefficients are shown in grey value while the others in black (the illustration follows grouping orders).An approximately block-diagonal structure is obtained which can be easily discovered by the dominant eigenvectors.
Compared to the ideal form shown in Figure 3(d), the accuracy of derived adjacency matrix is up to 89%.Spectral clustering is then applied on the constructed graph Laplacian, resulting in four cluster groups corresponding to four planar roof structures as shown in Figure 3(e).
We use a robust estimation algorithm such as RANSAC on each segmented points to determine plane model parameters.The closed plane boundary is constructed by consequently connecting pair of end-points in each scan line.An optional polyline simplification, as described in Section 2.1, can also be used to obtain more abstract boundary in Figure 4.
Figure 5. Experiments on testing areas.The first column shows true orthoimages.The second column is the extracted plane structures on the true orthoimage.The third column illustrates the clustering result of line segments of a local area and the last column presents the extracted planar structures.

EXPERIMENTAL RESULTS
We have evaluated the proposed algorithm on International Society for Photogrammetry Remote Sensing (ISPRS) Commission WG III/4 benchmark LiDAR data in Vaihingen, Germany.On purpose of urban object classification and 3D building reconstruction, the data has 4-7 pts/m² resolution collected by Leica ALS 50 system.The entire three test sites are involved for evaluation and ground truths are detailed (LoD2) 3D models of the building roofs derived by manual stereo plotting.The first column of Figure 5 shows the true orthoimages of the three test areas respectively.Area 1 is characterized as dense development site consisting of historic buildings with complex shapes.Area 2 contains a few highrising residential buildings with mixed flat and sloped roof types.Area 3 is a purely residential area with detached buildings.For more details of site and data description, we refer to Rottensteiner et al. (2012).

Evaluation metrics
The evaluation was performed to quantify the segmentation quality of rooftop in terms of completeness and correctness (Rutzinger et al., 2009) in object level.In addition, we also ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey evaluate the geometric property of the extracted roof faces.The height difference between the extracted planar structures and the reference surfaces, denoted as ^_ ^` is computed.The other one checks the consistency of the model to the raw data, denoted as ^_ `.To be more specific, the ^_ ` measures consistency distance of every point to its extracted plane.

Results and discussion
The experimental results over three datasets are shown in the second column of Figure 5 respectively.The plane structures are represented as the 2D polygons visualising together with the orthoimage.Table 1 lists the parameter settings used in the experiment.Visual inspection illustrates that the proposed plane extraction method works well on all testing areas.The third column in Figure 5 depicts the clustering result of straight line segments from the highlighted buildings.The colour encodes the clustering status.The forth column shows the extracted planar structures from the buildings.Please note that the boundaries are derived from individual groups only without further modelling.
The extracted results were sent to the ISPRS benchmark test evaluation group, and were evaluated independently.187, 54 and 139 planes have been extracted in the three test areas, respectively.The comparison results are tabulated in Table 2.The proposed algorithm achieved over 80% completeness (Cm) for Area 1 and Area 3, while it is lower in Area 2, with a completeness of 71.0%.The performance is better over the roofs with area larger than 10m².The completeness for all areas increases up to above 90%, with over 95% in Area 1 and Area 3. The higher completeness in Area 1 and Area 3 indicates the algorithm performs better in residential areas.The evaluation shows the correctness (Cr) for all test areas is above 97%.For larger size roofs, the correctness for Area 2 and Area 3 is 100%, indicating reliable results have been achieved by the proposed algorithm.Over-segmentation (1:N) is low for all testing areas (maximum 3), while the under-segmentation (M:1) is the major problem in Area 1 and Area 3 with 40 and 44 cases respectively.The low under-segmentation in Area 2 is mainly because some roof planes with small size are not detected (low completeness).The mixed error (M:N) is also low in all testing area.The RMS errors of the height differences to the reference planes are all smaller than 0.45m.Higher completeness rate in Area 1 and Area 3 results in smaller RMS errors than in Area 2.
The RMS errors of the consistence with data are all smaller than 0.3 m.Overall, our method achieves acceptable plane extraction results.Based on noble performance on large roof planes, higher accuracy can be expected on LiDAR data with higher density.By comparing our method with state-of-art plane detection methods shown in benchmark result (Rottensteiner et al., 2012), our method achieved the top performance in terms of completeness, correctness and segmentation quality.
Some errors remain in extraction results.Clearly, some short line segments on small roof structures are eliminated in filtering step, as indicated in the circle 1 in Figure 5.A few flat and dense vegetation areas are incorrectly detected as plane.An example is given in Figure 5 (circle 2).We also notice that small planar structures with less than three line segments are hardly detected (circle 3 in Figure 5).The reason is that the plane has too few coefficients in affinity matrix to be distinguished.

CONCLUSION
This paper presents a novel algorithm for planar structure extraction from airborne LiDAR data.Rather than using a parameterized plane model for data fitting, the task of planar feature extraction is transformed into spectral subspace clustering by exploring the similarities between the straight line segments derived from the LiDAR scan lines.Therefore, line segments rather than redundant raw LiDAR points are employed as the basic elements which convey more geometric information for successive processing while significantly reducing the data volume in clustering analysis.The extracted straight line segments are projected to spectral space through graph Laplacian transformation and spectral decomposition.This enhances the correlation of high-similarity points while penalizes dissimilar points allowing more efficient plane detection via mean-shift clustering algorithm in feature space.
The proposed method has been experimentally evaluated on the ISPRS benchmark LiDAR data over three test sites.The independent evaluation by the ISPRS benchmark test group showing on average 80% of completeness with over 98% of correctness has been achieved.The errors are observed on dense vegetation with flat top which are incorrectly detected as roofs.
In addition, buildings with small size are not detected due to insufficient straight line segments on roofs.This is indicated in the evaluation where the completeness increases up to 92% over large size of buildings.It is also noted that the correctness is very high over all test areas.The close to 100% correctness on large roofs demonstrates the robustness and high reliability of the proposed method.
Current research centres on the refinement the developed algorithm and modelling of the extracted planes which will be further processed for automated 3D building reconstruction in our project.
Figure 1 (a) shows a portion of LiDAR raw data in one scan line.Figure 1 (b) presents the line simplification result and coarse ground level shown in white dashed lines.Figure 1(c) demonstrates the filtering result without irrelevant line segments.The line colours encode number tags.

Figure 1 :
Figure 1: Straight line segment extraction.(a) Point cloud from a portion of a scan line; (b) initial straight line segments and coarse ground level; (c) Straight line segments on roofs.

Figure 4 :
Figure 4: Boundary dection results.(a) and (b) Detected boundaries of extracted planes and final structure representation; (c) and (d) simplified boundaries and planar structure.

Table 1
Parameters setting for experiment

Table 2
Statistical evaluation of the plane extraction results