Classiﬁcation of ALS Point Cloud with Improved Point Cloud Segmentation and Random Forests

: This paper presents an automated and effective framework for classifying airborne laser scanning (ALS) point clouds. The framework is composed of four stages: (i) step-wise point cloud segmentation, (ii) feature extraction, (iii) Random Forests (RF) based feature selection and classiﬁcation, and (iv) post-processing. First, a step-wise point cloud segmentation method is proposed to extract three kinds of segments, including planar, smooth and rough surfaces. Second, a segment, rather than an individual point, is taken as the basic processing unit to extract features. Third, RF is employed to select features and classify these segments. Finally, semantic rules are employed to optimize the classiﬁcation result. Three datasets provided by Open Topography are utilized to test the proposed method. Experiments show that our method achieves a superior classiﬁcation result with an overall classiﬁcation accuracy larger than 91.17%, and kappa coefﬁcient larger than 83.79%.


Introduction
Commercial Airborne Laser Scanning (ALS) systems emerged in the mid-1990s for bathymetric and topographic applications.With the aid of direct geo-referencing technique, laser scanning equipment installed in the aircraft collect a cloud of laser range measurements for calculating the 3D coordinates (xyz) of the survey area [1].In contrast to the 2D remote sensing imagery, an ALS point cloud is a swarm of points with XYZ coordinates [2], and thus describes the 3D topographic profiles of natural surfaces.Moreover, ALS point clouds have other benefits such as no effects of relief displacement, penetration of vegetation, and insensitivity to lighting conditions [1].Therefore, ALS technique has been effectively used for ground point detection [3][4][5][6][7], topographic mapping [8], 3D city modelling [9][10][11][12][13], object recognition [14][15][16], solar energy estimation [17], etc.
Over the last two decades, significant contributions to the consolidation and extension of ALS data processing methods have been witnessed [1].Among these processing methods, classifying the ALS data into categorical object instances is the first and most critical step for further data processing and model reconstruction [18].Based on the granularity of basic processing units, these existing classification strategies can be categorized into three groups, i.e., point-based classification [18][19][20], segment-based classification [21][22][23][24], and multiple-entity-based classification [25].A brief description of these existing methods is provided as follows.

