Effective Planar Cluster Detection in Point Clouds Using Histogram-Driven Kd-Like Partition and Shifted Mahalanobis Distance Based Regression

: Point cloud segmentation for planar surface detection is a valid problem of automatic laser scans analysis. It is widely exploited for many industrial remote sensing tasks, such as LIDAR city scanning, creating inventories of buildings, or object reconstruction. Many current methods rely on robustly calculated covariance and centroid for plane model estimation or global energy optimization. This is coupled with point cloud division strategies, based on uniform or regular space subdivision. These approaches result in many redundant divisions, plane maladjustments caused by outliers, and excessive number of processing iterations. In this paper, a new robust method of point clouds segmentation, based on histogram-driven hierarchical space division, inspired by kd-tree is presented. The proposed partition method produces results with a smaller oversegmentation rate. Moreover, state-of-the-art partitions often lead to nodes of low cardinality, which results in the rejection of many points. In the proposed method, the point rejection rate was reduced. Point cloud subdivision is followed by resilient plane estimation, using Mahalanobis distance with respect to seven cardinal points. These points were established based on eigenvectors of the covariance matrix of the considered point cluster. The proposed method shows high robustness and yields good quality metrics, much faster than a FAST-MCD approach. The overall results indicate improvements in terms of plane precision, plane recall, under-, and the over- segmentation rate with respect to the reference benchmark methods. Plane precision for the S3DIS dataset increased on average by 2.6pp and plane recall- by 3pp. Both over- and under- segmentation rates fell by 3.2pp and 4.3pp.


Introduction
Plane detection is a key aspect of a comprehensive point cloud analysis and segmentation. It is particularly important for indoor scenes, most of which can be easily decomposed into planar primitive shapes [1]. Point cloud processing (or its other form-depth map) is widely applied in such areas as reconstruction of human-made items, making inventories of architectural interiors, the documentation of cultural heritage structures [2], roofs detection in outdoor scans [3], or even for driver drowsiness estimation [4] and autonomous robot control [5]. It is also successfully used in civil engineering for deformation analysis [6]. Plane detection might also be applied to LIDAR-based rooftop reconstruction [7,8] and for primitive compression purposes [9], by storing shapes in a form of mathematical formulas rather than dense sets of points. The segmented points represent, among others, planes sampled during the scanning process. The problem of accurate and efficient plane detection, or more precisely, of extracting planar point clusters from an unorganized point cloud, has normally selected as a median along the considered dimension. It was used by Liu et al. in [27], who combined a kd-tree partition with a self-organizing fuzzy k-means clustering. They applied kd-tree to search for the nearest neighbours in order to model point set outliers based on the average point-to-point distance and the standard deviation. The authors made use of the clustering method and Ransom Sampling Consensus (RANSAC) to find appropriate planar patches. In fact, in [27] kd-tree was used for its dedicated purposes-nearest neighbours extraction. However, the authors present neither qualitative nor quantitative results of their studies. The validation dataset was not described either.
Concerning space partitioning, structure of a kd-tree itself might not be an optimal choice, unless division respecting points distribution is conducted. Although equinumerous divisions provides a well-balanced structure, selecting median value to construct division plane is not always a favourable choice since it may lead to many redundant subdivisions.

