An Improved Ball Pivot Algorithm-Based Ground Filtering Mechanism for LiDAR Data

Automatic ground filtering is an essential step for Digital Elevation Model (DEM) generation, which has significant application value. However, extraction and classification of ground points from the Light Detection and Ranging (LiDAR) data, especially in multitudinous terrain situations, is a challenging task because it is difficult to determine the set of optimal parameters for removing various non-ground features. In this paper, a new ground filtering technique based on an improved Ball Pivot Algorithm (BPA) is proposed. At the beginning, the LiDAR point cloud dataset was divided into different subsets based on the 2D regular grid. The lowest point in each grid was selected as the seed point to build a single-layer surface. After that, the improved BPA was executed to remove points on the higher location. Then, the rest of the points were calculated and selected as a new seed point according to the spatial relationship with the initial surface. Finally, non-ground points were filtered by means of improved BPA traversing all the grids. Our experimental results on the Benchmark dataset provided by the International Society for Photogrammetry and Remote Sensing (ISPRS) Working Group III/3 showed high accuracy (with a mean kappa coefficient over 80%) in terms of completeness, correctness, and quality for DEM generation. The experimental results demonstrated the proposed method is robust to various terrain situations, as it is more effective and feasible for ground filtering.


Introduction
Digital Elevation Model (DEM) generation is an important research problem in the remote sensing area, with good application prospects in the fields of smart city, 3D mapping, and virtual reality. Compared with traditional photogrammetric methods, Light Detection and Ranging (LiDAR) shows its advantages, especially in height change detection, so that it is gradually being adopted as the primary technique to produce DEM. The laser-scanned three-dimensional cloud of points not only includes bare ground, but also contains buildings, vehicles, trees, etc. DEM generation requires a filtering step, which separates the points that belong to the ground from the ones that do not. Therefore, the quality of the ground filtering algorithm will substantially impact the accuracy and visual effects of the DEM.
Nowadays, increasingly more researchers have focused on ground filtering issues for DEM generation [1][2][3]. In fact, the methods used to solve this problem vary from machine-learning-based filters and Triangulated Irregular Networks (TIN) ones to filters that rely on morphological properties. In addition, the integration of LiDAR data and other data sources for land cover classification has gained popularity [4]. However, most of the existing works either require a large number of samples for model training or need to Remote Sens. 2019, 11,1179 3 of 16 points from non-ground ones, by selecting the lowest voxels with a value of 1 as ground seeds and then labeling them and their 3D-connected set as ground voxels [32]. The surface-based filters were widely used in DEM generation. Nonetheless, the surface method tends to fail in preserving terrain details and small non-ground object classification. Overall, these new methods have proved their competency in building reliable DEMs, even under various challenging terrains. However, limitations still exist. The ground filtering methods are usually applied to large and diverse areas. Most of the approaches proposed require proper parameters or thresholds for deriving DEMs under different terrain situations. Their utility is rather diminished in DEM generation by the complexity of the parameter adjustment.
Our goal in this paper is to develop a ground filtering model based on 3D alpha shapes (α-shape) with fixed and robust parameters to apply to different environments and conditions, which will significantly reduce labor cost and increase practicability. The Ball Pivot Algorithm (BPA) and TIN construction are two main ways to implement the α-shape, which has been widely used to find the outline of an unorganized set of data points [33][34][35][36]. Recently, the TIN construction method has been popular in LiDAR data processing and ground filtering. Demir use the α-shape algorithm to vectorize the roof outlines based on LiDAR data [37]. Paavo produced an alternative TIN model for stone detection [38]. However, there are still limits. On the one hand, it is difficult to find a perfect value of alpha to exclude the points that reflect the detailed surface. For this problem, Wu proposed an adaptive α-shape algorithm [39]. However, the algorithm can only consider a 2D profile, which loses the 3D spatial information. On the other hand, even if the alpha shape radius is chosen suitably, it is still hard for the TIN construction method to divide the ground part from the 3D point clouds due to the hollow formed by TIN construction in bare ground regions, which severely impedes the application of the α-shape algorithm in 3D space.
In this paper, a novel filtering method is proposed to enhance the efficiency and accuracy of ground filtering under varied terrain conditions, which combines the BPA and spatial sorting methods. The BPA is further improved so that only the bottom points are considered in the calculation as the ground points and no 3D TIN construction is required. The proposed method not only effectively reduces the computation, but also avoids the loophole caused by the single-layer point cloud inability to construct the tetrahedron in 3D TIN. The applicability of the method is proved in the Benchmark dataset provided by the ISPRS Working Group III/3, and results show a moderate performance in terms of the average total error. Our contribution in this paper can be summarized as follows: (1) A novel automatic method for ground filtering is proposed, which does not rely on the model training or complex parameters adjustment. Once verified by standard data, the proposed method can keep working in a relatively stable manner, even when the environment changes, as it is more effective and feasible for ground filtering.
(2) An improved BPA is designed in the proposed filtering model, which directly extracts bottom boundary points without the process of constructing a 3D TIN model, and effectively solves the problem of output loss in bare ground regions.
The organization of this paper is as follows. Section 2 introduces the details of the main method and the limitations of BPA. The area of interest and data set is covered in Section 3. Section 4 reports the experimental setup and discusses the experimental results. In Section 5, the conclusions and future work are presented.

