IMPROVED TORNADO METHOD FOR GROUND POINT FILTERING FROM LIDAR POINT CLOUDS

: In this paper, an improved Tornado method for filtering LiDAR (Light Detection and Ranging) point clouds is presented. The original method uses a vertical cone with a downward vertex and an upward base to remove the points within it as non-ground points. The remaining points are ground points. The cone moves on the ground surface over the entire region of the point cloud. In this work, the regions of the objects are predicted by extracting the vertical features that have points in the vertical plane or vertical column. Therefore, the tornado method is only used in regions that contain objects. In addition, our improved method uses a specific height for a tornado to reduce the Type I error in mountainous areas. Also, a cylinder surrounding the cone is used to reduce the distance calculations between the cone and the point cloud. The results show that this method is very effective and fast compared to the original method. It also has promising results for the Type I error. In addition, this method was tested on the International Society for Photogrammetry and Remote Sensing (ISPRS) datasets and produced outstanding results. The results show that this method achieves high filtering accuracy. Moreover, the proposed method achieves an overall average error of 6.83%, which is lower than most other methods.

Ground point filtering is an important research field in remote sensing and photogrammetry.It can be done in a variety of methods.These methods can be classified into three categories.The first one is slope based methods (Vosselman, 2000, Sithole and Vosselman, 2001, Meng et al., 2009, Susaki, 2012).These methods are based on the ground surface consisting of several smooth, curved surfaces.Also, the probability of a massive height difference in the terrain is small.If the elevation difference between two adjacent points changes abruptly, the two points could be on different areas of the terrain surface or objects.These methods use a height difference threshold that allows us to determine if the points are ground points or non-ground points by comparing the heights of two nearby points.The second one is related to mathematical morphology methods (Zhang et al., 2003, Arefi and Hahn, 2005, Li et al., 2013, Pingel et al., 2013, Hui et al., 2016).The basic assumption of these methods is to apply morphological operations to the LiDAR point cloud, such as closing, opening, dilation, and erosion, to extract features from raster data.They are used to increase or decrease the size of * Corresponding author objects in raster data.The area type does not affect these methods, but the size of the objects to be filtered does.Consequently, the window size used is important in these methods.The third category is surface methods (Zhang and Lin, 2013, Zhang et al., 2016, Nie et al., 2017, Qin et al., 2017, Shi et al., 2018, Cheng et al., 2019).The concept behind these methods is to use a parametric surface to simulate the bare earth stepwise.The parametric surfaces used in the approaches include the thinplate spline model (Cheng et al., 2019), weighted iterative leastsquares interpolation (Qin et al., 2017), triangulated irregular network (TIN) model (Shi et al., 2018), active shape model (Elmqvist, 2002), and the fabric simulation model (Zhang et al., 2016).Seed points are detected (ground points) and repeatedly compacted to generate an output DEM that continuously refines the ground surface based on specific criteria.The effectiveness of these algorithms depends on the grid size used to select the seed points and the elevation and angle thresholds used to add additional ground points.
In addition to the aforementioned methods, there is an unconventional method that simulates a natural Tornado, which uses a vertical cone to remove object points located within the cone (Mahphood and Arefi, 2020a).This method has been used to obtain good results for flat areas and Type I errors.However, the method has some drawbacks in terms of Type I error in mountainous areas and processing time for large point clouds.
In this paper, an improved Tornado method is proposed to filter the ground points from the LiDAR point cloud.The drawbacks of the Tornado method are eliminated by various improvements.Our contribution to this paper is as follows: 1-Using a specific height of the vertical cone in proportion to the heights of the objects.Thus, it will reduce the number of removed ground points in the steep mountainous areas.2-Using a cylinder that surrounds the cone to decrease the distance computations between the cone and the point cloud.Thus, the processing will be faster.3-Predicting the object regions by extracting the vertical features.4-Moving the cone over the object regions in order not to remove ground points in the mountainous areas.Thus, the processing time and the Type I error will be reduced.