Plane Model Fitting
Hitherto mentioned point cloud subdivision methods produce points clusters, which should undergo planar regression. Taking into account the way the regression is processed, three main groups can be distinguished.
The first group involves transforming the data into model parameters space, which is referred to as a Hough transform [28]. In its pure form, in its generalization [29], or variants [30], Hough transform is a time demanding procedure, which does not deal well with noise [31]. Furthermore, the parameters space grid grows exponentially with a number of parameters.
Random Sampling Consensus (RANSAC) is a well known and widely used method for robust model fitting [10,11,32]. This approach is robust when considering the outliers, however, it does not handle noise (disturbances in points which may fit to a model). Furthermore, a fixed value of tolerance for inliers may be difficult or impossible to determine in advance; therefore, as Nurunnabi et al. observed, fixed threshold is one of the major problems of RANSAC based solutions [33].
Several methods relying on an effective estimation of the centroid and the covariance matrix for plane model fitting may be found in the literature. A solution of eigen-decomposition of a covariance matrix Σ, (Equation (1)) supplies three vectors, where the two longest indicate directions in which the plane spreads, whereas the shortest may be thought of as a plane's normal vector.
where n is a cardinality of a point cloud D = {p 1 , p 2 , p 3 , . . . , p n }(||D|| = n), µ is a centroid of the point set D, defined as in Equation (2).
To obtain a full plane equation, a point belonging to that plane has to be selected. In a perfect case, where all points are indeed sampled from a perfect plane, any point may be chosen for that purpose . However, because of the presence of noise in real cases, an average point is calculated as "the most representative" of a plane-centroid µ = [x, y, z], (Equation (2)). This strategy is usually referred to as the Principal Component Analysis (PCA). However, PCA works under strong assumptions: 1.
All points are sampled from the same ideal plane (which means the absence of any outlying points), and 2.
Noise is symmetric so that its impact on the solution is minimised.
Whereas the second assumption may be justified with the central limit theorem [34] (which says that a number of random distributions will tend to form a bell curve-which is characterized by a symmetric shape), the first one is difficult to defend since a considered set usually contains outlying points. This results in extra operations performed prior to the analysis.
In order to deal with the outlying elements, Hubert et al. introduced a robust version of PCA (ROBPCA) [35]. The authors determine a set of outliers by means of a modified Stahel-Donoho outlyingness [36,37] using a centroid and a covariance matrix rather than the median and the median absolute deviation.
As the authors pointed out, their method is extremely time demanding even for small datasets. In addition, the method proposed by Hubert et al. searches for global outliers, which as indicated by Nurunnabi et al. in [33] may cause unreliable and unexpected output in the case of point clouds.
To overcome this problem, Nurunnabi et al. proposed to reject outliers locally (per point) to make the vicinity of each point accurate and representative for the plane fitting purpose [33]. In order to achieve this goal, Nurunnabi et al. introduced an algorithm for determining the Maximum Consistency Set (MCS) [33]. It starts with the random selection of three points and fits a plane with PCA-based approach on that subset. Then, having calculated a point-plane distances for all neighbouring points, another subset of points with the smallest values of the Euclidean distances is selected and another plane is fitted. This procedure is carried out iteratively in a manner similar to the concentration step (C-step) in the FAST-MCD algorithm by Rousseeuw and Driessen [38], which minimizes the covariance determinant by iterative selection of a subset of points with the shortest Mahalanobis distance. Final plane parameters are those, which result in the smallest eigenvalue of the covariance matrix of the considered points subset.
Based on parameters of the plane computed for MCS, [33] proposed two methods to identify the outlying points. The first uses the measure of the Robust z-score (Rz-score) [33]. The second relies on the so-called, Robust Mahalanobis Distance (RMD), where rejection of the outliers is performed by screening out the points lying farther than 97.5-th quantile of the Chi-square distribution. Robustness of the methods is an effect of computing both a covariance matrix and a centroid based on the Maximum Consistency Set (a set with removed outliers) rather than the entire collection. Robust Mahalanobis Distance (RMD) is used solely for a slight reduction of the noise in the considered set of points. The authors reported accuracy exceeding both ROBPCA and RANSAC algorithms. In spite of quality improvements and reduction of the calculation time, the method presented in [33] leaves some areas for improvement. The first problem concerns determination of the vicinity of the points. Whereas the authors suggested using the kNN strategy, they did not indicate how to determine a value of k. This is a crucial issue, since the size of the vicinity influences determination of the MCS. Moreover, selection of the points whose vicinity contributes to more than one plane may produce an unexpected Maximum Consistency Set ( Figure 1). As may be seen therein, for point clusters representing plane cross-sections, both methods selected points from the expected set. However, in case of the MCS (similarly to RANSAC), three points were detected, which in the case of noise, badly approximate real plane (see estimated plane marked in dark grey in Figure 1). On the other hand, FAST-MCD considered more core points, better approximating estimated plane model ( Figure 1), but it takes much more time than MCS. An attempt to solve the issue of an adequate neighbourhood choice may be found in the studies of Eckart et al. in [13] or this of Xu et al. [39]. The authors of [13] dealt with point cloud registration rather than plane detection, they presented an interesting way of space partition which favours planar fragments. Eckart et al. employed a top-down recursive point cloud partition, making use of Gaussian Mixture Models (GMM) with an early stopping heuristic, which terminates procedure when the curvature of the points in a node C λ , (Equation (3)) is beneath a threshold. This procedure adaptively sets proper division scale to produce efficient partition reducing oversegmentation.
where λ 1 , λ 2 , and λ 3 are eigenvalues of the covariance matrix of the considered points. On the other hand, Xu et al. applied octree-based voxelization coupled with region growing. Such voxelization imposes a trade-off between details preservation and processing performance [39].
Yet the other strategy is represented by an algorithm of Dong et al. [22]. Their method processes partition and a hybrid region growing in order to reduce a number of candidate planes. The final stage is a global energy optimization procedure, which in fact is constructed by two nested optimization strategies: simulated annealing and an extension of the α-expansion algorithm introduced in [ [35] or MCS [33], Dong et al. focused on high granularity of space partition rather than on rejection of the statistical outliers in order to construct robust planar fragments. This, in turn, significantly increases the processing time increase, accompanied by an excessive oversegmentation. Nonetheless, [22] reported the best results validated on the popular and recognized dataset (S3DIS), which makes it a reference method for precision comparison.
In contrast to the current leading approaches, the motivation for the presented strategy is to exploit MD for a selection of the core points rather than rejection of the outliers. Eigenvectors-driven shifts of a default centroid, supported with a selection of the MD core points, should result in a better and faster fitting of the plane model. According to our hypothesis, core points selected semi-deterministically, with the shortest MD, are more efficient than randomly selected model constructed by reducing the noise with the MD or FAST-MCD algorithm.
From the above review of the state-of-the-art, just a few authors [10,22] considered space partitioning simultaneously respecting point set characteristics and the fitting of the plane model. These authors recommended other approaches than presented in this paper. In most cases, just one plane detection stages were thoroughly studied. However, none of the authors explicitly addressed the problem of excessive point cloud oversegmentation. Methods conducting space subdivision find nodes comprising corners and edges difficult to split efficiently. They are treated as many others non-planar, cluttered nodes and subdivided into smaller patches.
Our contribution relies on a better fitting of a plane model within non-planar point clusters, by means of a MD-based, selection of highly planar core points around six, adaptively shifted centroids. Preceding it, kd-tree like histogram-driven space partitioning method (hd-kd-tree), assures planarity-oriented point cloud partition and results in lower over-segmentation, which, in turn, increases time efficiency.

Proposed Method
The method of a point cloud planar clusters detection can be divided into four stages ( Figure 2): initial point cloud alignment (preprocessing); 2.
planar patches robust refinement; 4. refined planar patches aggregation and post division;

Initial Point Cloud Alignment (Preprocessing)
Let a set D of n points D = {p 1 , p 2 , p 3 , . . . , p n } with cardinality (||D|| = n) be defined. Each point p i ∈ D is defined as a tuple of three coordinates Often, a point cloud D may be misaligned with respect to an object inherent coordinating system (Figure 3a) due to several reasons. For instance, imperfect laser device placement during scene capture. If the dataset is oriented better with respect to the axes, it is very likely that at least some of the planes might be extracted without an unnecessary subdivision. Certainly, the benefits of this stage depend strongly on the point cloud, but in the case of interiors, it may be assumed that many elements are mutually orthogonal. Markiewicz used computer vision algorithms, based on features matching for successful point cloud orientation [42]; however, this method needs a reference model with respect to which orientation will be determined.
In the problem considered in this paper, such a reference model is not known in advance. Hence, it was decided to align the point cloud by minimizing the product of its dimensions.
Let a sample point cloud represents the interior of a room. Very often, the set is composed of mutually parallel or orthogonal planes, like walls, floor, ceiling, tables, and others. This knowledge might be used to optimize a way the space is partitioned. In fact, the dataset is supposed to be transformed so that the major planes are aligned with the coordinate system axes. It can be seen from Figure 3a that an indoor misaligned set usually has a larger Axis Aligned Bounding Box (AABB) than the properly oriented point cloud (Figure 3b) [43]. In view of the above, a proper alignment procedure can be formulated as an optimization task, where the objective function takes into account the AABB dimensions (Equation (4)).
whereR is the resulting rotation matrix; W x , W y , and W z are dimensions of the AABB (x max − x min , y max − y min , z max − z min ) along the x, y and z axes respectively. Because it is a trivial optimization task, one of the simplest, fast, gradient-free Nelder Mead method [44] was applied. As a result of this procedure, the optimal rotation matrixR is calculated. This matrix transforms a point cloud so that the volume of the AABB is minimized and a point cloud is better aligned with the coordinate system axes. This should enable extracting larger planar fragments, especially for indoor scenes.

Histogram-Driven Point Cloud Partition
The aim of this stage is to subdivide the considered point cloud D into major, inherently planar groups O j : j ∈ {1, 2, 3, . . . , m} for which change of curvature (Equation (3)) is low and one additional group (O out ) representing remaining points.
Kd-tree structure seems to be better-suited for extraction of the planar fragments than the octree. Therefore, it was selected for space partitioning. A basic kd-tree recurrently produces equinumerous sets of points resulting in a well-balanced tree. However, in the ordinary kd-tree building procedure, the division plane is supposed to be constructed based on a median value along the successive dimensions. Having obtained better oriented point cloud (in the stage 1 from Section 3.1), subdivision separation values can be selected making use of histogram peaks of point coordinates in order to reduce the number of redundant subdivisions and to force favourable planar points concentration.
To do so, three histograms H x , H y , H z (H x/y/z for short), describing the distribution of the point coordinates along three axes x, y, z, are defined. It must be noticed that the distribution of points might be different for the three coordinating system directions.
Assuming that numbers of bins for x, y and z axes take values b x , b y , and b z respectively, three sets of corresponding bin border values, B x , B y , B z (B x/y/z for short), are calculated as in Equation (5).
where d is an empirical coefficient of inlier plane tolerance threshold (Figure 4), discussed later in this section, determining the range of points constituting individual histogram bin.
Histograms H x/y/z may be formally expressed as shown in Equation (6).
where || · || represents the cardinality of a point set constituting the selected bin.
A crucial value of each histogram is a number of bins. It should be sufficient to differentiate close parallel planes, yet large enough to encompass most of the planes described by the histogram bins. If d is selected as an inlier-plane tolerance (Figure 4), then the number of bins b x , b y , b z (b x/y/z for short), considered separately for each dimension x, y, z, might be defined using Equation (7). The symbol · represents the floor function (integer part). An example of a histogram along the z axis for the selected (Area1/office-19) point cloud is presented in Figure 5. Histograms should be examined for bin height values as they reveal the concentration of points throughout the whole set D. For instance, considering an interior scene, if the first and the last bin of one of the axes histograms contain the most points, this suggests that the large planar fragments are contained therein. Note that some high values of histograms' bins along other dimensions might also imply the presence of the bigger planes. Therefore, they might be considered as foundations for constructing point cloud division planes. It does not matter if the point cloud is rotated by 90 • and a ceiling will be vertical whereas the walls will be horizontal because the semantic analysis is not applied for plane. Furthermore, even in the case of a poorly defined point cloud, the larger part of the planar fragments can be detected. Histogram construction is of linear complexity (O(n)).
Peaks of the histograms could indicate regions where larger planes are contained even though real planes might not necessarily be orthogonal to the coordinate system axes. With this in mind, the division planes are selected based on the inner border values (B x/y/z ). The highest bin border value, closer to a value representing half of a point set dimension (min D x/y/z + W x/y/z 2 ), along the corresponding axis, is selected. Assuming a plane is defined as in Equation (8) and H x/y/z [bin] is the highest histogram for the selected z/y/zcoordinate, a division plane is calculated based on the associated value B x/y/z [bin] as presented in Equation (9). π : n · p = −ρ (8) where n is a unit-length normal vector, p is an arbitrary point belonging to the plane, and ρ is the plane intercept, defined as a distance between that plane and the origin of the coordinate system.
Division of a point set with a plane (Equation (9)) produces two child nodes containing points being on the left (or exactly on) or right-hand side of that division plane. Planarity check is carried out by the analysis of the covariance matrix of points contained in each node. A good estimator of the planarity is a low curvature change, expressed with Equation (3). This determinant has been widely used by researchers [22,33,45,46].
Each subsequent plane subdivision is calculated based on the highest (or the first highest value if there are more), yet not considered bin, among all three histograms H x , H y , H z . Subdivision is continued until the stop condition is reached. Even though many planar point clusters in a point set D may not be orthogonal to the coordinate system axes, applying a histogram for points space partitioning reflects their distribution and may improve the result because of the extraction of larger planar fragments instead of their deterministic division.
There are two stopping conditions for the recurrent subdivision. The first stops recursive calls if a number of points, in a current node, is lower than the minimum (n = 20 [22]-user-defined minimum cardinality of a valid node). The second is the planarity condition C λ < C (Equation (3)), where C is an assumed threshold, which says how many points are spread along the candidate normal vector ( C = 0.002 in the presented method) When the points planarity test in a node was passed, the set of points was appended to the set of patches O j satisfying the planarity condition (C λ < C ). If not, the division plane is selected as indicated earlier in this section-based on the currently highest histogram peak. Original histograms are exploited throughout the whole partitioning process, in spite of the fact that the sub-histograms of the selected nodes (subsets of points histograms) may suggest different subdivision planes.
It may happen that all border values B x , B y , B z were already used to construct division planes but some points remained because of failure of the planarity test. In this situation, the space subdivision follows the ordinary kd-tree partition procedure. It involved selecting as a division plane one based on the median of coordinate values along the current axis.
For clarity, the partition procedure is presented in Algorithm 1.  Assuming the histogram is calculated once prior to the subdivision procedure, the procedure complexity is of the order O(n · log n).

Algorithm 1 Histogram-driven partition algorithm
Obviously, the superiority of the proposed method is strictly tied with the presence of the orthogonal and parallel planes. If planes of other orientations exist, they will not be segmented entirely in a single iteration, but they are supposed to be divided into parts according to the assumed threshold of the histogram bins.

Planar Patches Refinement
The previous stage of the method gives a set of patches, retrieved from the histogram-driven kd-tree (hd-kd-tree) nodes O j : j ∈ {1, 2, 3, . . . , m} satisfying the planarity condition, and the group of remaining points O out . The patches are roughly planar because the measure of curvature (Equation (3)) is a relative value rather than an absolute one. Hence, it may happen that the planarity condition is satisfied but some outlying points can still be contained therein (Figure 1). If so, these planar patches need additional processing to fit a plane properly.
Mahalanobis distance (MD) of a point set, whose sample covariance matrix is Σ and sample meanµ, is defined as in Equation (10).
Usually, a standard MD-based rejection of outliers is performed by screening out points lying farther than the 97.5-th quantile of the Chi-square distribution. Since, MD indeed follows Chi-square distribution (χ 2 ) with the number of degrees of freedom set to 3 (because there are three variables: x, y, and z), outliers are the points being farther than χ 2 3,0.975 = 3.057 from the centroid µ [33,35,47]. This works under the assumption of a Gaussian distribution of points, but this assumption may not be valid globally. Nevertheless, at a sufficiently small scale, it might be a good local estimate of a surface, as indicated by Stoyanov et al. in [26].
Even though a covariance matrix and a centroid are robustly estimated with FAST-MCD [35] or MCS [33] procedures, nodes containing points close to corners or edges may yield improper estimation where the crucial issue is to find points which actually belong to the main plane. Even though the randomness of MCS and FAST-MCD forces a result of the plane fitting to tend to a global optimum but to ensure it, many iterations may be required. Hence, in the proposed method, the stochastic space exploration, which is applied in the MCS and FAST-MCD procedures, was replaced by performing C-step with respect to, so called, cardinal points ψ i . Eigenvectors become predominant directions of alternative centroid candidate (cardinal points) locations as they represent inherent distribution of the point variations. The cardinal points are calculated by shifting the default centroid of a points cluster along its positive and negative eigenvectors by distances resulting from an analysis of the variance. This leads to a fast and accurate result of plane fitting.
The first cardinal point (ψ 1 ) is a default centroid µ, which is a centre of the mass of a point set. However, the C-step calculated with respect to µ may converge to core points distributed across several planes in case of their unfavourable points distribution. Hence, C-step should be carried out with respect to relocated centroids (cardinal points) to verify convergence around other possible centroids.
The most promising directions, which indicate possible centroid shifts, are those indicated by the eigenvectors of the covariance matrix Σ. Even though Σ may be disturbed by outliers, eigenvectors would point either towards these outliers or towards candidate centroids of a set of core inliers (see Figure 6). As a result, one of the shifted centroids (cardinal points) might be better than the default one. Hence, the shifts are defined as appropriately scaled eigenvectors e 1 , e 2 and e 3 of the covariance matrix Σ and their opposites (−e 1 , −e 2 and −e 3 ). This results in seven cardinal points ψ 1,...,7 (in a 3D case there are six points shifted along positive and negative eigenvectors plus the default centroid). Next, the issue of a proper shift range has to be solved. It should be radical change of a default centroid position in order to force C-step to verify also boundary regions of the sub-clusters. Let a portion h = 50%, which provides the highest breakdown point [48] be used for C-step. Then, knowing that eigenvalues are variances along the eigenvectors [13], the other half of a set-according to a univariate Gauss distribution-is contained outside the ranges along the corresponding eigenvectors e 1 , e 2 , and e 3 . Hence, seven cardinal points may be defined as shown in Table 1. Table 1. Seven cardinal points ψ i . µ is the centroid of the contaminated dataset; e 1 , e 2 and e 3 are dataset eigenvectors corresponding to λ 1 , λ 2 and λ 3 eigenvalues.

Point Aggregation
The successive stage of the proposed method aims at assigning points p i across the entire point cloud D to appropriate planes (π j ) determined in Section 3.3. The stage begins with calculating core points for the local vicinity of each point p i making use of Algorithm 2. Based on those core points, the normal vector is estimated by fitting a plane to them. Vicinity of a point p i , from which core points are selected, was set to the 7 nearest neighbours according to the method presented in [49].
Next, the global procedure of point-plane assignment is performed for the entire set D. This is done by checking an angular deviation between the refined planes normal vectors n π j (Section 3.3) and each points' normal vectors (n p i ) estimated based on their vicinity. Points whose normal vector deviates not more than assumed θ are assigned to the given plane (Equation (11)). arccos(n π j · n p i ) ≤ θ (11) where n π j is the robustly estimated normal vector of a j-th plane (π j ) and n p i is the robustly estimated normal vector of a point p i . while applying angular deviation threshold, many parallel planes could be joined together. To overcome this problem, density-based clustering was employed (for instance, k-means based density clustering [50] or HDBSCAN/DBSCAN [51]). HDBSCAN was used in the presented approach because of its superior efficiency, reported by [51]. As an output of the method, a set of detected planar points clusters, assigned globally to robustly selected planes π j , is obtained.

Datasets
The notion of a benchmark dataset for planes detection task is not well established in the literature. Almost every reported method used a different dataset, including artificially generated ones. In [22] and [52], S3DIS [41] dataset was used. However, it contains points labelled with respect to objects' adherence (i.e., a chair, a table, a lamp) rather than the individual planes. On the other hand, Li et al. [10] made use of a laser scan (Room-1) from Rooms UZH Irchel dataset [53], in spite of the fact that it does not contain any labelled data. The S3DIS dataset shows much sparser point clouds density than in the case of rooms from UZH Irchel dataset. These datasets differ significantly in terms of accuracy, noise, scan shadows, cardinality, and scene complexity. Therefore, it was decided to use representatives from both of them to verify the proposed method on point clouds of varying nature.
The present study uses point clouds of the S3DIS dataset and the Room-1 point set [53] used by Li et al. [10] (Table 2).
Because for the Room-1 dataset no ground truth segmentation was provided, it was labelled manually. The ground truth segmentation of the S3DIS dataset, in turn, was manually modified to represent individual planes. An example of six point clouds from S3DIS and Room-1 datasets are presented in Table 2.
For the space partition juxtaposition, four values were presented: the division tree spreadness (number of all nodes in division tree), the final number of groups, the number of points, which remains after the partition process was accomplished, and the partition time. All space partition methods were tested on the S3DIS dataset with the same setup (Table 3).  (12)), plane recall (Equation (13)) and over-as well as under-segmentation rates (Equations (14) and eqrefeq:usr) were used as the validity measures of the entire procedure. 13) where N C stands for a number of correctly segmented planar clusters (with maximum overlapping strategy, 80%), N S represents the total number of planar clusters obtained as the algorithm output, and N G is a number of ground truth planar clusters.
where N U is a number of resulting planar clusters that overlap the multiple ground truth planar clusters. N O , in turn, is a number of ground truth planar clusters overlapping more than one resulting planar cluster.