Overview
The main idea of the proposed method is identifying bottom points directly based on the improved BPA. Firstly, the raw point cloud was divided into different sub datasets by using 2D regular grids, and a pre-operation was designed to remove outlier points. Then, the boundary points were extracted by improved BPA. The details of the proposed filtering algorithm are illustrated in  Firstly, the raw point cloud dataset was divided into different sub datasets by 2D regular grids according to the horizontal position of each point. Then, the points in a same grid are sorted by elevation in ascending order. Once the outlier points are removed, the lowest point in each grid is collected to build a single-layer surface (initial surface) by using the strategy of 2D TIN. By improving the BPA in a relatively large scale, a point cloud at a higher location can be excluded from the surface. Meanwhile, the initial surface is adjusted according to the deletion of relative higher points, and the distance between each point and its nearest facet is calculated. Then, points with small distances from the initial surface are selected as the new seed point. The improved BPA is applied to search ground points from those starting points iteratively. Finally, non-ground points are filtered by traversing all the grids.

BPA and Its Limitation for Point Cloud Filtering
The BPA is one of important ways to implement alpha-shape, which treated as a linear approximation of the original shape was designed for shape reconstruction from an unorganized finite set of points [34]. The main idea of BPA is putting a circle on the point set and rolling it along the outside of the points in 2D space. If two points in dataset are tangential to the circle and no other points in the circle, the two points are considered as boundary points. In 3D space, the circle is replaced by a ball and the two points become three.
In terms of its definition, α-shape (∅) is the boundary shape consisting of selected points from the original point set (S). The resulting boundary shape is determined by a positive parameter (α). The formula can be described as follows.
where α is the reciprocal of the planimetric Euclidean distance between the current point and its next neighbor.
where dist S i , S j is a function to calculate the planimetric Euclidean distance between S i and S j [35]. The BPA starts by putting a ball in contact with three points (forming thus a seed triangle). The ball with a radius of r is rotated around two of these three points; it remains in contact with the two points and rotates around the axis they form until it touches another point to form another seed triangle [40]. The alpha-shape algorithm works properly in searching the surface anchor points for relative even point clouds, but it is difficult to find a perfect alpha value α (or a ball with a certain radius r) that does not exclude too many surface points that reflect the detailed surface. Besides, the alpha adaptive shape without qualified seed ground points will consume a lot of computation to adapt the size of the ball repeatedly [39]. Therefore, the alpha shape algorithms for ground point identification and segmentation in 3D space are limited by the computational complexity. Even if the alpha shape radius is chosen properly, it is still hard for the traditional algorithm to realize ground filtering from the 3D point clouds because of the hollow problem in the bare ground area, such as the hollow formed by TIN construction in the bare ground region (Figure 2b).
Remote Sens. 2019, 11, x FOR PEER REVIEW 5 of 16 two points and rotates around the axis they form until it touches another point to form another seed triangle [40]. The alpha-shape algorithm works properly in searching the surface anchor points for relative even point clouds, but it is difficult to find a perfect alpha value α (or a ball with a certain radius ) that does not exclude too many surface points that reflect the detailed surface. Besides, the alpha adaptive shape without qualified seed ground points will consume a lot of computation to adapt the size of the ball repeatedly [39]. Therefore, the alpha shape algorithms for ground point identification and segmentation in 3D space are limited by the computational complexity. Even if the alpha shape radius is chosen properly, it is still hard for the traditional algorithm to realize ground filtering from the 3D point clouds because of the hollow problem in the bare ground area, such as the hollow formed by TIN construction in the bare ground region (Figure 2(b)).
(a) (b) Figure 2. A raw point cloud and resulting shape of the 3D TIN structure. (a) Prime raw point cloud demonstrated by its class, ground (brownness points), or non-ground (white points). (b) The resulting shape of the 3D TIN construction by using a traditional 3D alpha shape (α = 3).

