Major Orientation Estimation-Based Rock Surface Extraction for 3 D Rock-Mass Point Clouds

In the fields of 3D modeling, analysis of discontinuities and engineering calculation, surface extraction is of great importance. The rapid development of photogrammetry and Light Detection and Ranging (LiDAR) technology facilitates the study of surface extraction. Automatic extraction of rock surfaces from 3D rock-mass point clouds also becomes the basis of 3D modeling and engineering calculation of rock mass. This paper presents an automated and effective method for extracting rock surfaces from unorganized rock-mass point clouds. This method consists of three stages: (i) clustering based on voxels; (ii) estimating major orientations based on Gaussian Kernel and (iii) rock surface extraction. Firstly, the two-level spatial grid is used for fast voxelization and segmenting the point cloud into three types of voxels, including coplanar, non-coplanar and sparse voxels. Secondly, the coplanar voxels, rather than the scattered points, are employed to estimate major orientations by using a bivariate Gaussian Kernel. Finally, the seed voxels are selected on the basis of major orientations and the region growing method based on voxels is applied to extract rock surfaces, resulting in sets of surface clusters. The sub-surfaces of each cluster are coplanar or parallel. In this paper, artificial icosahedron point cloud and natural rock-mass point clouds are used for testing the proposed method, respectively. The experimental results show that, the proposed method can effectively and accurately extract rock surfaces in unorganized rock-mass point clouds.


Introduction
Feature surface extraction in point clouds is a requisite part of many computer graphics, image processing and computer vision, including 3D reconstruction [1], object recognition [2,3], virtual reality [4], etc. Remote sensors, such as Light Detection and Ranging (LiDAR) and Differential SAR Interferometry (DInSAR), can capture high resolution and accurate 3D information of object surface from a considerable distance [5,6].The rapid development and popularization of remote sensors promotes the research of surface extraction and 3D modeling in unorganized point clouds.Nevertheless, the existing surface-extraction methods can't deal with the growing size of the point clouds.
Currently, the three common methods of surface extraction in point clouds are Hough Transform (HT) [7], Random Sample Consensus (RANSAC) [8] and Region Growing (RG) [9,10].The Hough Transform is a traditional approach of detecting parameterized shapes and objects in 2D images.It is typically used for lines and circles detection [11].The 3D Hough Transform can also extract feature surfaces, cylinders and spheres in 3D point clouds [12].However, the traditional 3D Hough Transform is complex, time consuming and the detected surfaces may be disconnected.To detect a surface with RANSAC, three points are chosen randomly to calculate the surface parameters defined by them.A score function is used to determine whether the model fits the remaining points [13].
The algorithm stops when it reaches stability.RANSAC may detect spurious feature surfaces or discontinuous surfaces.Region Growing performs a local search to identify and expand the regions with same characteristics [14].Region Growing may lead to holes and over-segmentation issues.
In recent years, the prevailing methods of surface extraction have focused on urban scenes or indoor scenes.While the urban point clouds consist of regular geometric surfaces (such as walls and floors, etc.), the rock-mass point clouds are shaggy and mostly in the form of high and steep slope [15].Stability analysis is the core of geotechnical engineering such as aboideau, tunnel and highwall [5,16,17].Presently, the various methods of stability analysis in rock mass (such as block theory [18], discontinuous deformation analysis [19] and numerical manifold method [20], etc.) mostly be performed on the accurate numerical models of rock mass.The rock surface extraction and modeling in rock-mass point clouds becomes a requisite part of geotechnical engineering.The main driver of our methodology is to extract rock surfaces and reconstruct the watertight 3D model of rock mass itself.
In this paper, the proposed method is designed based on the rough surface and irregular shape of rock mass.By using the spatial grid for fast voxelization, the approximate coplanar points are classified.At the same time, the different coplanar criteria are compared and advice for adjusting the parameters of coplanar criteria are pointed out.To estimate major orientations and overcome the shortage of computational complexity, we cast votes for the normal vectors [21] of every approximate coplanar voxels by using a bivariate Gaussian Kernel [22].Eventually, according to the results of major orientation estimation, the voxels of region growing units are selected [23].The region growing based on voxels guarantees the connectivity of detected rock surfaces and improves efficiency, which results in sets of coplanar clusters and the largest connected sub-surfaces of each cluster.Notice that the discontinuity sets are now actual planes but not surfaces (such as folds, etc.).To its analysis, it is common to accept them as planes, but for that consideration, it must be assumed that the region of analysis has homogeneous characteristics.
The remaining part of this paper is structured as follows: after providing an overview on related work, the proposed method of rock surface extraction is thoroughly discussed in Section 3. Section 4 shows and analyzes the experimental results.In the final Section 5, we conclude this paper and point out the future work.