Space Partitioning Results
All methods used the same parameters values (see Table 3). For VCCS, remaining parameters (Min r, Max r, ∆r) were taken from the studies of Dong et al. who used the method for the S3DIS dataset [22]. The symbol Ω represents an average distance between a point and its closest neighbourhood.
The results of the experiments are presented in Figure 7. Figure 7a dhows that the proposed method (hd-kd-tree) yields, on average, similar tree spreadness as the octree (close to 7500 nodes); however, with a higher variance. Both kd-tree and, especially PCP, produce much more spread trees. In case of VCCS, tree structure is always very flat (up to 10) due to constraints put on the partition (minimum supervoxel size-Min r, and maximum supervoxel size-Max r). However, in spite of the flat structure of the VCCS tree, Figure 7b shows that this tree has a very wide structure and produces many more clusters ( 260,000) than the hd-kd-tree strategy ( 100,000). The octree partition results in even more clusters than VCCS. Figure 7c clearly shows that the proposed method preserves significantly more points (on average, more than 95%) than the other methods. In this comparison, also ordinary kd-tree is superior versus octree, PCP, or VCCS. Figure 7d, in turn, shows time required to perform partition. All methods, except VCCS, show similar time demand (2.5 s). VCCS usually requires substantially more time (17.0 s).