Improved Ball Pivoting Algorithm based on Spatial Sorting
In the 3D scene, the BPA is usually used to reconstruct the 3D surface of the stereo object represented by point clouds. The failure becomes more likely as the damage of the 3D TIN construction due to unachievable tetrahedron establishment in the bare-ground region (Figure 2(b)). To solve this problem, an improved ball-pivoting algorithm based on spatial sorting was proposed. In this algorithm, we consider three basic types of point, each of which defines a characteristic: seed point, auxiliary point, and candidate points (shown in Figure 3). We now explain these basic types: (1) Seed point is the original ground point for searching nearby points and also the important reference for determining the other points type; (2) candidate point refers to a point which that be selected as a ground point if a certain spatial relationship with a seed point is satisfied; (3) auxiliary point is the third point of the spatial triangle to support the 3-D sphere, which has a minimum distance and the minimum vertical angle relative to the seed point. As shown in Figure 4, the ground points searching strategy can be described as the following steps. The specific implementation procedure is described in Algorithm 1.

Improved Ball Pivoting Algorithm Based on Spatial Sorting
In the 3D scene, the BPA is usually used to reconstruct the 3D surface of the stereo object represented by point clouds. The failure becomes more likely as the damage of the 3D TIN construction due to unachievable tetrahedron establishment in the bare-ground region ( Figure 2b). To solve this problem, an improved ball-pivoting algorithm based on spatial sorting was proposed. In this algorithm, we consider three basic types of point, each of which defines a characteristic: seed point, auxiliary point, and candidate points (shown in Figure 3). We now explain these basic types: (1) Seed point is the original ground point for searching nearby points and also the important reference for determining the other points type; (2) candidate point refers to a point which that be selected as a ground point if a certain spatial relationship with a seed point is satisfied; (3) auxiliary point is the third point of the spatial triangle to support the 3-D sphere, which has a minimum distance and the minimum vertical angle relative to the seed point. As shown in Figure 4, the ground points searching strategy can be described as the following steps. The specific implementation procedure is described in Algorithm 1. Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 16  Output: ground points and none-ground points .
While is not empty: (1) Identify the lowest point or select points with less distance than a threshold to a base surface in the set , add these points to the initial points collection . (1) Identify the lowest point or select points with less distance than a threshold ϕ to a base surface in the set Q, add these points to the initial points collection G.
(2) For iterations G i ∈ G, initialize a set N: Select nearby points by 3D distance of 2r, then add to the set N.
If N is empty, continue the iteration Step (2). Else If N is not empty: (4) Sort points n i ∈ N in ascending order. (5) For iterations m i in the sorted N as an auxiliary point, and temporarily exclude m i from N: (6) Sort N by the following three quantifications in ascending order. Define a BOOL flag _hasFind_ with default value 0. (7) For iterations p i in the resorted N as a candidate point: Build a triangle ϑ by G i , m i , and p i in 3D space; (9) Check the validities of ϑ.
If validities==1: (11) Set the value of _hasFind_ with 1. Remove p i from N, and add m i and p i to the set G, replace m i with p i , and continue the iteration Step (7). (12) End If (13) If _hasFind_== 0 (14) Return m i to the set N. continue the iteration Step(5). (15) End If (16) End If (17) Reture G and E.  Output: ground points and none-ground points .
While is not empty: (1) Identify the lowest point or select points with less distance than a threshold to a base surface in the set , add these points to the initial points collection . A point is chosen as the searching seed in the point dataset. This seed point is the lowest point or with less distance than threshold ϕ to the base surface. Once no further seed point is searched from the primary seed point, this primary seed will be excluded from the ground.
The nearby points are selected by a 3D distance 2r from the seed point (or by the quantity n). If the nearby points are selected successfully, they will be sorted in ascending order by the following two quantifications (Step (4)): (a) The distance d between the seed point and nearby points; (b) the profile angle v between the vertical direction and the direction of the seed point towards nearby points.
The auxiliary point is detected by the minimum distance and the minimum vertical angle. Meanwhile, a candidate points list is sorted by the following three quantifications in ascending order (Step (6)): (a) The horizontal angle h between the direction of the seed point towards the auxiliary point and the direction of the seed point towards nearby points; (b) the distance d between the seed point and nearby points; (c) the profile angle v between the vertical direction and the direction of the seed point towards nearby points.
The spatial triangle is built by the seed point, auxiliary point, and candidate point, and is returned by the iterator to the sorted list. The validity of the spatial triangle can be checked, including (Step (9)): (a) The non-collinearity of the three points; (b) the profile angle θ between the predefined vector and the upward normal vector S of the spatial triangle facet must be less than a given limitation, for instance 55 • ; (c) all other nearby points could be checked if the points are located under the 3D triangle facet inside of the ball; d) no points should be in the half bottom of the ball that was cut in half.
Once the candidate point passes all these validity checks, it will be considered as a new auxiliary point, and this process will return to Step (5). When the process goes through the whole list, all the selected points can be treated as the new seed points and become Step (2) iteratively, until there is no more new seed points.
All these processes and parameters can be represented by this formula: where F s is the whole algorithm, including the two filtering stages. The result ∅ g is the final ground points set. In this algorithm, θ are usually set as constant values. The parameter of r is determined by the largest object above ground in the study area and in terms of double average point-spacing. The parameter of ϕ is dynamic for different terrain environments. The seed point is selected according to whether its distance to the initial surface is lesser or equal to the threshold ϕ, unless there is no initial surface yet. The profile angle can be calculated from the normal direction of one facet in the initial surface, as shown in Figure 3. The spatial triangle consists of the seed point, auxiliary point, and candidate point. Under the spatial triangle, if the dark region formed by 3D spheres and triangles contains no other points, all three points of the spatial triangle are considered ground points, and serve as the seed point for the next search.