Hough Transform
The 2D Hough Transform is a standard method proposed by Hough [7] which could detect geometric primitives in the field of digital image processing.The 3D Hough Transform is proposed on the basis of 2D Hough Transform, which is able to detect the feature surfaces in 3D point clouds.Every single point in point set {P} is converted to a sinusoidal surface in Hough space and voting for the accumulator cells which intersect with the sinusoidal surface (i.e., Hough Voting).After voting for all points of {P}, the local maximum values are selected from the Hough space (i.e., peak detection).Then the peaks are converted to the surface equations and the sample points that satisfy these equations are found out.The time complexity of Hough Transform is O(|P|N θ N φ ), where |P| is the number of sample points, N θ is the number of accumulator cells in direction of θ and N φ in direction of φ, respectively.The time complexity and space complexity of HT depend on the discretization degree of the accumulator.The higher the discretization degree, the higher computational complexity.
Borrmann et al. [24] summarized the traditional Hough Transform and designed a new spherical accumulator.Different variants of the Hough Transform were compared and analyzed through experiments, including the standard Hough Transform (SHT), probability Hough Transform (PHT), random Hough Transform (RHT) and so on.Vosselman et al. [21] proposed a two-step procedure for Hough Transform, which used normal vectors to confirm the spatial position of the parameter.The two-step procedure reduces the complexity of mapping from the measurement space to the parameter space.Huang and Brenner [25] proposed a joint multi-plane detection strategy to improve the performance of the Hough Transform, which had good performance on the roof with ridge lines.Hulik et al. [26] proposed an optimization method for continuous surface extraction in point clouds.The integral image normal computation of the Point Cloud Library (PCL) [27] was used for estimating the local normal vectors of each point.When accumulating to the Hough space, a hierarchical structure was used to ensure high-resolution would not consume additional memory.A buffer strategy was proposed for Gaussian smoothing and a sliding-window technique was applied in Peak detection.Fernandes and Oliveira [28] proposed the Kernel-based Hough Transform (KHT) for real-time detection of straight lines in two-dimensional images.Due to the lack of clear information of the adjacent samples, it was not applicable to unorganized point clouds.In this regard, Limberger and Oliveira [22] proposed an improved Octree-based clustering strategy for 3D unorganized point clouds (i.e., 3D-KHT).Octree and Principal Component Analysis (PCA) were used to segment approximate coplanar samples.Hough voting for each sample point was carried out on the basis of Gaussian Kernel, which significantly increased the efficiency.

Region Growing
The use of region growing in point clouds to extract structural feature surfaces is a 3D simulation of region growing in the field of Image Processing [23,29].By selecting a seed, it searches for local to identify and expand regions with the same characteristics.
Poppinga et al. [30] used priority queues and depth images to quickly extract the adjacent neighbors.The calculation methods of incremental covariance matrices and mean square error (MSE) through mathematical formulas were introduced.When a new point was added, there was no need to construct a new covariance matrix, which improved the efficiency of region growing.Xiao et al. [23] proposed the strategy of using sub-windows as seeds for region growing.A sub-window-based region growing (SBRG) for structured environment and a hybrid region growing (HRG) which combined sub-windows and spatial points for unstructured environment were introduced respectively.Vo et al. [31] proposed an octree-based region growing method for urban point cloud segmentation.The algorithm was divided into two stages based on a coarse-to-fine concept.In the first stage, the octree-based region growing was performed for major segmentation.The second stage was a refinement process.This method was at least an order of magnitude faster than conventional region growing algorithms.