THE PROPOSED METHOD
Figure 1 illustrates the flowchart of the proposed method.
Initially, the noise points are removed.Then, the vertical features are extracted from the point cloud.After that, the vertex points of the cone are extracted from the vertical features.After that, the cone is moved over these points using a cylinder to reduce the computation processing for each cone.Finally, pre-processing step starts to remove the remaining small objects.
Figure 1.Workflow of the proposed method.

Overview and problems
This paper improves the tornado method for ground point filtering (Mahphood and Arefi, 2020a).This method simulates a realistic tornado to remove objects (buildings, trees, etc.) from the earth's surface.A vertical cone models the tornado with a down vertex and an up base (Figure 2).The cone moves on the ground surface using a series of points selected as vertices for the vertical cone.These points are the virtual last pulse (VLP) points (Mahphood and Arefi, 2020b).These pulses are extracted using the virtual first and last pulse (VFLP) method.It extracts the lowest point within a local neighborhood with dres* dres dimensions.The VFLP method resamples the point cloud using a sampling distance dres using equations ( 2) and (3).Thus, the vertex points are distributed over the entire study area.Then, points located inside the cone are removed as non-ground points and the remaining points are ground points (Figure 2(c)).This process is based on the cone equation ( 1).This method has some problems and disadvantages.The use of all vertex points and the filtering procedure have many disadvantages.For using all vertex points, assuming there is an area with no object.Therefore, using vertex points for the entire area will lead to implementing the tornado method in areas that do not contain any object.That is, many cones will be used without obtaining results.Therefore, this will only increase the processing time.On the other hand, in mountainous areas (Figure 3(a)), using cones with great θ angles will lead to removing the ground points and thus increase the Type I error.For handling these problems, the improved method will expect the objects' locations by extracting the vertical features.Then, only the VLP of the vertical features will be extracted.This improvement will be introduced in the tornado motion section.The disadvantage of the filter method is that the method is applied to all points for each cone, which leads to unnecessary calculations of the distance measurement.This only increases the processing time.In addition, in mountainous areas (Figure 3(b)), when cones with large θ-angles are used, the ground points are removed even if the cone is far from the mountains.This increases the Type I error.To deal with these problems, the improved method will use a height threshold for the cone and a cylinder model to constrain the search space.This improvement will be presented in the tornado modeling section.

Pre-processing
In this step, the noise points are removed because they strongly affect the result of the proposed method.This method depends on the points that are inside the imaginary column; therefore, the underground noise points affect the accuracy of the final results if they are not removed.This is because they are selected as vertices, while the ground points are located inside the tornado and are removed as object points.CloudCompare software is used to remove the noise points.It uses the statistical outlier removal (SOR) method, where the distances between neighboring points are averaged to remove the noise points.

Tornado modelling
As mentioned earlier, the tornado is modeled using the cone equation ( 1).Then, the points located inside the cone are removed as non-ground points.For this purpose, two distances are compared.The first is the distance between the point and the cone, and the second is the cone radius at the point height z.The first distance should be equal to or less than the second to remove the point.In other words, equation ( 4) must be satisfied.However, this filtering method has some problems, as seen in the previous section.To deal with the problem in mountainous areas, a height threshold hc is used for the cone so that this height is related to the highest object in the area under study.In other words, we need to know approximately the height of objects, especially buildings and trees.The height threshold hc should be equal to or greater than the maximum object height to ensure that all objects are removed.Therefore, an additional conditional equation ( 5) is added to equation (4) to remove the point.In this way, the mountainous areas will not be affected when the cone is used in the flat areas.This procedure reduces the number of removed ground points in steep mountainous areas.

𝑧 − 𝑧 𝑣 <= ℎ 𝑐 (5)
To handle the search space problem, a cylinder is generated so that it has the maximum cone radius rc related to the cone height hc.This cylinder surrounds the cone (Figure 4(a) and (c)) and has the same vertical axis.It detects the points that should be classified as ground or non-ground points using equations ( 4) and ( 5) (Figure 4(d)).This process reduces the research space for each cone (Figure 4(e)).The radius rc is computed using equation ( 6).The cone vertex is used as the center of the surrounding cylinder to detect the points.Then the points are detected in the 2D space using the radius rc (Figure 4(d)).