Area of Interest and Data Set
The performance of the underside reconstruction filtering was tested by using the Benchmark dataset provided by the International Society for Photogrammetry and Remote Sensing (ISPRS) Working Group III/3 to quantitatively evaluate the performance of different filters and identify developing directions for further research [16]. These samples are located in the Vaihingen/Enz test field in south Germany and Stuttgart city center. They cover diverse feature objects and contain 15 samples, with average point-spacing of 1.0-1.5 m in urban areas (from sample 11 to sample 42) and 2.0-3.5 m in rural areas (from sample 51 to sample 71) [41]. As shown in Table 1, the main terrain characteristics of those samples can be summarized by two groups, named A and B. The outputs, points classified as bare-earth, and object type after the classification process were compared with the true primary classes. All points were labelled as correct bare-earth, type I error, type II error, or correct object. For comparison with other methods, four accuracy indicators for classification algorithms, namely omission error, commission error, total error, and kappa coefficient, were still adopted to access the filtering performance of the proposed method in this paper [9].

Experimental Results
The aforementioned filtering algorithm has been implemented using Visual C++ programming language and all the parameters can be specified by the user on the command line. The sample dataset was divided into two groups in terms of the average point-spacing, as the dense group included samples from sample 11 to sample 42 and the sparse group included samples from sample 51 to sample 71.
The BPA is executed twice in the proposed method. The first time we have filtered the large buildings and objects on the ground. The value of r 1 is set according to the approximate size of the largest non-ground object in the experimental area (r 1 = 20 m, 30 m, 40 m). The second time we have filtered the details on the ground. So, the value of r 2 is set based on the average point distance of the raw point dataset. In samples 11-42, the average point distance is 1 to 1.5 m. Thus, the value of r 2 is twice the middle value of the distance (r 2 = 2.5 m). In a similar way, r 2 is equal to 5 m in samples 51-71, whose average point distance is 2 to 3.5 m. With limited prior knowledge of terrain features for all of these samples, the θ is determined by the maximum angle of topographic inclination (θ = 60 • ) and the constant parameters are fixed as ϕ = 0.5 m (half a meter). The results from this setting are shown in Table 2. From the results in the Table 2 we can see that when we set r 1 = 30 m, the algorithm had the biggest average Kappa coefficient (78.10%) and the lowest average total error (6.94%).
In order to test the parameter of ϕ, the parameters of ϕ from 0.2 to 1.7 with intervals of 0.3 were used to run the proposed filtering process for those samples, and the results of Kappa coefficients and Total errors, in Tables 3 and 4, were employed to demonstrate the effect of ϕ for the algorithm.    As shown in Figure 5, the total means of the Kappa coefficient had a slower increase, while the parameter ϕ was set to be larger and reach its climax after 1.1 m increase of ϕ. Similarly, Group B, which consisted of rough terrain samples, had the same trend in terms of average Kappa coefficient, nevertheless, obtained relative lower accuracies than total accuracies. The difference between total accuracy and the accuracy of Group B could be diminished by increasing the parameter ϕ. On the other hand, Group A, composed of relatively flat samples, with a tiny peak at 0.8 m in terms of accuracy, however, received stable accuracies for all the values of ϕ. The trend of mean Kappa coefficients for corresponding is illustrated by a line chart with a blue line, which is represented as a slower increase for Kappa coefficient with a larger parameter of . The average Kappa coefficients of Group A and Group B are linked by the green line and the red line, respectively.
The detail accuracy information was organized into Table 5, with the parameter combination as r = 30m, = 1.1 , and = 60˚. The errors and the Kappa coefficients are shown in Table 5. The trend of mean Kappa coefficients for corresponding ϕ is illustrated by a line chart with a blue line, which is represented as a slower increase for Kappa coefficient with a larger parameter of ϕ. The average Kappa coefficients of Group A and Group B are linked by the green line and the red line, respectively.
The detail accuracy information was organized into Table 5, with the parameter combination as r 1 = 30 m, ϕ = 1.1 m, and θ = 60 • . The errors and the Kappa coefficients are shown in Table 5. The six samples named 12, 21, 31, 42, 51, and 54 in Group A have an average Kappa coefficient near 90% (89.91%) and an average total error of less than 5% (4.38%). Other samples belonging to Group B had an average kappa coefficient over 70% (73.71%) and the average total error was 7.45%.