Planar Patches Extraction Results
This section presents the results of the proposed method in comparison with the outcome of the approach of Dong et al. [22] and Li et al. [10], which are the most competitive state-of-the-art procedures.
The approach of Dong et al. partially exploits the advantages of a hierarchical partition; however, due to constraints put on the cell sizes at each level, it still introduces many redundant subdivisions.
Time complexity of these partition algorithms is presented below. The way Dong et al. divided the space suggests at least O(n 2 ) time complexity because of the application of the region growing k-means algorithm, assuming supervoxels sizes to be constant and fixed in advance. The proposed partition algorithm has a complexity of the order O(n log n).
An outcome of the proposed method, for the test point clouds presented in Table 2, is shown in Figure 8. Note, that most of the planes were correctly detected both inside and outside the room.
Tables 4 and 5 present the results and comparison with those obtained by Dong et al. (Table 4) and by Li et al. (Table 5). The results of the SMD-based method results and the method of Li et al. [10] results are compared in Table 5. The results were obtained for Room-1 dataset [53].
In terms of PP and PR metrics, the experiments reveal that the SMD-based method outperforms (by 2.6pp, 3.0pp) the reliably documented method of Dong et al. Moreover, the newly developed points cloud histogram-driven partition is characterized by much lower over-and under-segmentation rates (by 3.2pp, 4.3pp). In comparison to the method of Li et al., the proposed method reaches higher plane precision (by 1.5pp) and slightly higher plane recall (by 0.2%). OSR and USR were not reported in [10].  The limitation of the proposed method concerns non-planar object detection that is decomposed into planar sets, like a trash bin (Figure 9), but such objects should be considered individually with other methods. The presented results clearly confirm the contribution stated at the beginning of this paper. Robust plane refinement, based on SMD coupled with hd-kd-tree points cloud partition and its simple aggregation, yields better results than the current state-of-the-art methods. SMD method tends to converge to actual planar fragments, like FAST-MCD and MCS (Figure 10), yet within a much shorter time (Table 6). Partition with hd-kd-tree subdivides the space in such a way that fewer points are rejected compared to other methods ( Figure 8) and the output is supplied in reasonable time virtually equal to the octree, which had the best time so far. (c) Figure 10. Comparison of (a) FAST-MCD, (b) MCS, and (c) SMD methods of core points determination (core points marked with green) for exemplary, roughly planar points cluster (blue points) containing a fragment of other planar group (orange points). Fitted plane is depicted in dark gray. Table 6. Time comparison for FAST-MCD, MCS and the proposed SMD method for core points selection for points presented in Figure 10.