Tornado motion
As mentioned earlier, the tornado method uses the VLP as the vertex point for the motion.These points are extracted by VFLP method so that the lowest point is selected for each dres*dres region (Figure 5(a) and (b)).These points are scattered throughout the studied area.The use of these points has problems, as mentioned in section 2.1.To handle these problems, only vertex points related to objects will be extracted.In other words, the virtual last pulse points of the vertical features VLPvf will be extracted.Vertical features have points located in a vertical plane or vertical column (such as trees, walls, etc.) (Mahphood and Arefi, 2020b).Therefore, these points will ensure that the filtering procedure is applied only to objects.This process will reduce the possibility of removing the ground points.The processing speed will also increase significantly because the number of vertex points is much smaller.1-In the first order, the rearrangement is in descending order according to the x-axis.2-In the second order, the rearrangement is in descending order according to the y-axis.3-In the third order, the rearrangement is in descending order (for the VFP matrix) and ascending (for the VLP matrix) according to the z-axis.
Then, equation ( 7) is used.Where pi and pi+1 represent the point and the next point, respectively.Finally, to extract the points of the VFLP, the equations ( 8)-( 10) are applied to detect primary points pi.After that, the following point pi+1 is removed (10) The resulting matrices are the VFP matrix (Figure 6(b)) and the VLP matrix (Figure 6(c)).The count of points in two matrices is the same, and the points have the same order.Finally, the vertical features can be extracted using these matrices so that two height difference thresholds (hmin and hmax) are used to detect the points in this range of height difference.hmin threshold is used to restrict the extracted vertical features by avoiding the low-height features that may relate to the ground (such as small rocks, earth surface ripples, etc.) (equation ( 11)).At the same time, this threshold is related to the lowest objects in the studied area, such as cars and low vegetation.Thus, if there are cars or low vegetation, this threshold is set to 0.5 or 1 m.If there are no low-height objects, it is set to 2.5 m (minimum possible height of buildings).It is worth noting that this parameter can be set to zero.hmax threshold restricts the extracted vertical features by avoiding the unwanted high features (equation ( 12)).This parameter is significant for our method, especially for areas with steep or vertical slopes.These areas will be excluded from the filtering procedure of the tornado method.But this threshold should be equal to or greater than the maximum object height to ensure that all vertical features of the objects will be detected.In other words, this parameter is equal to the cone height threshold: hmax = hc.
where zVFP and zVLP: the z coordinates of VFP and VLP points.
Thus, the VFP and VLP points of the vertical features can be extracted using equations ( 11) and ( 12).The extracted VFP points of the vertical features are denoted as VFPvf (red points in Figure 6(d)), and the extracted VLP points of vertical features are represented as VLPvf (blue points in Figure 6(d)).For our improved tornado method, the VLPvf is used for tornado motion instead of VLP. Figure 6 shows that the VLPvf points are very few compared to the VLP points.This will reduce the time processing and the number of misclassified ground points in the mountainous area.Thus, Type I error will be better.

Post-processing
After implementing the method and getting the results, several sparse points and segments have not been removed in the proposed method.Where a few tiny regions are leftover from plants or objects that aren't classified as objects.As a result, an area threshold Ta is utilized to eliminate these regions.Thus, the region growing method is used to segment the resulting points using a distance threshold r related to the point cloud density of the studied area.

Data Description
The used point cloud is the ISPRS dataset.It is divided into two categories: urban and rural areas and has a wide range of features (open fields, vegetation, small and large buildings, roads, trains, rivers, steep and vertical slopes, low vegetation, gaps, ramp, bridges, discontinuity, cars, etc.).The density in the urban areas is 0.67 points/m 2 , whereas, in the rural area, it is 0.18 points/m 2 , corresponding to point spacing of 1-1.5m and 2-3.5m, respectively.