Discussion
The distribution of error points derived from the proposed algorithm is illustrated in Figure 6. The blue points represent type I error (omission error: a ground point is misjudged to be a non-ground point.), which there were obviously more of than the red points that represent type II error (commission error; a real ground point is misjudged to be a non-ground point.). As shown in Figure 7, the results in Sample 12, Sample 42, and Sample 71 indicate the high accuracy of our method. In each picture the darker the green, the higher the relative elevation.
In the last decade, several researchers have proposed their new filtering algorithms and used the 15 samples published by ISPRS to evaluate the performances. The average total errors of these filtering algorithms are summarized in Table 6. Without setting intricate thresholds and with no training samples required, the average total error of the proposed method is 6.22%, which is approximately the same as the error in previous studies [23,42], and is 3.48% higher than the lowest one (4.11%).
The distribution of error points derived from the proposed algorithm is illustrated in Figure 6. The blue points represent type I error (omission error: a ground point is misjudged to be a non-ground point.), which there were obviously more of than the red points that represent type II error (commission error; a real ground point is misjudged to be a non-ground point.). As shown in Figure 7, the results in Sample 12, Sample 42, and Sample 71 indicate the high accuracy of our method. In each picture the darker the green, the higher the relative elevation. In the last decade, several researchers have proposed their new filtering algorithms and used the 15 samples published by ISPRS to evaluate the performances. The average total errors of these filtering algorithms are summarized in Table 6. Without setting intricate thresholds and with no training samples required, the average total error of the proposed method is 6.22%, which is approximately the same as the error in previous studies [23,42], and is 3.48% higher than the lowest one (4.11%). Table 6. Comparison of average total errors of the proposed method and eleven other filtering algorithms proposed in the last decade.