Random Sample Consensus
Another important parameterization method for surface extraction in point clouds is Random Sample Consensus (RANSAC) [8].It iteratively selects three points to calculate the surface equation and count the number of points located in this surface.The time complexity of RANSAC is O(I|P|), where I is the number of iterations.RANSAC is robust to noise.However, due to its randomness, the time complexity is uncertain.
Yang and Förstner [32] presented an approach of surface extraction by integrating RANSAC and Minimum Description Length (MDL).This method could avoid detecting false feature surfaces caused by the complex geometry of point cloud.To avoid generating spurious-surfaces (i.e., consisting of points from several coplanar surfaces), Xu et al. [33] proposed the weighted RANSAC method according to the difference in error distribution between proper and improper surface hypotheses.Li et al. [34] proposed an improved RANSAC method based on Normal Distribution Transformation (NDT) cells.In each iteration, a surface NDT cell was selected as the minimum sample to ensure the correctness of the sampling on the same surface.

Other Methods
Holz et al. [35] presented a method to cluster the sample points in normal space and spherical coordinates simultaneously.By using 3D semantic awareness, this method could conduct obstacle detection and collision avoidance.Yamazaki et al. [36] utilized topology method and graph-cut method to extract surfaces.Without constructing a triangular mesh or others, the surfaces were segmented into different feature representations.
In the field of rock-mass engineering, some automatic or semi-automatic tools for surface extraction and analysis of discontinuities are publicly available to engineers and researchers (i.e., COLTOP 3D [37], Split-FX [38], DiAna [39], PlaneDetect [15], DSE [6] and so on).Riquelme et al. [6] presented a method for automatic characterization of rock mass discontinuities.This method analyzed the coplanarity of adjacent points to identify and define the algebraic equations of different rock surfaces.The Kernel Density Estimation (KDE) was used to find the main directions and the density-based scanning principle was used to identify the clusters.The authors published the Discontinuity Set Extractor (DSE) open-source software and this software also calculated the normal spacing and the discontinuity persistence [40,41].By considering the Hough Transform and region growing simultaneously, Leng et al. [42] proposed a multi-scale surface-detection algorithm.Hough Transform was used to extract the seeds and feature surfaces of rock mass were extracted by region growing.The rock surfaces of different spatial scales were detected separately by introducing the concept of multi-scale.

Methodology
On the basis of summarizing the conventional surface extraction methods in point clouds, this paper proposes an effective rock surface extraction method in rock-mass point clouds.The summary of this method is shown in Algorithm 1.In this method, the spatial grid is used for fast voxelization and clustering approximate coplanar points.After voting for each coplanar voxel in normal space, the major orientations are obtained and the region growing based on voxels is performed.Finally, multiple sets of rock surfaces with the largest connected sub-surfaces are obtained.The proposed method avoids the phenomenon of "riband" [42] and facilitates subsequent operations (i.e., surface merging, concave calculation and 3D reconstruction).