Conclusions
The experiments described in this paper were focused on searching for an effective method of a point cloud subdivision process and efficient clustering of the set of planar points. As a result, a new efficient method exploiting histogram-driven kd-tree structure was introduced. It processes recurrent and orthogonal point cloud subdivision. A point cloud was divided efficiently into a set of planar point fragments. The experiments have clearly shown that the proposed partition method is superior over the benchmark partition approaches presented in the literature, in terms of the number of resulting clusters and portion of the preserved points. It extracts a lower number of larger planar clusters than other state-of-the-art methods.
The conducted experiments have revealed that semi-deterministic strategy relying on the node core points selection, based on Shifted Mahalonobis Distance (SMD), enables precise, reliable, robust, and fast estimation of the plane parameters.
The proposed method was verified for sparse and noised datasets like these from the S3DIS database. The method correctly detects planes with the average plane precision (PP) score of 93.0% and the average plane recall (PR) at the level of 94.4%. It exceeds the examined reference methods by 2.6pp and 3.0pp respectively. Both over-and under-segmentation rates were relatively low and better than in the state-of-the-art methods (for the proposed method they fell to 4.4% and 3.9% respectively). The over-segmentation occurs mainly for curved objects that could not be approximated well by a single plane (like a trash bin) and should be treated separately.
Selective experiments to determine efficiency of the method were also conducted for the dense point cloud-Room-1. Plane precision reached 86.4% whereas plane recall increased to 98.5%.
In the case of a laser scan, the over-segmentation rate is higher than in the case of sparse datasets because of disturbances in point distribution in areas where the laser beam falls under very small or very large angle (ceiling above the laser head or floor just below the laser head). Li et al. did not report OSR nor USR rates.
Conducted surveys have demonstrated that application of the SMD-based plane fitting procedure allows for the selection of clusters of core points that are more representative than the other methods and the resulting planes better fit the considered clusters. Furthermore, application of a histogram-driven kd-tree partition yields a more balanced partition than other state-of-the-art partition procedures. This research also reveals a noticeable potential for further generalization, especially for non-planar, semantic model-based detection in point clouds segmentation tasks.