Evaluation Systems
The qualitative and quantitative evaluations of the results are described in this section.The state of the removed objects will be discussed for the qualitative evaluation.The proposed method evaluates the metrics: Type I (equation ( 13)), Type II (equation ( 14)), and Total error (equation ( 15)) for quantitative evaluation.
Type I =   +  (13) where a is the number of ground points which are correctly detected as ground.b is the number of ground points that are incorrectly detected as an object.c is the number of object points that are incorrectly detected as ground.d is the number of object points that are correctly detected as an object.

Improved Tornado Method evaluation
Results of the improved method with some advantages and disadvantages are shown in Figure 7, Table 1, Table 2, and Table 3.For qualitative evaluation, Figure 7 shows that the proposed method successfully removed all types of objects.It removed the bridges (samples 21, 22, and 71), large buildings (samples 22, 23, and 31), buildings on steep slopes (sample 11), low-resolution buildings (sample 54), vegetation on steep slopes (samples 11 and 51), and other.Also, the gaps and discontinuity are not for our method.From the previous results, we can observe that the improved method achieves promising results.But there are some problems.Some buildings are not removed completely in the areas with steep slopes (squares in Figure 7 (b)) because the roofs are very near the ground surface from one of their sides.Therefore, the buildings will not be located inside the cone and will not be removed.Other disadvantages are shown inside the black circle in Figure 7 (c), where some ground points are removed in the post-processing step.As these points are not connected with the rest of the ground points.Also, from samples 51, 52, and 53, vegetation is not removed completely because these samples have steep slopes and the used angles θ are small, so the ground points are not removed.All samples of the ISPRS datasets are used.These datasets cover different types of areas and terrain features.From Table 1, all the results of the Type I error are perfect.The best result is 0.44% which is related to an area with a steep slope (sample 51).The average of the results is 4.48%.This result is excellent.Most of the results are close to each other.But samp11 has the worst result, which is 12.44%.This result is due to ramps between the ground and buildings (Figure 7 (c)).So, the regions of the ramps were partly removed.Samples 21, 31, 51, 61, and 71 have the best results.However, they contain various features: slight slopes, flat areas, steep slopes, and ramps.The method obtains 21.39% for the average Type II error, which is significantly higher than the Type I error.Sample 53, which contains vertical slopes and ground surface discontinuities, has the worst result, 83 %.Due to the vertical slopes, a small cone angle (20 degrees) is employed.Therefore, the object points are not removed completely.The same problem appears with sample 52.Also, the improved method achieves good results for samples 12, 41, and 42.These samples have different objects: large buildings, low vegetation, vegetation, cars, and others.For the Total error, the average result is 6.83% which is very good.The best result is 2.52%, and the worst is 17.55% which is related to sample 11, which has bad results for Type I error and Type II error.In general, our improved method gives excellent results for Type I error but at the expense of Type II error.Also, the Total error results are very good.