Conclusions and Future Work
In this study, a novel filtering method based on the improved ball-pivoting algorithm and spatial sorting is proposed. Compared to previous filtering algorithms, the parameters are easier to comprehend and set. An improved BPA extracted bottom boundary points directly without the process of 3D TIN model construction and effectively solved the problem of output loss in the bare ground region. Furthermore, the performance of this algorithm tends to be stable based on the simple parameter selection method. The experimental results show that this method features reliable performance and robustness. The interpolation error or uncertainty can be effectively avoided on account of this characteristic. In the future work, we will enhance the filtering method   [32] 5.13 Proposed algorithm 6.22 Due to the interference of water area, the lowest points from each grid cell may be clustered and isolated by rivers or lakes. The fundamental step of the improved BPA is spatial searching. If the searching distance is shorter than the size of a river it could miss those low points in other clusters isolated by water. Therefore, to obtain a certain amount of low points in clusters as searching starts in the next step, the lowest points in each grid cell are roughly filtered by improved BPA with a relatively large searching distance. The radius of the pivoting ball is fixed as the average of the width and height of the dataset extent. At this stage, the reference direction for calculating the profile angle was determined as the absolute upward direction with the unit vector V(0, 0, 1). However, in the following stages, this directional vector was specified as the normal vector of the spatial triangle facet constructed in the last stage, which was closest to the relative seed point.
The set of the lowest point in the grid cell may contain those points that are actually building points or vegetation canopy points. Accordingly, the improved BPA based on spatial sorting is applied, beginning with filtered points in stage one, to remove those relatively high points representing building or vegetation at a larger scale. If the estimated length of the largest object in a study area is 1m, the parameter of the ball radius r 1 in this process would be set as 1m (r 1 = 1m). Ground points are selected through the extraction process from the selected seed points in stage one, and then to build the initial TIN surface. For this and the next stage, the profile angle threshold was halved to limit the searching of spatial triangle facets involved with non-ground points.
Another parameter of distance threshold ϕ is used to find seed points with a smaller distance to the nearest facet in the initial surface than the threshold ϕ, or equal to the distance. In one cell, the second improved BPA process is conducted on a relatively small scale. Therefore, the last parameter is introduced as the second ball radius r 2 for the small scale process. The parameter of r 2 is suggested to be set with a value of two or three times the average point-spacing and must be smaller than r 1 . On the woodland, points reflected by ground would be more scattered than non-vegetation areas, thus the parameter r 2 in the condition of dense vegetation would be larger than in non-vegetation environments. Similarly, the parameter of ϕ would be greater under the same conditions. An underlying parameter, which can be determined as two times r 2 , is the size of regular cells with the same width and height. To solve the problem of possible failure of querying no facets in the margin of dataset extent, the reference direction which served as the base direction for searching facets was specified as the unit vector V(0, 0, 1) in those cells which lie at the edge of grid net.
The results in Table 5 showed that the proposed algorithm has a moderate performance in terms of the average total error. In line with the majority algorithms, performances of relatively flat samples received higher accuracy than those in relatively rough or deep areas. While the distance threshold ϕ was increased, relatively flat samples showed tiny difference of performance compared with relatively rough samples. The subtle different trends between the two types of terrain environment may originate from the slope of the ground surface. Owing to those validity rules, above ground points fail to be included into the ground point set beginning with the new selected seed points, which are actually non-ground points, due to a larger distance threshold. Furthermore, a robust performance of total accuracy could be derived from a bigger distance threshold than 0.8 m.
With other unadjusted parameters, the best accuracy of the proposed algorithm occurred at the setting of r 1 = 30 m. This optimal value for the parameter r 1 coincided with the manually measured size of the largest buildings or trees in the studyied area. It is suggested that the application of this proposed algorithm requires the prior acknowledgement of the size of the above-ground objects.
As is the case with any filtering algorithm, there were also some poor results in our experiment. The worst outcomes among Group A and Group B were in Sample 11 and Sample 24 (Figure 7). From these results, we can see that error points tend to occur around steep cliffs.

Conclusions and Future Work
In this study, a novel filtering method based on the improved ball-pivoting algorithm and spatial sorting is proposed. Compared to previous filtering algorithms, the parameters are easier to comprehend and set. An improved BPA extracted bottom boundary points directly without the process of 3D TIN model construction and effectively solved the problem of output loss in the bare ground region. Furthermore, the performance of this algorithm tends to be stable based on the simple parameter selection method. The experimental results show that this method features reliable performance and robustness. The interpolation error or uncertainty can be effectively avoided on account of this characteristic. In the future work, we will enhance the filtering method and further reduce the count of parameters in this algorithm. Some optimizations for algorithmic logic can be adopted to improve the efficiency and the performance of this method as well.