Point-Based Classification
This kind of classification has attracted the most majority of research works in contrast to the other two kinds of classification strategies.In the process of point-based classification, ALS features of individual points [19] are firstly extracted.Then a classifier such as JointBoost [18] is trained using a number of selected training samples.Finally, the input ALS point cloud is classified via the trained classifier and the extracted features.
Additionally, to compute the features of individual points, a respective neighborhood definition is required to describe the local 3D structure around each individual point.Generally, there are three kinds of neighborhoods, i.e., spherical neighborhood [26], cylindrical neighborhood [27], and k-closest neighborhood [28].Among the three neighborhoods, the scale parameter, either a fixed radius or a constant value k, is required.Due to the variation of local 3D structures and point densities, the constant scale parameter often fails to describe the local structural configurations.Thus, more and more studies such as [18,[29][30][31][32][33][34] focus on seeking an optimal neighborhood size for each individual point.Unfortunately, these neighborhood optimization methods require repetitive calculations of eigenvectors and eigenvalues for each point, therefore they are rather time-consuming [35], which is the main disadvantage of this kind of classification.

Segment-Based Classification
Point cloud segmentation has been involved in ALS point cloud classification since its emergence.Generally, segment-based classification methods first perform segmentation on the point cloud after removing the ground points [21].Then, the non-ground points are segmented into a number of segments, and features are extracted for each segment.Finally, a fuzzy model classifier [21,36] or several classification rules [22,24] are utilized to classify the segments.However, most of these studies are for non-ground points, and none of them uses Random Forests (RF) for feature selection and classification.
In addition, segment-based classification relies heavily on its employed segmentation method.A variety of point cloud segmentation methods have been proposed, which can be roughly classified as model-fitting-based methods, region-growing-based (RG-based) methods, and clustering-featurebased methods [37].However, these existing methods segment input 3D point clouds into only one type of geometric structure.Actually, point clouds consist of a variety of geometric structures, such as planes, smooth surfaces and rough surfaces.In a complex 3D scene, there may exist regular and irregular man-made objects, and natural objects.Regular man-made objects such as buildings are composed of planar surfaces and smooth surfaces, while irregular man-made objects such as cars and natural objects like trees are composed of rough surfaces.
Therefore, segmenting point clouds into only one type of geometric structure is unreasonable.For example, existing planar segmentation methods segment all the points in an input point cloud into planes.If points are on building roofs, these methods are logical and perform well, however, if the points are on trees or cars, these methods which roughly segment these points into false planes are illogical.To obtain a superior classification result, we should consider a query point's geometric structure, and then segment it into a planar surface, smooth surface, or rough surface.
Although the aforementioned limitations exist, segment-based methods still have two main benefits in contrast to point-based classification methods, i.e., (i) segments are helpful to compute geometric features which relieve the dependence on neighborhood optimization [18,34] methods, and (ii) segments give several new attributes which are helpful to employ semantic rules.

Multiple-Entity-Based Classification
Multiple-entity-based classification [25] is considered as a combination of the segment-based and point-based classification.To solve the problem that a complex 3D scene is difficult to be characterized by only individual points or one kind of segments, this method utilizes three kinds of entities, i.e., points, planar segments, mean shift segments.In the process of classification, the input ALS point cloud is first divided into ground points and non-ground points.Next, planar segments are extracted from the non-ground points, and the scattered points are remained.Then, the planar segments are classified into several classes.The remained points are point-wise classified based on the contextual information offered by the classified planar segments.Finally, in complex areas where vegetation covers building roofs, mean shift segments are extracted to classify these areas.
However, the process of this method is a hierarchical classification procedure, which involves many steps.Besides, the mean shift segments and planar segments are derived from different segmentation methods, which adds additional classification steps.To simplify the classification process, a point cloud segmentation method that is able to extract more than one kind of segments is required.
These above three strategies have two common elements, including feature extraction and classifiers.Therefore, we present a brief description of both them as follows.

Feature Extraction
There are three main groups of features for ALS point cloud classification, i.e., reflectance-based features, descriptor-based features, and geometric features.

•
The reflectance-based features are often related to the intensity [38] and echo [18] recorded by scanner systems.Therefore, the distinctiveness of this kind of features relies heavily on the quality of the scanner's signal.

•
Common geometric features are height-based features [19], eigenvalue-based features [19,36], projection-area-based features [18,39], surface-based features [18], etc.Specifically, the eigenvalue-based features derived from the 3D structure tensor which is represented by the 3D covariance matrix derived from the 3D coordinates of all points within a local neighborhood, are discriminative in a variety of classification approaches.In contrast, the geometric features are deeply studied and widely used by state-of-the-art methods.
Most existing studies often compute as many features as possible to obtain a superior classification result.When a large number of features are extracted, some of them may be redundant.These redundant features not only increase the computational burden, but also waste the memory space.Therefore, recent studies introduce a feature selection procedure [19,20,38,44] as an additional step between feature extraction and classification steps.

Classifiers for ALS Data Classification
In the classification stage, many studies have tried locally independent classifiers, such as Support Vector Machine (SVM) [45], Adaptive Boosting (AdaBoost) [46], Expectation Maximum (EM) [47], RF [48,49], JointBoost [18], etc.The fundamental idea is to train a classifier by using given training samples which is used for prediction later [20].Specifically, due to the excellent performance, the RF classifier [50] has received increasing attention [51].Some studies [19,20,34,38] have looked into the potential of the RF classifier to improve urban objects classification and select uncorrelated features for ALS point clouds.
However, the integration of RF and the segment-based classification is rarely studied, as well as the importance analysis of segment features.In addition, the robustness of the classification methods is rarely analyzed when noises exist in the extracted features.
In this paper, we focus on the segment-based classification due to its advantages over the point-based classification.To address the aforementioned problems, we design a segment-based classification framework.This framework has three improvements compared to the existing methods: (1) A novel point cloud segmentation method is proposed.This method clusters the points with regular neighborhoods into planar and smooth surfaces, and the points with scattered neighborhoods into rough surfaces.(2) RF is integrated with the segment-based classification to select features and perform classification.
The integration of RF and the segment-based classification improves the robustness of ALS point cloud classification.(3) Semantic rules are employed to optimize the classification result.The semantic rules are more convenient to be detected when we process segments.
The outline of this paper is shown as follows.Section 2 presents the methodology of our proposed classification framework which contains a novel point cloud segmentation method, feature extraction based on segments, the integration of RF and segment-based classification, and post-processing based on semantic rules.The experiments and discussions are presented in Section 3, followed by Section 4 which summarizes the uncertainties, errors and accuracies of the proposed classification framework.The research conclusions are presented in Section 5.

Methodology
The proposed classification framework is composed of four stages as shown in Figure 1.First of all, a step-wise point cloud segmentation method which is able to cluster points with different neighborhoods into different geometric structures is proposed (see Section 2.1).Next, a segment rather than an individual point is considered as the basic processing unit for feature extraction (see Section 2.2).Then, we employ RF to select uncorrelated features based on a backward elimination method [52], and improve the robustness of ALS point clouds classification (see Section 2.3).Finally, we utilize semantic rules to optimize the classification result in the post-processing stage (see Section 2.4).
Remote Sens. 2017, 9, 288 4 of 34 (1) A novel point cloud segmentation method is proposed.This method clusters the points with regular neighborhoods into planar and smooth surfaces, and the points with scattered neighborhoods into rough surfaces.(2) RF is integrated with the segment-based classification to select features and perform classification.The integration of RF and the segment-based classification improves the robustness of ALS point cloud classification.(3) Semantic rules are employed to optimize the classification result.The semantic rules are more convenient to be detected when we process segments.
The outline of this paper is shown as follows.Section 2 presents the methodology of our proposed classification framework which contains a novel point cloud segmentation method, feature extraction based on segments, the integration of RF and segment-based classification, and postprocessing based on semantic rules.The experiments and discussions are presented in Section 3, followed by Section 4 which summarizes the uncertainties, errors and accuracies of the proposed classification framework.The research conclusions are presented in Section 5.

Methodology
The proposed classification framework is composed of four stages as shown in Figure 1.First of all, a step-wise point cloud segmentation method which is able to cluster points with different neighborhoods into different geometric structures is proposed (see Section 2.1).Next, a segment rather than an individual point is considered as the basic processing unit for feature extraction (see Section 2.2).Then, we employ RF to select uncorrelated features based on a backward elimination method [52], and improve the robustness of ALS point clouds classification (see Section 2.3).Finally, we utilize semantic rules to optimize the classification result in the post-processing stage (see Section 2.4).

Step-Wise Point Cloud Segmentation
The proposed point cloud segmentation method is a RG-based one, and it clusters the points into planar surfaces, smooth surfaces, or rough surfaces.
Our segmentation procedure consists of three steps: region growing (RG) with RANdom SAmple Consensus (RANSAC), scattered points clustering, and small segments merging (see Figure 2).The first step extracts planar and smooth surfaces, and recognizes scattered points from the input point cloud.Then, the scattered points clustering step extracts rough surfaces from the scattered points.At last, an optional step is performed to merge small segments.

Step-Wise Point Cloud Segmentation
The proposed point cloud segmentation method is a RG-based one, and it clusters the points into planar surfaces, smooth surfaces, or rough surfaces.
Our segmentation procedure consists of three steps: region growing (RG) with RANdom SAmple Consensus (RANSAC), scattered points clustering, and small segments merging (see Figure 2).The first step extracts planar and smooth surfaces, and recognizes scattered points from the input point cloud.Then, the scattered points clustering step extracts rough surfaces from the scattered points.At last, an optional step is performed to merge small segments.
Remote Sens. 2017, 9, 288 5 of 34 It is notable that rough surfaces are separated from the other types of segments by the aforementioned steps.However, planar surfaces and smooth surfaces are not distinguished when we extract them in RG with RANSAC step.For the application of semantic recognition, planar surfaces and smooth surfaces are easily to be distinguished by their curvatures.In this study, the feature extraction, and the RF-based feature selection and classification stages perform without distinguishing these three types of segments.

RG with RANSAC
The plane-based RANSAC algorithm is widely used in point cloud segmentation tasks [53].However, there is rare study employing RANSAC to improve RG-based segmentation methods.
To extract planar surfaces, two difficulties should be overcome for a RG-based method, i.e., nonoptimal segmentations around edges where two surfaces meet [54], and the detection of small or narrow planes [55].The integration of RANSAC and RG is able to solve both problems.In addition, we utilize a smooth RG procedure, which is able to extract smooth surfaces simultaneously.
There are three substeps in the RG with RANSAC step, i.e., normal estimation, RG with redefined constraints, and small segment elimination.
(1) Normal estimation We employ RANSAC-Normal [56] to the RG-based method.Our previous approach [56] has validated that the RANSAC-Normal is efficient to extract a suitable plane from a complex neighborhood with intersecting surfaces.This procedure first determines the neighbors of the -th query point , then fits a local plane based on the RANSAC algorithm, and finally, defines the normal of the fitted plane as the RANSAC-Normal.
In addition, during the normal estimation procedure, a number of scattered points are detected, and they are stored in a scattered point set .The pseudo code which shows details of this procedure, is presented in Algorithm 1.The parameter is utilized to determine how many neighbors of a query point will be detected in Row 4. The parameter is the threshold of the plane-based RANSAC algorithm, which is utilized in Row 5.It is notable that rough surfaces are separated from the other types of segments by the aforementioned steps.However, planar surfaces and smooth surfaces are not distinguished when we extract them in RG with RANSAC step.For the application of semantic recognition, planar surfaces and smooth surfaces are easily to be distinguished by their curvatures.In this study, the feature extraction, and the RF-based feature selection and classification stages perform without distinguishing these three types of segments.

RG with RANSAC
The plane-based RANSAC algorithm is widely used in point cloud segmentation tasks [53].However, there is rare study employing RANSAC to improve RG-based segmentation methods.
To extract planar surfaces, two difficulties should be overcome for a RG-based method, i.e., non-optimal segmentations around edges where two surfaces meet [54], and the detection of small or narrow planes [55].The integration of RANSAC and RG is able to solve both problems.In addition, we utilize a smooth RG procedure, which is able to extract smooth surfaces simultaneously.
There are three substeps in the RG with RANSAC step, i.e., normal estimation, RG with redefined constraints, and small segment elimination.
(1) Normal estimation We employ RANSAC-Normal [56] to the RG-based method.Our previous approach [56] has validated that the RANSAC-Normal is efficient to extract a suitable plane from a complex neighborhood with intersecting surfaces.This procedure first determines the k n neighbors of the i-th query point p i , then fits a local plane based on the RANSAC algorithm, and finally, defines the normal of the fitted plane as the RANSAC-Normal.
In addition, during the normal estimation procedure, a number of scattered points are detected, and they are stored in a scattered point set SP.The pseudo code which shows details of this procedure, is presented in Algorithm 1.The parameter k n is utilized to determine how many neighbors of a query point will be detected in Row 4. The parameter d r is the threshold of the plane-based RANSAC algorithm, which is utilized in Row 5. Five sets (NOR, I N, PR, SP, RP) are generated in Algorithm 1, and they are useful for the RG with redefined constraints.
(2) RG with redefined constraints This procedure is similar to the RG method presented in [57].However, two constraints (local connectivity and surface smoothness [57]) are redefined based on the plane-based RANSAC algorithm.

• Constraint 1: local connectivity
The points in a segment should be locally connected.In contrast to literature [57], we utilize the inliers I N i to optimize this constraint.If p i is in I N i , the points in I N i are local connective to it.Otherwise, there is no point local connective to p i .

• Constraint 2: surface smoothness
The points in a segment should locally make a smooth surface, whose RANSAC-Normals do not vary "too much" from each other.This constraint is expressed through dot product between normals: where th θ is the threshold of the constraint, nor s is the normal of current seed, and nor n is the normal of a point which is local connective to the current seed.
The pseudo code which shows details of the RG with redefined constraints procedure, is presented in Algorithm 2. The parameter th θ is utilized to restrict the dot product between normals in the Row 12 and its usage is presented in Formula (1).

8:
Find Current inliers of SC j I N j ← {IN} (constraint 1)

9:
For k=0 to size(I N j ) do 10: Neighbor point index pn k ← I N j (k) (3) Small segment elimination If the point density of trees is dense, there may be some small segments (planar and smooth surfaces) in the tree areas.Therefore, we should remove these small segments from FP, and add the points in these small segments to SP.There is a parameter in this procedure, i.e., the minimum size threshold s min .s min is expressed via the number of points in a segment.Specifically, The RG with RANSAC step divides the input 3D point cloud into two point sets.The first set is the regular point set RP, and the second set is the scattered point set SP. Points in RP are clustered into planar and smooth surfaces, and points in SP are clustered into rough surfaces in the subsequent step.In the RG with RANSAC step, there are two procedures for dividing regular points and scattered points, which are shown as follows: • For a query point p i in Algorithm 1, we first find k n neighbors of it, and then plane-based RANSAC is performed to determine inliers PI i which is on the local fitted plane.If p i is not in PI i , it will be recognized as a scattered point.

•
When the Algorithm 2 is performed, regular points are recognized and clustered into a number of segments containing planar and smooth surfaces.There may be misjudgment if a planar or smooth surface is small enough.Therefore, the small segment should be removed and points in it will be recognized as scattered points.
The result of the RG with RANSAC step is presented in Figure 3a.Points rendered in black are scattered points and other points are regular points.

Scattered Points Clustering
In this step, the scattered points in are first clustered into a number of initial patches, and then rough surfaces are extracted by growing these initial patches.The two substeps are detailed as follows: (1) Initial patch construction This substep is an iterative procedure.An initial patch is iteratively extracted from the scattered point set , until all the points in have been traversed.There are two parameters in this procedure, i.e., the maximum size threshold and the maximum distance threshold .The number of points in an initial patch has to be smaller than .Moreover, an initial patch has to satisfy the follow constraint:

Scattered Points Clustering
In this step, the scattered points in SP are first clustered into a number of initial patches, and then rough surfaces are extracted by growing these initial patches.The two substeps are detailed as follows: (1) Initial patch construction This substep is an iterative procedure.An initial patch is iteratively extracted from the scattered point set SP, until all the points in SP have been traversed.There are two parameters in this procedure, i.e., the maximum size threshold k p and the maximum distance threshold d max .
The number of points in an initial patch has to be smaller than k p .Moreover, an initial patch has to satisfy the follow constraint: (2) where (x i , y i , z i ) is a point in an initial patch, and (x c , y c , z c ) is the centroid point of the initial patch.The details of this procedure is shown in Algorithm 3. To stabilize this procedure, the curvature of each point in {SP} is estimated.These points in {SP} are sorted according to their curvatures from minimum to maximum, and then are traversed in order.As shown in Algorithm 3, the parameter k p is utilized in Row 2 and 6 to determine how many neighbors of a query point will be detected for kd_tree search.The parameter d max is utilized in Row 7 to restrict the distance between a neighbor and the query point.Only neighbors with the distance smaller than d max are extracted to construct an initial patch.After all the initial patches are extracted, their covariance matrices are computed (Row 17).Let p i = (x i , y i , z i ) for i = 1, 2, . . ., N, be the points in an initial patch, the covariance matrix M is defined as: where p u is the mean vector of all the points in the patch.
After all the covariance matrices are determined, these initial patches in PT are sorted by the determinants of their covariance matrices.Finally, each covariance matrix is normalized by its determinant.The result of the initial patch construction is shown in Figure 3b.
(2) Patch growing The patch growing is similar to the RG with redefined constraints procedure.Therefore, we do not present the pseudo code in the following.In this procedure, the growing unit is an initial patch rather than an individual point.Each initial patch is considered as a seed and grew in the order obtained in the initial patch construction procedure (Row 18 of Algorithm 3).Two constraints for patch growing are defined as follows: Only adjacent patches of the seed patch can be added into the current segment.
• Constraint 2: geometrical similarity The geometrical difference D pa of two patches in a segment has to be smaller than a threshold th sp .D pa is defined by using log Euclidean Riemannian metric [58] which measures how close two covariance matrices are.Given two covariance matrices M 1 and M 2 , D pa is computed as follows: where log(•) is the matrix logarithm operator and ||•|| F is the Frobenius norm.
After these two steps have been performed, points are clustered into three kinds of segments, i.e., planar, smooth, and rough surfaces.If small segments with a size smaller than the minimum size threshold s min still exist, we merge them with their nearest neighbor segments.This step, which is an optimizing procedure, will iterate until all the segments are traversed.The result of the step-wise point cloud segmentation method is shown in Figure 3c.

Employed Features and Their Calculation
The difference of our method from the others is that we extract features of segments rather than individual points.Herein, we focus on four groups of geometric features, namely projection-area-based ones [18], eigenvalue-based ones [19], elevation-based ones [18], and other ones.

Projection-Area-Based Features
Projection-area-based features are first proposed in [39], and then are applied to 3D point cloud classification [18].In this paper, we borrow the idea from literature [18] which is shown in Figure 4.However, the difference of our method is that the basic processing unit is a segment rather than an individual point.A segment has no fixed size compared to the neighborhood of an individual point.Therefore, a larger segment has larger projection area than a smaller segment.This problem affects the distinctiveness of this kind of features.We define two ratios to overcome this problem.The first one is the tangent plane projection ratio PR t , and the second one is the horizontal projection ratio PR h .

•
Tangent plane projection ratio PR t First of all, a covariance matrix is computed from all points in a segment.The normal vector is determined via the eigenvector corresponding to the lowest eigenvalue and a new 2D coordinate system with two coordinate axis corresponding to the largest and middle eigenvectors is constructed.Next, we project all the points in the segment along the normal direction into the coordinate system.The maximum and minimum coordinates in the coordinate system are determined, and a 2D grid is constructed by a given gird bin size.Then, we determine the number of non-bare bins that have projection points as the tangent plane projection area PA t .Finally, the ratio PR t is defined as: where PAA t is the number of all bins in the 2D grid.

• Horizontal projection ratio PR h
This feature extraction procedure is similar to that of the tangent plane projection ratio, but we select the direction parallel to the Z-axis as the normal direction.Next, we construct 2D grid in XY-plane, and then accumulate the number of non-bare bins as the horizontal projection area PA h .Finally, the ratio PR h is defined as PR h = PA h /PAA h , where PAA h is the number of all bins.Figure 5a depicts this feature.In this paper, the bin size is set to 0.2 m for both the two projection-area-based features.
The difference of our method from the others is that we extract features of segments rather than individual points.Herein, we focus on four groups of geometric features, namely projection-areabased ones [18], eigenvalue-based ones [19], elevation-based ones [18], and other ones.

Projection-Area-Based Features
Projection-area-based features are first proposed in [39], and then are applied to 3D point cloud classification [18].In this paper, we borrow the idea from literature [18] which is shown in Figure 4.However, the difference of our method is that the basic processing unit is a segment rather than an individual point.A segment has no fixed size compared to the neighborhood of an individual point.Therefore, a larger segment has larger projection area than a smaller segment.This problem affects the distinctiveness of this kind of features.We define two ratios to overcome this problem.The first one is the tangent plane projection ratio , and the second one is the horizontal projection ratio .

 Tangent plane projection ratio
First of all, a covariance matrix is computed from all points in a segment.The normal vector is determined via the eigenvector corresponding to the lowest eigenvalue and a new 2D coordinate system with two coordinate axis corresponding to the largest and middle eigenvectors is constructed.Next, we project all the points in the segment along the normal direction into the coordinate system.The maximum and minimum coordinates in the coordinate system are determined, and a 2D grid is constructed by a given gird bin size.Then, we determine the number of non-bare bins that have projection points as the tangent plane projection area .Finally, the ratio is defined as: where is the number of all bins in the 2D grid.

 Horizontal projection ratio
This feature extraction procedure is similar to that of the tangent plane projection ratio, but we select the direction parallel to the Z-axis as the normal direction.Next, we construct 2D grid in XY-plane, and then accumulate the number of non-bare bins as the horizontal projection area .Finally, the ratio is defined as , where is the number of all bins.Figure 5a depicts this feature.In this paper, the bin size is set to 0.2 m for both the two projection-area-based features.

Elevation-Based Features
In this part, we define four elevation-based features, the latter three borrow the idea from the height-based features for individual points presented in [18].

•
Relative elevation H r We first determine the adjacent segments for a query segment S q in the XY-plane.Neighbors of all the points in S q are searched in XY-plane by a kd_tree with a given distance threshold.
The segments S i a for i = 1, . . ., N which contain the searched neighbors are the adjacent segments of S q .Next, a point p i a which is closest to S q is determined in S i a .Then, a point p i q which is closest to S i a is determined in S q , and p i a and p i q construct a neighboring point pair which represents the relationship between S i a and S q .Finally, the relative elevation H r is defined as: where Z p i q is the z-coordinate of p i q , and Z p i a is the z-coordinate of p i a .H r is shown in Figure 5c.

• Elevation variance H v
This feature H v is the variance of elevation values of all points in each segment.H v is computed as: where N is the total number of points in the segment, H ave is the average elevation of all points in the segment, H i is the elevation of the i-th point in the segment.

• Elevation difference H di f f
This feature is the difference between the highest elevation H highest and the lowest elevation H lowest of a segment.In other words, the elevation difference H di f f is computed as This feature is the elevation difference between the centroid point and the lowest point of a segment.The normalized elevation H n is defined as: where H ave is the elevation of the centroid point, and the H lowest is elevation of the lowest point.

• Area O ar
The number of points in a segment is defined as area feature O ar , which reflects the area of the segment.

• Slope O sl
This feature O sl is computed as the included angle between the normal of a segment and a vertical vector.O sl is shown in Figure 5d.
Overall, the extracted four groups of features are list in Table 1.The RF classifier [50] is an ensemble of a set of decision trees.These trees in RF are created by drawing a subset of training data through a bagging approach.The bagging randomly selects about two thirds of the samples from a training data to train these trees.This means that the same sample can be selected several times, while others may not be selected at all [51].Then, the remaining samples are used in an internal cross-validation technique for estimating performs of RF.In addition, the Weighted Random Forest [59] method is utilized for solving the imbalanced sample problem in RF.
Two parameters, i.e., the number N tree of trees and the number N f eas of features, are required for using a RF classifier.Then, each tree in RF is independently produced without any pruning.The number N f eas of features is used for training each tree.Each node in a tree is split by selecting N f eas features from the d-dimensional input feature space at random.The splitting function usually uses Shannon entropy or Gini index as a measure of impurity.In prediction, each tree votes for a class membership for each test sample, the class with maximum votes will be considered as the final class.

Feature Selection
The objective of feature selection is to identify a small set of discriminative features that can still achieve a good predictive performance [19].The RF provides a measure V I i of variable importance based on averaging the permutation importance measure of all the trees which is shown to be a reliable indicator [60].The permutation importance measure is based on Out-Of-Bag (OOB) errors, and is utilized to select features.
Herein, we use the variable importance measure in RF and the backward elimination method [52] to select features.The backward elimination method removes the most relevant features by iteratively fitting RF.In our approach, only one feature with the lowest importance value is eliminated at each iteration, and then a new forest is built by the remained N rm f features.At the end of each iteration, we compute the mean decrease permutation accuracies [19] V I i for i = 1 . . .N rm f and rank the remained features by them.To measure the importance of the remained features at each iteration, the overall mean decrease permutation accuracy OD k is computed by averaging all the remained features' importance values.OD k is computed as follow: where k corresponds to the iteration times, and N bi is the total number of iterations.The iterative procedure stops and all the RFs are fitted when N rm f is equal to N f eas .After all the RFs are fitted, we computed the range I r between the maximum and the minimum overall mean decrease permutation accuracies (OD max = max(OD k ) f or k = 1 . . .N bi , and OD min = min(OD k ) f or k = 1 . . .N bi ).Then, we select a critical point according to the variation tendency of OD k .In this paper, the principle of the selection is that the variation of OD k caused by the eliminated features should be lower than 20%, i.e., the critical point which divides the range I r by a ratio of at least 8:2 is selected.At last, we select the most important features according to the critical point and the backward elimination results.

Supervised Classification
After the feature selection, the training samples with selected feature set are utilized to train a RF classifier, and the classifier is used to predict the labels of unlabeled segments.

Post-Processing
There may be some misjudgments in the above initial classification results.Actually, both natural and artificial objects generally have been associated with specific semantic information.For example, cars occupy small areas and are parked on ground, most of building roofs are composed of planar faces and have relatively large area, and wires are elevated over ground.Moreover, semantic information is more convenient to be detected via segments than individual points.Therefore, semantic information is utilized to define several rules for optimizing the initial results.Note that the post-processing stage cannot correct all the misjudgments.
To define the semantic rules, two types of neighborhoods are determined.The first one Ngbr 1 is determined in 3D space, and the second one Ngbr 2 is determined in XY-plane.In the post-processing stage, if we find a misjudged segment based on these rules, we first detect labels of its adjacent segments based on Ngbr 1 , then we relabel it as the class which arises the most times.Herein, we only list the useful rules for our approach.The rules for each class are shown as follows.

•
In the neighborhood Ngbr 1 of a query ground segment S gs , S i ga1 for i = 1, 2, . . ., N s is the i-th adjacent segment.The elevation difference between each segment pair containing S gs and S i ga1 is in a small interval Int gs .

•
In the neighborhood Ngbr 2 of a query ground segment S gs , there is no adjacent segment S i ga2 , i = 1, . . ., N s whose maximum elevation much lower than the minimum elevation of S gs .

•
Considering an extreme case, if we want to obtain a high precision DTM from a complex mountainous region, we can combine the initial ground segments with the segment-based PTD (progressive TIN densification) filtering method [61] or a progressive graph cut method [62].First, the ground segments in the initial classification result are considered as latent ground segments.Next, a 2D grid is constructed with a bin size equal to the maximum building size.Then, a latent ground segment which contains the lowest point in each bin is recognized as a ground seed segment.Finally, points in the ground seed segments are utilized to construct a TIN as the initial ground surface, and ground points will be extracted by a PTD or a progressive graph cut method.
(2) Rule for building roofs In the neighborhood Ngbr 1 of a query building roof segment S bs , it has the same characteristics as ground.However, in the neighborhood Ngbr 2 , we can find an adjacent segment S i ba2 , i = 1, . . ., N s whose maximum elevation is much lower than the minimum elevation of S bs , the elevation difference threshold of this elevation difference is denoted as D mm .

•
A car segment should have a certain range of area: [Ar vmin , Ar vmax ].

•
In the neighborhood Ngbr 1 of a query car segment S vs , there is at least one adjacent segment S i va1 , i = 1, . . ., N s labeled as ground.

•
In the neighborhood Ngbr 1 of a query wire segment S ps , there is no adjacent segment S i pa1 , i = 1, . . ., N s with a label of ground or car.

•
The 3D shape of a wire segment can be linear or planar, however, cannot be volumetric.

•
Vegetation segments are often confused with vehicle segments.Besides, the misjudged segments in the initial classification result tend to arise in high vegetation rather than low vegetation.Therefore, the misjudged segments in high vegetation are able to be corrected by the rules for vehicles simultaneously., the segment will be relabeled as vehicle.The thresholds of these ranges are cited from the approach in literature [24].Besides, the height, width and length of a segment are computed in a local coordinate framework composed of v 1 , v 2 and v 3 .The v 1 , v 2 and v 3 are eigenvectors of the covariance matrix M constructed by points in the segment as formula (3).

Experiments and Discussions
We developed a protype framework for the proposed segment-based point cloud classification method using C++ language and Point Cloud Library (PCL) [63].We also implemented segment-based point cloud classification using SVM and the RG segmentation [57], in which the open source libSVM [64] is used for the implementation of SVM and PCL is used for the implementation of RG segmentation.
The experiments are conducted on a workstation running Microsoft Windows 7 (×64) with two 16-Core Intel Xeon E5-2650, 64GB Random Access Memory (RAM) and 3TB hard disk.
There are two parts in the experiments and discussions.The first part is experimental setting which includes study areas and evaluation metrics.The second part is results and discussions which are presented in Sections 3.2-3.6.
In the second part, there are three improvements compared with existing classification methods, i.e., the step-wise point cloud segmentation, the integration of RF and the segment-based classification, and the employment of semantic rules in the post-processing stage.Therefore, we first discuss the impact of the three improvements and analyze the advantages of them in Sections 3.2-3.4.Then, the classification results and accuracies of the protype framework are shown in Sections 3.5 and 3.6.Note that the classification results should be presented after the discussion of the three improvements.Because the classification results are obtained by the parameters which are determined by the discussion in Sections 3.2-3.4.

Study Areas
Three study areas are involved in our experiments.The first one is selected from a publicly available ALS dataset which is obtained by the University of Iowa in 2008 [65].The data are collected to survey the Iowa River Flood along the Iowa River and Clear Creek Watershed.The data collection is funded by NSF Small Grant for Exploratory Research (SGER) program.Area 1 is shown in Figure 6a and it contains 1,512,092 points with an average point spacing of 0.6 m.In Figure 6a, the point cloud of Area 1 is colored by elevation.In this area, the ground is flat and smooth.On the ground, there are some parking lots where many vehicles are parked.In addition, some vehicles are parked under The second one is selected from a publicly available ALS dataset which is collected in Sonoma County [66] between September 28 and November 23, 2013 by the Watershed Sciences, Inc. (WSI).The dataset is provided by the University of Maryland and the Sonoma Country Vegetation Mapping and Lidar Program under grant NNX13AP69G from NASA's Carbon Monitoring System.Area 2 is shown in Figure 6b and it contains 1509228 points with an average point spacing of 1.0 m.In figure 6b, the point cloud of Area 2 is colored by elevation.In this area, there is a mountain which is full of trees and ornamented by several houses.Buildings are overlapped with tree canopies significantly.Moreover, building elements such as skylights have complex 3D structure.Wires go across the trees The second one is selected from a publicly available ALS dataset which is collected in Sonoma County [66] 6b and it contains 1509228 points with an average point spacing of 1.0 m.In Figure 6b, the point cloud of Area 2 is colored by elevation.In this area, there is a mountain which is full of trees and ornamented by several houses.Buildings are overlapped with tree canopies significantly.Moreover, building elements such as skylights have complex 3D structure.Wires go across the trees have high elevations and are intersected with each other.A large number of vehicles are parked under tree crowns or near low vegetation.After the step-wise point cloud segmentation performed, 10697 segments are extracted from Area 2. 414 segments are selected as training samples for the RF classifier.Details of the training samples in Area 2 are shown in Table 2. Specifically, the number of training samples for the class ground is smaller than other classes.Although a mountain exists in Area 2, a majority of ground regions are flat and smooth.Therefore, points on ground in Area 2 are clustered into a smaller number of segments than other classes.The third one is selected from the same dataset of Area 2. Area 3 is shown in Figure 6b and it contains 685870 points with an average point spacing of 1.0 m.In Figure 6b, the point cloud of Area 3 is colored by elevation.In this area, all the points are on a mountainous ground which is full of trees and ornamented by several houses.Area 3 is utilized to test the transplanting of the proposed classification framework to mountainous areas.There are only three types of objects, i.e., ground, buildings and vegetation.The ground of Area 3 is rugged with step edges.All the buildings are surrounded by vegetation.After the step-wise point cloud segmentation performed, 6658 segments are extracted from Area 3. 165 segments are selected as the training samples for the RF classifier.Details of the training samples in Area 3 are shown in Table 2. Specifically, the number of training samples for the class ground in Area 3 is larger than those in Area 1 and 2 because of the complex topographies.
In this paper, we classify Area 1 and 2 into five classes, i.e., ground, building, vegetation, vehicle and wire, and Area 3 into three classes, i.e., ground, building and vegetation.Area 1 and 2 are utilized to analyze the impact of the step-wise point cloud segmentation, the integration of RF and segment-based classification method, and the post-processing stage.Area 3 is employed to analyze to the classification result of a mountainous area with complex topographies.
To quantitatively analyze the classification accuracy, we obtain ground true for the three study areas by manual labelling.The details of the ground true for the three study areas are shown in Table 3.Note that quantitative analysis is derived by individual points rather than segments, because the hypothesis that no error exists in segments is unreasonable.Therefore, we present the ground true using individual points rather than segments.

Evaluation Metrics
For evaluation, we employ the confusion matrix and consider five commonly used measures: overall accuracy OA, Kappa coefficient KA. , precision P, recall R, and F 1 -score.They are computed according to the confusion matrix as follows: where tp i is the main diagonal element in i-th row, f p i is computed from the sum of i-th column, excluding the main diagonal element, f n i is the sum along i-th row, excluding the main diagonal element, m is the number of classes, and N is the number of all the points in an input point cloud.

Impact of the Step-Wise Point Cloud Segmentation
Our proposed step-wise point cloud segmentation method can extract three kinds of segments, i.e., planar, smooth and rough surfaces.We compare our method with the RG method [57], which has been published in PCL.The PCL supplies two kinds of processes based on the RG method.The first one is designed for plane extraction (RG-based plane segmentation) and the second one is designed for smooth surface extraction (RG-based smoothness segmentation).Six parameters in RG method published in PCL are used, i.e., the number k n of neighbors for normal estimation, the number k g of neighbors for growing, the smoothness threshold th θ , the curvature threshold th c , the residual threshold th r , and the minimum size threshold s min .

Qualitative Comparison
To visually compare our segmentation method with the RG methods, we select a small part from Area 1 to present the segmentation results.The parameter setting of our segmentation method is shown in Table 4 and the parameter setting of the RG method is shown in Table 5. Figure 7 presents the comparison of these segmentation results.We can find that our segmentation method not only can detects three kinds of segments, but also can clusters the ground points into a small number of segments, especially a single segment.Besides, intersecting planes, building roofs covered by vegetation, and small or narrow planes can be extracted by our segmentation method (see Figure 7a).The plane segmentation has the following shortcomings: (1) The segmentation result of the intersection between two planes is insufficient; (2) Small objects such as cars cannot be extracted; (3) Most majority of scattered points such as tree points cannot be segmented.The smoothness segmentation has the following shortcomings: (1) Intersected planes cannot be divided; (2) Small objects such as cars cannot be extracted, although it is implemented to extract smooth surfaces; (3) The building points and tree points could be clustered into the same region if they overlap with each other.

Quantitative Comparison
To quantitatively analyze the advantages of our segmentation method, the entire data of Area 1 and Area 2 are utilized.The time costs of our proposed segmentation method for processing Area 1 and 2 are 3.7 min and 4.4 min, respectively.We compare our method with the RG methods in terms of the classification measure score.In order to facilitate an objective comparison, all results based on the RG methods are averaged over 20 runs.The parameter setting of the step-wise point cloud segmentation method is shown in Table 6.The -score comparison is shown in Figure 8.The red bar describes the -score of our segmentation method, while the green bar describes the -score of the RG-based smoothness segmentation method, and the blue bar describes the -score of the RG-based plane segmentation.All classes' -score values of our segmentation method are larger than those of the RG methods for both Area 1 and 2. Specially, the -score values of the class vehicle derived by our method are significantly higher.The reason is that our segmentation method clusters all the points in a small-scale object into a few segments.

Quantitative Comparison
To quantitatively analyze the advantages of our segmentation method, the entire data of Area 1 and Area 2 are utilized.The time costs of our proposed segmentation method for processing Area 1 and 2 are 3.7 min and 4.4 min, respectively.We compare our method with the RG methods in terms of the classification measure F 1 score.In order to facilitate an objective comparison, all results based on the RG methods are averaged over 20 runs.The parameter setting of the step-wise point cloud segmentation method is shown in Table 6.The F 1 -score comparison is shown in Figure 8.The red bar describes the F 1 -score of our segmentation method, while the green bar describes the F 1 -score of the RG-based smoothness segmentation method, and the blue bar describes the F 1 -score of the RG-based plane segmentation.All classes' F 1 -score values of our segmentation method are larger than those of the RG methods for both Area 1 and 2. Specially, the F 1 -score values of the class vehicle derived by our method are significantly higher.The reason is that our segmentation method clusters all the points in a small-scale object into a few segments.

Parameter Tuning for Random Forests
The two parameters in the RF classifier, i.e., and are utilized not only in the classification procedure, but also in the feature selection task.In the feature selection task, is set to 4 according to the existing studies [19,38], which is a default setting and turned out to be a good choice of OOB rate [19].However, different studies have different settings of .Therefore, we have to unfold a test before feature selection for finding an appropriate value of .Besides, this test is also meaningful for supervised classification tasks in which the setting of is still stochastic [23].
To find an appropriate value of for a stable classification, we test the classification procedure with varying from 100 to 1000 and from 1 to 9. For evaluating the classification results, we utilize overall accuracy and Kappa coefficient to analyze the overall performance.Figure 9a,b show the variation tendency of overall accuracy and Figure 9c,d show the variation tendency of the Kappa coefficient.To make them more concise, we compute the mean and the standard deviation of the derived overall accuracy and kappa coefficient values under different settings of (see Figure 10a-d).The stability of the classification result will increase with an increasing value of .It is surprising that the overall accuracy and Kappa coefficient increase rapidly and the standard deviation decreases markedly with the augment of until it reaches 400.This case occurs in both Area 1 and Area 2. It may be the reason that 500 is the default value of in the R package for random forests [51].Therefore, setting to 400 is reasonable in our approach, especially for decreasing the computational burden.The two parameters in the RF classifier, i.e., N tree and N f eas are utilized not only in the classification procedure, but also in the feature selection task.In the feature selection task, N f eas is set to 4 according to the existing studies [19,38], which is a default setting and turned out to be a good choice of OOB rate [19].However, different studies have different settings of N tree .Therefore, we have to unfold a test before feature selection for finding an appropriate value of N tree .Besides, this test is also meaningful for supervised classification tasks in which the setting of N tree is still stochastic [23].
To find an appropriate value of N tree for a stable classification, we test the classification procedure with N tree varying from 100 to 1000 and N f eas from 1 to 9. For evaluating the classification results, we utilize overall accuracy and Kappa coefficient to analyze the overall performance.Figure 9a,b show the variation tendency of overall accuracy and Figure 9c,d show the variation tendency of the Kappa coefficient.To make them more concise, we compute the mean and the standard deviation of the derived overall accuracy and kappa coefficient values under different settings of N tree (see Figure 10a-d).The stability of the classification result will increase with an increasing value of N tree .It is surprising that the overall accuracy and Kappa coefficient increase rapidly and the standard deviation decreases markedly with the augment of N tree until it reaches 400.This case occurs in both Area 1 and Area 2. It may be the reason that 500 is the default value of N tree in the R package for random forests [51].Therefore, setting N tree to 400 is reasonable in our approach, especially for decreasing the computational burden.

Parameter Tuning for Random Forests
The two parameters in the RF classifier, i.e., and are utilized not only in the classification procedure, but also in the feature selection task.In the feature selection task, is set to 4 according to the existing studies [19,38], which is a default setting and turned out to be a good choice of OOB rate [19].However, different studies have different settings of .Therefore, we have to unfold a test before feature selection for finding an appropriate value of .Besides, this test is also meaningful for supervised classification tasks in which the setting of is still stochastic [23].
To find an appropriate value of for a stable classification, we test the classification procedure with varying from 100 to 1000 and from 1 to 9. For evaluating the classification results, we utilize overall accuracy and Kappa coefficient to analyze the overall performance.Figure 9a,b show the variation tendency of overall accuracy and Figure 9c,d show the variation tendency of the Kappa coefficient.To make them more concise, we compute the mean and the standard deviation of the derived overall accuracy and kappa coefficient values under different settings of (see Figure 10a-d).The stability of the classification result will increase with an increasing value of .It is surprising that the overall accuracy and Kappa coefficient increase rapidly and the standard deviation decreases markedly with the augment of until it reaches 400.This case occurs in both Area 1 and Area 2. It may be the reason that 500 is the default value of in the R package for random forests [51].Therefore, setting to 400 is reasonable in our approach, especially for decreasing the computational burden.

RF-Based Feature Selection
To select features, we utilize the backward elimination method and iteratively fit a RF with the aforementioned parameter setting.In each fitted RF, the importance of each feature is estimated.Figure 11 shows all the feature importance values of Area 1 and 2 which are estimated in the first time fitting.

RF-Based Feature Selection
To select features, we utilize the backward elimination method and iteratively fit a RF with the aforementioned parameter setting.In each fitted RF, the importance of each feature is estimated.Figure 11 shows all the feature importance values of Area 1 and 2 which are estimated in the first time fitting.

RF-Based Feature Selection
To select features, we utilize the backward elimination method and iteratively fit a RF with the aforementioned parameter setting.In each fitted RF, the importance of each feature is estimated.Figure 11 shows all the feature importance values of Area 1 and 2 which are estimated in the first time fitting.accuracy is 0.2215, and the minimum overall mean decrease permutation accuracy is 0.03991.When only 10 features are remained, the overall mean decrease permutation accuracy is 0.07588, which divides the range with a ratio no less than 8:2.However, the tendency of Figure 12b decreases rapidly from 10 to 11 remained features, therefore, we select 11 features for Area 2 as a compromised solution.According to the aforementioned analysis, the number 11 of independent features is reasonable for both Area 1 and 2.   The backward elimination method first fits a RF from the feature set S f eas , then the feature with the lowest importance value is eliminated from S f eas .The two steps are iterated until the number of features in S f eas is equal to 4. Subsequently, we select features based on the eliminating order.
To determine the number of independent features, we present the variation tendency of the overall mean decrease permutation accuracy OD k in Figure 12.For Area 1 (see Figure 12a), the maximum overall mean decrease permutation accuracy OD max is 0.1834, and the minimum overall mean decrease permutation accuracy OD min is 0.03426.When only 11 features are remained, the overall mean decrease permutation accuracy OD i is 0.05864, which divides the range I r with a ratio no less than 8:2.Besides, the tendency in Figure 12a decreases rapidly until the remained feature number reaches 11.For Area 2 (see Figure 12b), the maximum overall mean decrease permutation accuracy OD max is 0.2215, and the minimum overall mean decrease permutation accuracy OD min is 0.03991.When only 10 features are remained, the overall mean decrease permutation accuracy OD i is 0.07588, which divides the range I r with a ratio no less than 8:2.However, the tendency of Figure 12b decreases rapidly from 10 to 11 remained features, therefore, we select 11 features for Area 2 as a compromised solution.According to the aforementioned analysis, the number 11 of independent features is reasonable for both Area 1 and 2.
According to the aforementioned procedure, the selected features for Area 1 are slope O sl , lowest eigenvalue λ 3 , elevation difference H di f f , scattering λ s , change of curvature λ c , eigenentropy λ e , omnivariance λ o , anisotropy λ a , highest eigenvalue λ 1 , normalized elevation H n and area O ar .According to the aforementioned procedure, the selected features for Area 1 are slope , lowest eigenvalue , elevation difference , scattering , change of curvature , eigenentropy , omnivariance , anisotropy , highest eigenvalue , normalized elevation and area .The selected features for Area 2 are scattering , area , elevation difference , lowest eigenvalue , omnivariance , relative elevation , slope , anisotropy , medium eigenvalue , eigenentropy and normalized elevation .The selected features are listed in elimination orders for Area 1 and 2 respectively.

Robustness of the Integration of RF and the Segment-Based Classification
To test this robustness and the stability of the integration of RF and the segment-based classification, we add some noisy vectors to the computed feature set.It is worth noting that we do not record which feature vector is noisy, and therefore, we do not know which one is the noisy vector in the feature set when we train a RF and make prediction.Then, we analyze the overall accuracy and Kappa coefficient values of the classification results with different numbers of noisy vectors.
Besides, we implement the integration of SVM and the segment-based classification presented in [36] and analyze its robustness.It is worth noting that both the two integration methods have no complemental step for dealing with noises.To impartially compare the robustness of the two integrations, the feature selection procedure in our method is not performed, due to the fact that feature selection is not employed to the integration in [36].
As shown in Figure 13, when there is no noisy feature vector, the two integration methods get the classification results with similar overall accuracy and kappa coefficient values.When the noisy vector number is increasing, the values of the integration presented in the literature [36] decrease.Besides, our method obtains a mean standard deviation 0.0014 of overall accuracy values (0.0017 for Area 1 and 0.0012 for Area 2), and a mean standard deviation 0.0024 of kappa coefficient values (0.003 for Area 1 and 0.002 for Area 2), with the noisy feature vector number from 0 to 10.
In a close-up theoretical inspection, when RF classifier splits a subset of features in a node, it finds a feature with the maximum entropy or Gini decrease.In this case, the noisy features can be eliminated and the best feature corresponding to the maximum decrease is selected.Therefore, our method which integrates RF with the segment-based classification is more stable and robust.The selected features for Area 2 are scattering λ s , area O ar , elevation difference H di f f , lowest eigenvalue λ 3 , omnivariance λ o , relative elevation H r , slope O sl , anisotropy λ a , medium eigenvalue λ 2 , eigenentropy λ e and normalized elevation H n .The selected features are listed in elimination orders for Area 1 and 2 respectively.

Robustness of the Integration of RF and the Segment-Based Classification
To test this robustness and the stability of the integration of RF and the segment-based classification, we add some noisy vectors to the computed feature set.It is worth noting that we do not record which feature vector is noisy, and therefore, we do not know which one is the noisy vector in the feature set when we train a RF and make prediction.Then, we analyze the overall accuracy and Kappa coefficient values of the classification results with different numbers of noisy vectors.
Besides, we implement the integration of SVM and the segment-based classification presented in [36] and analyze its robustness.It is worth noting that both the two integration methods have no complemental step for dealing with noises.To impartially compare the robustness of the two integrations, the feature selection procedure in our method is not performed, due to the fact that feature selection is not employed to the integration in [36].
As shown in Figure 13, when there is no noisy feature vector, the two integration methods get the classification results with similar overall accuracy and kappa coefficient values.When the noisy vector number is increasing, the values of the integration presented in the literature [36] decrease.Besides, our method obtains a mean standard deviation 0.0014 of overall accuracy values (0.0017 for Area 1 and 0.0012 for Area 2), and a mean standard deviation 0.0024 of kappa coefficient values (0.003 for Area 1 and 0.002 for Area 2), with the noisy feature vector number from 0 to 10.
In a close-up theoretical inspection, when RF classifier splits a subset of features in a node, it finds a feature with the maximum entropy or Gini decrease.In this case, the noisy features can be eliminated and the best feature corresponding to the maximum decrease is selected.Therefore, our method which integrates RF with the segment-based classification is more stable and robust.

Impact of the Post-Processing Stage Based on the Semantic Rules
After the supervised classification, an initial classification result is obtained, in which misjudged segments may exist.In the post-processing stage, we first define several semantic rules for each class, and then utilize them to optimize the initial classification result.Note that the utilized semantic rules may be different for different input data, and the post-processing stage cannot correct all the misjudgments, however is an optimizing procedure.The rules for ground and the rule 2 for vegetation are not employed when we optimize Area 1 while others are employed.The rule 3 for ground and the rule 2 for vegetation are not employed when we optimize Area 2 while others are employed.Some thresholds are utilized to constrain the optimization which are list in Table 7.These threshold values are determined based on the natural form of an object and the point spacing of the input point cloud.To qualitatively analyze the impact of the post-processing stage, we select two small areas in Area 1 and 2 respectively.Figure 14 shows the initial and the final results of the small areas.In the initial classification results (see Figure 14a,c), a number of segments are mislabeled as the class vehicle or wire.For example, building elements such as chimneys are often mislabeled as vehicle, tree segments with small size are often mislabeled as the class vehicle or wire.In the post-processing stage,

Impact of the Post-Processing Stage Based on the Semantic Rules
After the supervised classification, an initial classification result is obtained, in which misjudged segments may exist.In the post-processing stage, we first define several semantic rules for each class, and then utilize them to optimize the initial classification result.Note that the utilized semantic rules may be different for different input data, and the post-processing stage cannot correct all the misjudgments, however is an optimizing procedure.The rules for ground and the rule 2 for vegetation are not employed when we optimize Area 1 while others are employed.The rule 3 for ground and the rule 2 for vegetation are not employed when we optimize Area 2 while others are employed.Some thresholds are utilized to constrain the optimization which are list in Table 7.These threshold values are determined based on the natural form of an object and the point spacing of the input point cloud.To qualitatively analyze the impact of the post-processing stage, we select two small areas in Area 1 and 2 respectively.Figure 14 shows the initial and the final results of the small areas.In the initial classification results (see Figure 14a,c), a number of segments are mislabeled as the class vehicle or wire.For example, building elements such as chimneys are often mislabeled as vehicle, tree segments with small size are often mislabeled as the class vehicle or wire.In the post-processing stage, the second rule for vehicles can deal with the initial false results where other classes of objects are mislabeled as vehicle, the rules for wires can deal with the initial false results where tree segments are mislabeled as wire.In the final results presented in Figure 14b,d

Classification Results and Accuracies
The final results of Area 1 and 2 are shown in Figure 16.The time costs of our proposed classification framework for processing Area 1 and 2 are 21 min and 25 min, respectively.The feature extraction procedure is the most time consuming step which costs 9 min and 12 min for Area 1 and 2, respectively.To quantitatively analyze the accuracies of classification, we compute the confusion matrix and the aforementioned five measures in Tables 8 and 9 based on the reference data generated by manual labelling.In addition, there may be missing points after the step-wise point cloud segmentation performed.These missing points will exist in the finial classification results and but are not labeled.In the quantitative evaluation, if some missing points belong to a class according to the ground true, they will be considered for evaluation.The revealed missing points are shown in the confusion matrices.Note that quantitative analysis is derived by individual points rather than segments, because the hypothesis that no error exists in segments is unreasonable.
As shown in Tables 8 and 9, the proposed method achieves a mean overall accuracy of 0.95255 (0.9697 and 0.9374 for Area 1 and 2 respectively), and a mean kappa coefficient of 0.9231 (0.9442 and For quantitative analysis, we compare each class's F 1 -score for the initial results and the final results.Figure 15a,b show the F 1 -score comparison for Area 1 and Area 2. We can find that the F 1 -score of the class vehicle is improved obviously both in Area 1 and 2. Besides, the F 1 -score values of other classes are also improved at different levels.

Classification Results and Accuracies
The final results of Area 1 and 2 are shown in Figure 16.The time costs of our proposed classification framework for processing Area 1 and 2 are 21 min and 25 min, respectively.The feature extraction procedure is the most time consuming step which costs 9 min and 12 min for Area 1 and 2, respectively.To quantitatively analyze the accuracies of classification, we compute the confusion matrix and the aforementioned five measures in Tables 8 and 9 based on the reference data generated by manual labelling.In addition, there may be missing points after the step-wise point cloud segmentation performed.These missing points will exist in the finial classification results and but are not labeled.In the quantitative evaluation, if some missing points belong to a class according to the ground true, they will be considered for evaluation.The revealed missing points are shown in the confusion matrices.Note that quantitative analysis is derived by individual points rather than segments, because the hypothesis that no error exists in segments is unreasonable.
As shown in Tables 8 and 9, the proposed method achieves a mean overall accuracy of 0.95255 (0.9697 and 0.9374 for Area 1 and 2 respectively), and a mean kappa coefficient of 0.9231 (0.9442 and

Classification Results and Accuracies
The final results of Area 1 and 2 are shown in Figure 16.The time costs of our proposed classification framework for processing Area 1 and 2 are 21 min and 25 min, respectively.The feature extraction procedure is the most time consuming step which costs 9 min and 12 min for Area 1 and 2, respectively.
To quantitatively analyze the accuracies of classification, we compute the confusion matrix and the aforementioned five measures in Tables 8 and 9 based on the reference data generated by manual labelling.In addition, there may be missing points after the step-wise point cloud segmentation performed.These missing points will exist in the finial classification results and but are not labeled.In the quantitative evaluation, if some missing points belong to a class according to the ground true, they will be considered for evaluation.The revealed missing points are shown in the confusion matrices.Note that quantitative analysis is derived by individual points rather than segments, because the hypothesis that no error exists in segments is unreasonable.
As shown in Tables 8 and 9, the proposed method achieves a mean overall accuracy of 0.95255 (0.9697 and 0.9374 for Area 1 and 2 respectively), and a mean kappa coefficient of 0.9231 (0.9442 and 0.9020 for Area 1 and 2 respectively).The accuracy values of classes ground, building and vegetation are rather good.The accuracy values of vehicle and wire, are relatively lower than those of the other classes.The class vehicle achieves a mean F 1 -score of 0.8504 (0.8519 and 0.8488 for Area 1 and 2, respectively), and the class wire achieves a mean F 1 -score of 0.7531 (0.8155 and 0.6907 for Area 1 and 2, respectively).
The classification of small objects such as vehicles is the most challenging task.Generally, vegetation points make a 3D urban scene more complex and affect the classification of small objects.For example, wires often go across the trees with high elevation and vehicles are often parked near low vegetation or under tree canopies.Therefore, in the confusion matrices, vegetation points are often confused with another class points, and thus decrease the classification accuracy of other classes, especially for the classes of vehicle and wire.However, the accuracy of class vegetation achieves a mean F 1 -score of 0.9422 (0.9429 and 0.9415 for Area 1 and 2, respectively) which is superior.The classification of small objects such as vehicles is the most challenging task.Generally, vegetation points make a 3D urban scene more complex and affect the classification of small objects.For example, wires often go across the trees with high elevation and vehicles are often parked near low vegetation or under tree canopies.Therefore, in the confusion matrices, vegetation points are often confused with another class points, and thus decrease the classification accuracy of other classes, especially for the classes of vehicle and wire.However, the accuracy of class vegetation achieves a mean -score of 0.9422 (0.9429 and 0.9415 for Area 1 and 2, respectively) which is superior.

Transplantation to Mountainous Areas
To test the transplanting of the proposed method, extracted features and obtained parameters to mountainous areas, we classify a different 3D scene represented by Area 3. Area 3 is a mountainous area with step edges and rock faces.First, Area 3 is segmented by the proposed step-wise point cloud segmentation method with the same parameters as Area 2, and segment features are extracted.Next, a small feature set is selected by the obtained parameters of RF, i.e., = 400, and = 4. 9 features are selected by the backward elimination method presented in Section 3.4.2,and they are lowest eigenvalue , scattering slope , change of curvature , planarity , anisotropy , eigenentropy , tangent plane projection area .Then, the ALS point cloud of Area 3 is classified by the select feature set.None of the semantic rules are utilized in the post-processing stage.16 min are costed during the whole procedures among which the feature extraction procedure is the most time consuming step.The confusion matrix contained recall, precision, and -score values is shown in Table 10, and its corresponding segmentation and classification results are shown in Figure 17.

Transplantation to Mountainous Areas
To test the transplanting of the proposed method, extracted features and obtained parameters to mountainous areas, we classify a different 3D scene represented by Area 3. Area 3 is a mountainous area with step edges and rock faces.First, Area 3 is segmented by the proposed step-wise point cloud segmentation method with the same parameters as Area 2, and segment features are extracted.Next, a small feature set is selected by the obtained parameters of RF, i.e., N tree = 400, and N f eas = 4. 9 features are selected by the backward elimination method presented in Section 3.  10, and its corresponding segmentation and classification results are shown in Figure 17.As shown in Table 10, the proposed classification framework and obtained parameters of RF achieve an overall accuracy of 0.9117, and a Kappa coefficient of 0.8379.The class ground achieves a -score of 0.9026 which is lower than Area 1 and 2, because the ground of Area 3 is mountainous and rugged with step edges.The class building achieves a -score of 0.8760, and the class vegetation achieves a -score of 0.9327.

Uncertainties, Errors and Accuracies
The above experiments suggest that our proposed method obtains good results.There are three improvements in our proposed classification framework, which improve the accuracies of the classification results.However, there are still some errors in the classification results.We will list them according to missing points and the five aforementioned classes.
In the confusion matrices, the missing points appear more likely in the classes building and vegetation.For the class vegetation, the laser beam may penetrate the tree surface and collects a point in an internal branch.The internal point may be an isolated point, however, belong to the class vegetation.For the class building, the missing points often appear on building facades duo to the fact that building facades are incomplete and points in them are sparse in large-scale ALS point clouds.
For the class ground, uncertainties and errors are more likely to arise in the areas with mountainous and rugged topography.In Table 8, the class ground achieves a -score of 0.9939 due to the fact that the ground in Area 1 is flat and smooth.In Table 9, the -score of ground (0.9783) is lower than that of Area 1, because there is a hill that causes the topographical complexity increasing.As shown in Table 10, the proposed classification framework and obtained parameters of RF achieve an overall accuracy of 0.9117, and a Kappa coefficient of 0.8379.The class ground achieves a F 1 -score of 0.9026 which is lower than Area 1 and 2, because the ground of Area 3 is mountainous and rugged with step edges.The class building achieves a F 1 -score of 0.8760, and the class vegetation achieves a F 1 -score of 0.9327.

Uncertainties, Errors and Accuracies
The above experiments suggest that our proposed method obtains good results.There are three improvements in our proposed classification framework, which improve the accuracies of the classification results.However, there are still some errors in the classification results.We will list them according to missing points and the five aforementioned classes.
In the confusion matrices, the missing points appear more likely in the classes building and vegetation.For the class vegetation, the laser beam may penetrate the tree surface and collects a point in an internal branch.The internal point may be an isolated point, however, belong to the class vegetation.For the class building, the missing points often appear on building facades duo to the fact that building facades are incomplete and points in them are sparse in large-scale ALS point clouds.
For the class ground, uncertainties and errors are more likely to arise in the areas with mountainous and rugged topography.In Table 8, the class ground achieves a F 1 -score of 0.9939 due to the fact that the ground in Area 1 is flat and smooth.In Table 9, the F 1 -score of ground (0.9783) is lower than that of Area 1, because there is a hill that causes the topographical complexity increasing.In Area 3, all the ground points are located at a mountainous region which are rugged and full of step edges and rock faces.The F 1 -score of class ground in Area 3 is 0.9026 which is much lower than those in Area 1 and 2. However, according to the recall and precision values of ground in Area 3, our proposed classification framework achieves the average accuracy of existing ALS point cloud filtering methods.
For the class building, the F 1 -score (0.9368) in Area 1 is superior to those in Area 2 and 3. A close-up visual inspection shows that there are more buildings with vegetation confused in Area 2 and 3. Especially in Area 3, all the buildings are surrounded by vegetation and most of them are overlapped by tall trees.Although the regions where buildings and vegetation are confused are segmented correctly by the step-wise point cloud segmentation method, these regions still decrease the classification accuracy of the class building.
For the class vegetation in confusion matrix, vegetation points are often confused with other classes points.A close-up visual inspection shows that man-made objects such as buildings, vehicles and wires are often near vegetation.For example, buildings tend to be surrounded by trees, wires often go across tall trees.Therefore, vegetation makes a 3D scene more complex and affects the classification accuracy of other classes.However, the class vegetation achieves a mean F 1 -score of 0.9390 (0.9429, 0.9415 and 0.9327 for Area 1, 2 and 3, respectively) which is superior due to the large cardinal number.
For the class vehicle, the mean F 1 -score is 0.8504 (0.8519 and 0.8488 for Area 1 and 2, respectively) which is lower than the classes ground, building, and vegetation.Most misjudgments are caused by vegetation according to the confusion matrices, though the scene where vehicles are parked close to low vegetation is able to be correctly segmented.A close-up visual inspection shows that the mislabeled points always distribute randomly and irregularly.In an extreme case, it is inevitable that misjudgments exist in low vegetation which is in the same geometric form as a vehicle.However, the classification accuracy of the class vehicle is superior to the existing studies according to the analysis presented in Section 3.2.
For the class wire, the mean F 1 -score is 0.7531 (0.8155 and 0.6907 for Area 1 and 2, respectively).The wire is the class with the lowest accuracy in our experiments.The wires in our study areas are rather common low voltage electrical wires than some special parts of overhead electric power transmission corridors.The misjudgments arise in the areas where wires go across trees.In these areas, wires and trees have similar point densities and cannot be divided by the segmentation method.Therefore, there will be a number of mislabeled points in the classification results of these areas, which affect the accuracy of the class wire most.

Conclusions
In this paper, we classify ALS point clouds via a framework with four stages, i.e., (i) step-wise point cloud segmentation; (ii) feature extraction; (iii) RF-based feature selection and classification; (iv) post-processing.In the first stage, the step-wise point cloud segmentation extracts three kinds of segments, i.e., planar, smooth and rough surfaces.Planar and smooth surfaces are more easy to characterize piecewise planar objects, and rough surfaces are more easy to characterize objects with irregular shapes.In the second stage, we extract geometric features from the input ALS point cloud by considering segments as the basic computational units.In the third stage, we integrate RF with the segment-based classification method to classify ALS point clouds.Discriminative features are selected using the backward elimination method based on OOB errors, and an appropriate value (400) for the number of trees used in RF is determined.At last, we employ semantic information to define several rules for each class, and utilize them in the post-processing stage to optimize the classification results.
There are two contributions in the framework, i.e., step-wise point cloud segmentation, and the integration of RF and the segment-based classification method.In the step-wise point cloud segmentation, we utilize a RANSAC method to optimize the normal vector and neighborhood of each point, and next grow a region among the optimized neighborhood for each seed point.Then, scattered points are detected and initial patches are constructed.Finally, the log Euclidean Riemannian metric is utilized as a constraint to connect the initial patches to rough surfaces.Experiments validate that the step-wise segmentation is good at recognizing small-scale objects.To analyze the integration of RF and the segment-based classification, we first find a suitable parameter setting of RF, then select features based on these parameters, and finally analyze the robustness and show the benefits of the integration.
There is a limit existing in our method, i.e., objects with less geometric distinguishability cannot be recognized, such as flat roads which has similar geometric attributes with ground.In future work, we will take more features and complemental strategies into consideration to classify these kinds of objects.

Figure 1 .
Figure 1.The proposed framework for airborne laser scanning (ALS) point cloud classification.

Figure 1 .
Figure 1.The proposed framework for airborne laser scanning (ALS) point cloud classification.

Figure 2 .
Figure 2. The flowchart of the step-wise point cloud segmentation.

Figure 2 .
Figure 2. The flowchart of the step-wise point cloud segmentation.

Figure 3 .
Figure 3. Results of internal steps in the step-wise point cloud segmentation: (a) Intermediate result after the process of the region growing (RG) with RANdom SAmple Consensus (RANSAC) step.The scattered points are rendered in black; (b) Intermediate result after the process of the initial patch construction; (c) The result of the step-wise point cloud segmentation.

Figure 3 .
Figure 3. Results of internal steps in the step-wise point cloud segmentation: (a) Intermediate result after the process of the region growing (RG) with RANdom SAmple Consensus (RANSAC) step.The scattered points are rendered in black; (b) Intermediate result after the process of the initial patch construction; (c) The result of the step-wise point cloud segmentation.

Figure 6 .
Figure 6.The testing data, (a) shows the point cloud of Area 1, (b) shows the point clouds for Area 2 and 3.

Figure 6 .
Figure 6.The testing data, (a) shows the point cloud of Area 1, (b) shows the point clouds for Area 2 and 3.
between 28 September and 23 November 2013 by the Watershed Sciences, Inc. (WSI).The dataset is provided by the University of Maryland and the Sonoma Country Vegetation Mapping and Lidar Program under grant NNX13AP69G from NASA's Carbon Monitoring System.Area 2 is shown in Figure

Figure 7 .
Figure 7.Comparison of our segmentation method with the RG methods presented in Point Cloud Library (PCL): (a) The result of our method; (b) The result of the RG-based plane segmentation; (c) The result of the RG-based smoothness segmentation.

Figure 7 .
Figure 7.Comparison of our segmentation method with the RG methods presented in Point Cloud Library (PCL): (a) The result of our method; (b) The result of the RG-based plane segmentation; (c) The result of the RG-based smoothness segmentation.

Figure 8 .
Figure 8.All classes' -score values of different segmentation methods for Area 1 and Area 2: (a)score values of Area 1; (b) -score values of Area 2.

Figure 8 .
Figure 8.All classes' F 1 -score values of different segmentation methods for Area 1 and Area 2: (a) F 1 -score values of Area 1; (b) F 1 -score values of Area 2.

Figure 8 .
Figure 8.All classes' -score values of different segmentation methods for Area 1 and Area 2: (a)score values of Area 1; (b) -score values of Area 2.

Figure 9 .Figure 9 .Figure 10 .
Figure 9.The variation tendencies of overall accuracy and Kappa coefficient under different settings of N tree : (a,b) are for Area 1 and Area 2 respectively; (c,d) are for Area 1 and Area 2 respectively.

Figure 10 .
Figure 10.The variation tendencies of mean overall accuracy values and standard deviation values of overall accuracy and Kappa coefficient under different settings of N tree ; (a,b) are for Area 1 and Area 2 respectively; (c,d) are for Area 1 and Area 2 respectively.

Figure 11 .
Figure 11.Feature importance values based on the mean decrease permutation accuracy estimated by the first time fitted RF: (a) is for Area 1; (b) is for Area 2.Figure 11.Feature importance values based on the mean decrease permutation accuracy estimated by the first time fitted RF: (a) is for Area 1; (b) is for Area 2.

Figure 11 .
Figure 11.Feature importance values based on the mean decrease permutation accuracy estimated by the first time fitted RF: (a) is for Area 1; (b) is for Area 2.Figure 11.Feature importance values based on the mean decrease permutation accuracy estimated by the first time fitted RF: (a) is for Area 1; (b) is for Area 2.

Figure 12 .
Figure 12.Feature selection via iterative elimination for Area 1 and 2 based on the overall mean decrease permutation accuracy: (a) is for Area 1; (b) is for Area 2.

Figure 12 .
Figure 12.Feature selection via iterative elimination for Area 1 and 2 based on the overall mean decrease permutation accuracy: (a) is for Area 1; (b) is for Area 2.

Figure 13 .
Figure 13.Comparison between the two integrations under different numbers of noisy feature vectors, (a,b) depict the accuracy values of our method for Area 1 and 2 respectively, (c,d) depict the accuracy values of the literature [36] for Area 1 and 2 respectively.

Figure 13 .
Figure 13.Comparison between the two integrations under different numbers of noisy feature vectors, (a,b) depict the accuracy values of our method for Area 1 and 2 respectively, (c,d) depict the accuracy values of the literature [36] for Area 1 and 2 respectively.

Figure 14 .
Figure 14.Analysis of the impact of the post-processing stage based on semantic rules: (a) The initial result of the first small area; (b) The final result of the first small area: (c) The initial result of the second small area; (d) The final result of the second small area.

Figure 14 .
Figure 14.Analysis of the impact of the post-processing stage based on semantic rules: (a) The initial result of the first small area; (b) The final result of the first small area: (c) The initial result of the second small area; (d) The final result of the second small area.

Figure 14 .
Figure 14.Analysis of the impact of the post-processing stage based on semantic rules: (a) The initial result of the first small area; (b) The final result of the first small area: (c) The initial result of the second small area; (d) The final result of the second small area.

Figure 15 .
Figure 15.F 1 comparison between the initial and final results of Area 1 and 2: (a) The comparison of Area 1; (b) The comparison of Area 2.

Figure 16 .
Figure 16.The results of the proposed classification framework for areas Area 1 and 2, (a) is the segmentation result of Area 1; (b) is the classification result of Area 1; (c) is the segmentation result of Area 2; (d) is the classification result of Area 2.
3.2, and they are lowest eigenvalue λ 3 , scattering λ s , omnivariance λ o , slope O sl , change of curvature λ c , planarity λ p , anisotropy λ a , eigenentropy λ e , tangent plane projection area PA t .Then, the ALS point cloud of Area 3 is classified by the select feature set.None of the semantic rules are utilized in the post-processing stage.16 min are costed during the whole procedures among which the feature extraction procedure is the most time consuming step.The confusion matrix contained recall, precision, and F 1 -score values is shown in Table

Figure 17 .
Figure 17.The results of the proposed classification framework for Area 3, (a) is the segmentation result of Area 3; (b) is the classification result of Area 3.

Figure 17 .
Figure 17.The results of the proposed classification framework for Area 3, (a) is the segmentation result of Area 3; (b) is the classification result of Area 3.

Algorithm 1. Normal estimation. Input: Point cloud={P}. Parameters: k
n , d r 1: Regular point set {RP} ← ∅ , Scattered point set {SP} ← ∅ , Inlier Set {IN} ← ∅ , Proportion Set {PR} ← ∅ ,Specifically, when plane-based RANSAC algorithm fits a local plane, it divides k n neighbors into inliers PI i and outliers (Row 5).The points in PI i are on the fitted plane.If the query point p i is in PI i , we compute a proportion pr i of the inliers PI i to the k n neighbors (Row 7).

•
Considering an extreme case, we assume that some vehicle segments are misjudged as low vegetation segments.First, a misjudged query segment should have a neighbor labeled as ground, then, if the height, width and length of the segment are in certain

Table 2 .
Number of training samples (segments) per class for the three study areas.

Table 3 .
Number of points per class in ground true of the three study areas.

Table 4 .
The parameter setting of our method for a small part in Area 1.

Table 5 .
The parameter setting of RG method for a small part in Area 1.

Table 6 .
The parameter setting of our method for Area 1 and 2.

Table 6 .
The parameter setting of our method for Area 1 and 2.

Table 7 .
The threshold values for the post-processing stage used for Area 1 and 2.

Table 7 .
The threshold values for the post-processing stage used for Area 1 and 2.

Table 8 .
The analysis of Area 1.

Table 9 .
The accuracy analysis of Area 2. Area 1 and 2 respectively).The accuracy values of classes ground, building and vegetation are rather good.The accuracy values of vehicle and wire, are relatively lower than those of the other classes.The class vehicle achieves a mean -score of 0.8504 (0.8519 and 0.8488 for Area 1 and 2, respectively), and the class wire achieves a mean -score of 0.7531 (0.8155 and 0.6907 for Area 1 and 2, respectively).

Table 8 .
The accuracy analysis of Area 1.

Table 9 .
The accuracy analysis of Area 2.

Table 10 .
The accuracy analysis of Area 3.

Table 10 .
The accuracy analysis of Area 3.