Clustering Based On Voxels
In the preprocessing stage, we use the kd-tree module and the k-Nearest Neighbor (knn) algorithm of the Point Cloud Library [27] to create indexes and estimate the normal vector of each point.It takes a bit of time to load and pre-process raw datasets.In recent years, researchers have proposed a newer version of neighborhood algorithms with better performance [43,44].Similar to octree, the two-scale spatial grid is a data structure for describing the three-dimensional space.As shown in Figure 1, each node of the spatial grid represents a cube voxel and could be divided into 8 sub-cubes (i.e., sub-voxels).The main reason we choose two-scale spatial grid instead of octree is that we could rapidly and easily searching for adjacent voxels by using three-dimensional array.The pseudo-code of voxelization and clustering is described in Algorithm 2. The coplanar clustering is done by analysing the eigenvalues of covariance matrix.For a voxel A, the covariance matrix ∑ ∑ ∑ is computed and coplanar constraints are used to determine whether the sample points contained in A have the coplanar characteristics.We sort the eigenvalues of matrix ∑ ∑ ∑ in ascending order (i.e., λ 1 ≤ λ 2 ≤ λ 3 ).With respect to the setting of the approximate coplanar constraints, Limberger and Oliveira [22] proposed a criterion for comparing λ (see Equation ( 1)) and Riquelme et al. [6] used the curvature to determine the coplanarity (see Equation ( 2)), where S α , S β and η are thresholds. (1) In this paper, the coplanar constraints shown in Equation (3) are used, among which S α and ε are thresholds.Mean Squared Error (MSE), as shown in Equation ( 4), is a statistical method to measure the average error and the variation of data, which is commonly used for region growing.
After coplanar clustering, we have obtained hundreds (thousands) of approximately coplanar voxels, and the non-coplanar phenomenon still exists in coplanar voxels.On the one hand, since the coplanar voxels contain non-coplanar sample points, the normal vectors (i.e., "principal directions") obtained by PCA for each "coplanar voxel" is not accurate.On the other hand, it is difficult to find the "major orientations" from the hundreds (thousands) of "principal directions" for selecting growing seeds.Taking these two aspects into consideration, we use a bivariate Gaussian Kernel to vote for the "principal direction" of each coplanar voxels (see Section 3.2).In the end, we selected a few more representative "major orientations" from the hundreds of "principal directions".

Major Orientation Estimation Based On Gaussian Kernel
In the above, the point cloud of rock mass is segmented by coplanar clustering, and the approximate coplanar voxels are obtained.In this section, major orientation estimation based on Gaussian Kernel is mainly performed on the normal vectors of these coplanar voxels.

Hemispherical Surface Accumulator
Many geologists and engineers have proposed various plane-detection methods through using stereonets.These methods mainly concentrate on the analysis of rock-mass discontinuities.Since the main driver of our methodology is to reconstruct the watertight 3D model of rock mass itself, we do not use stereonets here.If we map all unit normal vectors of coplanar voxels to a sphere, the end point of each unit normal vector must fall on the spherical surface.So we consider using the accumulator ball as the accumulator design.Notice that only the surface of the ball is divided into accumulator cells and all of these accumulator cells are of the same size [24] (as shown in Figure 2a).As we all known, the direction of normal vector can be positive or negative.If we reverse the direction of these normal vectors that fall in the lower hemisphere (i.e., z < 0 in Figure 2a), the end point of each normal vector must fall in the upper hemisphere(i.e., z ≥ 0 in Figure 2a).According to these features of normal vectors, we eventually decide to use the hemispherical surface accumulator.The example of voting result is shown in Figure 2b.When it is close to the poles (i.e., φ is close to 0), the voting value is smaller than other regions, which is similar to the spherical accumulator.This ultimately results in the vertical surfaces being detected earlier than the horizontal surfaces, which does not affect the correctness of the results [22,24].

Computing Bivariate Gaussian Kernels
Assuming that K is an approximate coplanar voxel (see Section 3.1), whose covariance matrix is ∑ ∑ ∑ x,y,z and the center of gravity is µ is the corresponding eigenvalue.The surface passes through µ and it's normal vector is − → n = − → v 1 = (n x , n y , n z ) T , so that the Equation ( 5) can be rewritten as Equation (6).
Using spherical coordinates, Equation ( 6) can be rewritten as Equation (7), where The detail of computing ∑ ∑ ∑ φ,θ and voting thresholds of coplanar voxels is described in Algorithm 3.With the use of first-order uncertainty propagation analysis [45], ∑ ∑ ∑ φ,θ could be derived from ∑ ∑ ∑ x,y,z as where J is the Jacobian matrix: where w = p 2 x + p 2 y .