Comparison and Discussion
Our method is proposed to improve the tornado method by reducing the processing time as much as possible.At the same time, it improves Type I errors by reducing the number of the used vertex points.Therefore, in this section, we will compare the results of the original and improved methods in terms of the number of points used as cone vertices, the processing time, and the metrics of the results.At the same time, the proposed method will be compared with other methods.The results of the two methods are shown in To compare the original and improved tornado method in terms of the metric errors, the Type I error, Type II error, and Total error are computed using the same parameters for both methods.The values of the distance dres and the angle θ are shown in Table 2 and Table 3.The results of the errors are listed in Table 1.From Table 1, we find that the proposed method improved all the results of the Type I error.As all results are better than the corresponding ones in the original method.The average Type I error is 4.48% which is less than the average of the original method by 0.66%.For the Type II error, we note that all the results have become worse.These results are due to the fact that all points of the vertical features were not completely extracted.This problem is due to the low density of the used point clouds, which leads to incomplete extraction of the vertical features.Therefore, the objects are not surrounded by sufficient cone vertex points to remove these objects.As a result, the Total error was affected by these results.Therefore, the average of this error is 0.2% less than that of the original method.In any case, it is preferable to use the improved method in the density point clouds.
For comparison of the used vertex points, from Table 2, we find that the number of VLPvf points used as vertex points in the improved method is much less than the number of VLP points used in the original method.This difference is mainly related to the vertical features of the studied area.For example, the difference between the numbers of points is large for samples 53, 61, and 71 because they do not have many vertical features.In contrast, the difference is much less for samples 42 and 54.Also, we can see that the number of the used vertex points has been reduced by at least 90% for samples 53 and 61.Table 2 shows that the best reduction is 94.24%, while the worst is 31.27%.Also, the average reduction is 66.34% which is more than half of the points.As mentioned previously, these results follow the presence of vertical features.The above results certainly affect the processing time.From Table 3, we find that the best improvement is 97%.This improvement is excellent because the processing time has been reduced by 50 times for sample 53.In contrast, the worst is 60%, which is an excellent improvement.From these results, we find that our method gives promising and excellent results in improving processing time.This method is suitable for huge point clouds.
To compare with other methods, Table 5 shows that the Total errors are better than most methods.At the same time, only the Axelsson method is better than ours.Also, our method ranks first for two samples and second for five samples.The average Type I error achieved 4.48% (Table 4).This accuracy is better than the results of other methods.This indicates that the proposed method performs better than the other ones.In addition, our method ranks first for five samples and second for six samples.Unfortunately, from Figure 8, our method does not work well in the areas with ramps, discontinuities, and steep slopes.Where some of the samples (52, 53, and 61) have extremely bad results.Thus, these samples have an impact on the average Type II error.
1) where x, y, and z: the coordinates of the point cloud.xv, yv, and zv: the coordinates of cone vertex.θ: the aperture of the cone.a) The point cloud, b) the vertical cone, and c) removing the points.In (Mahphood and Arefi, 2020a)  res = round(/  ) ×   (2)  res = round(/  ) ×  res (3) where xres and yres: are the resampled coordinates.

Figure 3 .
The disadvantages of using a) all vertex points and b) the cone equation.The red lines are the mountainous areas.The blue lines are the filtered areas.
Figure 4. a) The cone and the surrounding cylinder, b) the cone, c) the point cloud and the surrounding cylinder, d) the points inside the surrounding cylinder, and e) implementing the tornado method with the detected points.

Figure 5 .
Figure 5.The VFLP method.a) The points inside the column before resampling, b) The points inside the column after resampling, and c) the VFP and VLP matrices extraction.In(Mahphood and Arefi, 2020b) To extract the vertical features, the virtual first pulse matrix VFP and the virtual last pulse matrix VLP should be extracted(Mahphood and Arefi, 2020b) (Figure5(c)).In the beginning, the points are resampled using equations (2) and (3) (Figure5(a) and (b)).The resampled coordinates inside the virtual column will have the exact horizontal coordinates (Figure5(b)).Thus, these points will be exactly above each other.After that, the resampled points' matrix is duplicated twice.The first one is for VFP matrix extraction, and the second one is for VLP matrix extraction ∆ ,+1 ≥ 0 for VFP extraction.∆ ,+1 ≤ 0 for VLP extraction.

Figure 6 .
VFP and VLP point extraction.a) The point cloud, b) the VFP points, c) the VLP points, and d) the VFPvf points (red points) and the VLPvf points (blue points).

Figure 7 .
Figure 7.The results of the improved method.a) The samples, b) the filtered ground points, and c) the removed non-ground points.The blue and red points represent the ground points and the non-ground points, respectively.

Table 1 .
Comparison of errors of the ITM with the TM.(The bold text refers to the best values, TM refers to the Tornado method, and ITM refers to the proposed method).

Table 2 .
The number of VLP points and VLPvf points for the samples.

Table 3 .
Time of the processing.