Estimating Major Orientations
For a coplanar voxel K, its normal vector − → n k can be converted into polar coordinates (φ k , θ k ) according to Equation (7).Voting for − → n k is to vote for accumulator cell of (φ k , θ k ) and its surrounding cells.The voting value can be calculated from a multi-dimensional Gaussian Kernel (see Equations (10) and ( 11)).Since the covariance matrix is symmetric and positive definite, the two-dimensional Gaussian distribution probability density equation [46] is shown as where µ denotes (φ k , θ k ), − → δ denotes the distance with respect to µ, and ∑ ∑ ∑ denotes the covariance matrix about (φ, θ).It is worth noting that, when voting for a voxel with a small variance, the voting area of this voxel is small and vice versa.
For a coplanar voxel K, voting is carried out on the accumulator's cells where the corresponding (φ k , θ k ) is located.Iteratively voting is performed around (φ k , θ k ) until the calculated value of Equation ( 11) is smaller than the voting threshold of K (see Section 3.2.2).In this paper, we set two standard deviation to calculate the voting thresholds (see line 8 in Algorithm 3), which can guarantee a correct rate of 95.4%.However, in our experimental tests, the rock surfaces could be extracted efficiently by using one standard deviation to calculate the voting thresholds, and the efficiency of major orientation estimation could be improved greatly.
Besides calculating the voting value by Equation (11), the spatial length of voxel and the number of points contained in voxel are also considered [22].A weighting function is defined as where w a + w d = 1, voxel size is the edge length of the voxel and mesh size is the edge length of the bounding box.|K| is the number of points contained in voxel and |P| is the total number of the point cloud.The voting value received by the accumulator cell is the product of Equations ( 11) and ( 12).Since we use the hemispherical surface accumulator, the number of accumulator cells is not very large and the Sliding-Window method is used for peak detection.A window (N × N) is defined to slide around the accumulator.In this window, only the maximum value is retained.After sliding the window around the hemispherical surface, the window has passed through each accumulator cell and the major orientations has been obtained.

Rock Surface Extraction
After major orientation estimation, a set of major orientations is obtained.Voxel-based region growing is then performed according to these major orientations.The details of region growing based on voxels are shown in Algorithm 4, where { − → N } refers to the major orientations and {Voxel cop } refers to the coplanar voxels (see Section 3.1)."t 1 " represents the maximum angle between normal vector of coplanar voxel and major orientation.The value of "t 1 " should be as small as possible."t 2 " is the maximum angle threshold for region growing.

Algorithm 4: Region growing based on voxels.
Input : { − → N }: major orientations; {Voxel cop }: the coplanar voxels; t 1 : maximum angle threshold for searching for seed voxels; t 2 : maximum angle threshold for region growing; Max dis : maximum distance between a point to a surface or a surface to another surface; Output : {Sur f aces}: result of rock surface extraction; As shown in Algorithm 4, for a normal vector − → n , we select the seed voxel in the set of coplanar voxels, carry on the voxel-based region growing and update the set of coplanar voxels in turn.These operations are done iteratively until no seed voxels are found within the t 1 threshold.In this case, a surface cluster with similar normal vector − → n and its giant component can be obtained simultaneously.
Then the above operations for next major orientation are repeated.In Algorithm 4, the symbol "i" of Sur f ace ij refers to the ID of surface cluster and the symbol "j" refers to the ID of parallel sub-surfaces.
When looking for the neighbor voxels of q, the information of voxelization in Section 3.1 is used.As shown in Figure 3a, voxel q is a big-voxel, voxel A and voxel B are q's surrounding sub-voxels.They are coplanar but not adjacent.Therefore, voxel A and voxel B are not the neighbor voxels of q.In Figure 3b, A, B, C and D are not the neighbor voxels of q as well.The real neighbor voxels of voxel q need to be discussed in a variety of cases.
(a) (b) Figure 3.The neighbor voxels of q, (a) q is a big-voxel, A and B are sub-voxels.The relationship between A (or B) and q is coplanar but not adjacent; (b) None of {A, B, C, D} is the neighbor voxel of q.
Algorithm 5 shows the pseudo code of whether the neighbor belongs to Sur f ace ij .In Algorithm 4, {Neighbor} could be coplanar voxels, sparse voxels or non-coplanar voxels.For sparse voxels and non-coplanar voxels, each point in the voxel is judged by coplanar constraints (i.e., point-based region growing, see lines 2-8 in Algorithm 5).
It is worth noting that, due to the rough and uneven surfaces of rock mass, the coplanar voxels are often not "really" coplanar (as shown in Figure 4a,b).Figure 4c shows the finally result without considering non-coplanar points in coplanar voxels and it is inaccurate.The parameters of MinScale and MaxScale are introduced to overcome this issue, where MinScale and MaxScale are the minimum and maximum percentage of t 2 , respectively.For coplanar voxels, the voxel-based region growing is performed when the angle between the normal vectors is less than t 2 × MinScale, and the point-based region growing is performed when the angle between the normal vectors is greater than t 2 × MinScale and less than t 2 × MaxScale (see lines 9-21 in Algorithm 5). Figure 4d illustrates the result of considering the non-coplanar points in coplanar voxels.information of all point clouds is shown in Table 1.The display of all point clouds is presented in Figure 5. Notice that, the joint faces of Dataset Rock4 are not clear as other rock-mass point clouds.

Evaluation Metrics
To evaluate the performance of the proposed method, the evaluation metrics of precision, recall, F1 and execution times [31,[48][49][50][51] are used to quantitatively evaluate the results of rock surface extraction at the point level.Precision represents the percentage of correctly detected elements and recall indicates the percentage of reference ground truth data that are correctly detected.The F1 score balances precision and recall.If a is an extracted segment and m is a reference segment, the true positive (TP) = a ∩ m, the false positive (FP) = a − m and the false negative (FN) = m − a.The evaluation metrics of precision, recall and F1 are computed as Equations ( 13)- (15).

Parameter Turing
The excellent performance of the proposed method depends on the appropriate configuration of all parameters.Table 2 lists the main parameters used in the proposed method.As shown in Table 2, these parameters are divided into three groups.In clustering stage, the parameters of knn, voxel size , S α , ε, min bigvoxel and min smallvoxel are used to segment the original point cloud into coplanar, non-coplanar and sparse voxels.In addition, the value of knn, min bigvoxel and min smallvoxel depends on the densities of rock-mass point cloud.
Due to the point cloud of rock mass is rough and uneven, the constraints of coplanar clustering are relatively strict or loose depending on the actual situation of the point cloud.In this paper, several guidelines for coplanar constraints are as follows: • Constraint 1: each "target area" contains at least one coplanar voxel.

•
Constraint 2: coplanar constraints should be as strict as possible.

•
Constraint 3: the resolution of voxel should be as large as possible.
Because only the sample points within the coplanar voxels could participate in the major orientation estimation (see Section 3.2.3)and the selection of seed voxels (see Section 3.3), the Constraint 1 is established.If no coplanar voxel is detected in coplanar region of the point cloud (i.e., the "target" area in Figure 6a), the final result of rock surface extraction will certainly not be able to extract this "target" area (as shown in Figure 6b).In our tests, the parameters of S α , ε and voxel size can all affect the finally results in this case, but here, we only adjust the parameters of S α and ε to explain the efficacy of Constraint 1 (as shown in Figure 6c,d).Constraint 2 is to ensure the coplanarity of coplanar voxels.Constraint 3, which related to the parameter of voxel size , is designed to reduce the number of voxels and improve efficiency.Due to the researchers have been focusing on Hough Transform for decades [24,52,53], we just use the default value of parameters in the stage of major orientations estimation.
In the process of region growing, the parameter of Max dis is the distance between a sample point to a surface or a surface to another surface, which related to the roughness of the point cloud.The parameter of Min pla denotes the minimum number of sample points contained in a spatial surface.It related to the sample densities of the point cloud and the spatial extent of feature surfaces at the same time.The parameter of t 1 is the maximum angle between major orientation and the normal vector of voxel.It related to the roughness of the surfaces.
Table 3 shows the configuration of main parameters used in our experiments.Since the surface of the icosahedron point cloud have a large spatial extension, the Min pla value for icosahedron point cloud is greater than others in Table 3.On the other hand, due to the surfaces of the icosahedron point cloud are much flatter than other 4 rock-mass point clouds, its t 1 value is much smaller than others (as shown in Table 3).

Results for Synthetic Icosahedron Point Cloud
The icosahedron point cloud used here is a synthetic 3D point cloud and would not be affected by real LiDAR problems (i.e., registration, uneven density, etc.).It composes of 20 identical equilateral triangular faces (see Figure 7a).The parameter settings for icosahedron point cloud are shown in Table 3 and the result of our method is shown in Figure 7h.
We compare our method with 4 conventional surface extraction methods and 2 specialized rock surface extraction methods for rock-mass point clouds to validate our method's performance.RHT is the random Hough Transform with a new accumulator design proposed by Borrmann et al. [24].RG is the typical region growing method of Poppinga et al. [30].The RHT method and RG method are implemented in C++ code with PCL library.RAN_PCL is the typical RANSAC method [27] and RAN_CGAL is the optimized RANSAC method of Schnabel et al. [54].The implementation of RAN_PCL and RAN_CGAL is provided by PCL library and CGAL library respectively.HT_RG is a method of combining Hough Transform and Region Growing to detect rock surfaces in rock-mass point clouds [42] and the implementation code is provided by the authors.DSE is an open-source MATLAB software for semi-automatic rock mass joints recognition published by Riquelme et al. [6].
The time consumption of seven methods for icosahedron point cloud is shown in Table 4 and the time unit is second.Obviously, our method is the most efficient.In Figure 7, The results of (a)-(h) illustrate that each method can extract the surfaces correctly.Since the proposed method has been applied using a synthetic icosahedron point cloud, the next step will be to verify with real rock-mass point clouds.Here we used a total of 4 sets of rock-mass point clouds for testing and comparing the results with other 6 methods.

Qualitative Comparison
The rock surface extraction results of Rock1 by different methods are shown in Figure 8.Notice that, we use red polygons to point out the deficiencies of each result.The results of RHT and RAN_CGAL have the phenomenons of over-segmentation and discontinuity (as shown in Figure 8b,e).Figure 8c shows that the local surfaces of RG are incorrect and over-segmented.The result of RAN_PCL has the problems of over-segmentation, discontinuity and bad boundary (as shown in Figure 8d).The HT_RG and DSE method have the phenomenons of over-segmentation and missing some rock surfaces (see Figure 8f,g).Visually speaking, Figure 8h shows that the result of our method has better continuity and planar boundary without over-segmentation.For statistical analysis, please refer to Section 4.5.2.
We also test all methods on 3 other Rock-Mass point clouds.As shown in Figure 9, the first row represents 3 different Rock-Mass point clouds, and other rows represent the results of rock surface extraction by each method on 3 datasets, respectively.The results are visual similar, but our method is significantly faster.

Quantitative Comparison
To quantitatively evaluate the results of rock surface extraction, we manually split the Rock1 point cloud into a reference point cloud containing 20 surfaces by using the CloudCompare software (see Figure 10).In Figure 10b, the black part represents sample points without coplanar feature and each of the other colors represents a separated surface.Since the rock mass itself is rough and bumpy, there may be cases that points with coplanar characteristics outside the reference surfaces are not classified as rock surfaces.We compare and analyze the time consumption and the correctness of rock surface extraction in the following part.
The time consumption of different methods for Rock1 is shown in Table 5.It can be clearly seen from the table that the efficiency of our method is much higher than other methods, especially for the method of HT-RT and DSE.Table 5 also shows the number of detected surfaces of our method is fewer than other methods.The main reason is that the results of other methods have the phenomenon of over-segmentation (see Figure 8).The reference point cloud was labeled with 20 surfaces, and our method finally detected 23 surfaces.First of all, this is because our method has detected all the surfaces being labeled (see Figure 8h).Secondly, since the surfaces of the rock-mass point cloud are rough and uneven, we only mark up the surfaces with obvious coplanar characteristics (see Figure 10).In this paper, the time distribution of our method for icosahedron point cloud and rock-mass point clouds is shown in Table 6.We can find that the major orientation estimation and region growing are more time-consuming.The efficiency of major orientation estimation can be greatly increased by using one standard deviation to calculate the voting thresholds (see Section 3.2.3).The main reason of time-consuming in region growing stage is that points in coplanar voxels may not fully exhibit coplanar features (see Sections 3.1 and 3.3).While all the results of rock surface extraction are visually similar, We make use of three evaluation metrics to quantify the accuracy of the surface extraction results.As shown in Table 7, the precision, recall and F1 scores of our method are all the highest.Notice that, the recall scores of the other 6 methods are much lower than ours.This means that the other methods excluded many points near the boundary of feature surfaces (see Figure 8).Figures 11 and 12 show that the F1 score and recall score for each surface of our method is the highest.In our experimental tests, the HT-RG method performs better than other conventional feature surface extraction methods of urban point clouds in terms of accuracy.But as a whole, the proposed method performs better than the HT-RG method in terms of efficiency and accuracy.

Conclusions
This paper focuses on rock surface extraction in unorganized point clouds of rock mass.An efficient and accurate rock surface extraction method was proposed to deal with the challenge of roughness and unevenness in rock-mass point clouds.We segment point clouds via a framework with three stages, i.e., (i) clustering based on voxels; (ii) estimating major orientations based on Gaussian Kernel and (iii) rock surface extraction.To avoid the problem of computational complexity, the two-scale spatial grid and a bivariate Gaussian Kernel are used to estimate major orientations for selecting seed voxels.The voxel-based region growing is performed to ensure the connectivity of rock surfaces.Eventually, sets of parallel surface clusters and their giant component are obtained with high efficiency.
The proposed method was benchmarked with one synthetic icosahedron point cloud and four rock-mass point clouds.Overall, a high level of detail with better boundary in the rock surface extraction was achieved.When compared to ground-truth data of Rock1 by using the evaluation metrics of precision, recall and F1, the each score of our method was more than 91%.Additional, the proposed method is several orders of magnitude faster than other techniques for rock surface extraction in rock-mass point cloud.The experiments showed the availability and efficiency of the proposed method.
However, there is a limitation exists in our method, i.e., clustering based on voxels adapts fixed resolution for each voxel.As a result, the voxels may not retain well the boundaries of the point cloud.In future work, we will try to solve this problem with adaptive resolutions (such as supervoxel [48,55,56]).

Figure 4 .
Figure 4.The coplanar voxels may be not "really coplanar".(a) The result of clustering; (b) The partial enlarged view of (a); (c) The finally result of surface extraction without considering the non-coplanar points in coplanar voxels; (d) The finally result of surface extraction with Algorithm 5.

Figure 6 .
Figure 6.The results of clustering and rock surface extraction with different parameters.(a,b) Results of clustering and rock surface extraction respectively (S α = 25 and ε = 0.05).(c,d) Results of clustering and rock surface extraction respectively (S α = 22.5 and ε = 0.05).

Figure 8 .Figure 9 .Figure 10 .
Figure 8.The rock surface extraction results of seven different methods for Rock1.(a) The reference point cloud of Rock1.(b-h) Results of RHT, RG, RAN_PCL, RAN_CGAL, HT_RG, DSE and ours, respectively.Each color represents a surface.Notice that, we use red polygons to point out the deficiencies of each result.Ellipse indicates over-segmentation, rectangle represents discontinuity, parallelogram represents incorrect result, rhombus represent bad boundary and trapezoid represents missing detection.

Figure 12 .
Figure 12.Recall scores of Rock1's feature surfaces for different methods.

Table 1 .
The information of 5 point clouds.

Table 2 .
Summary of the parameters.Minimum number of points in a small voxel Related to the sample densities of the point cloud and voxel size .Set manually.

Table 3 .
The parameter setting of our method.

Table 4 .
Execution time of 7 different methods for icosahedron (in seconds).

Table 5 .
Execution time of 7 different methods for Rock1 (in seconds).

Table 6 .
Time distribution of our method (in seconds).