Multi-Size Voxel Cube (MSVC) Algorithm – A Novel Method for Dense Point Clouds Terrain Filtering Using Deep Neural Network

1. Introduction

Point clouds are nowadays commonly used to describe object surfaces, including the surface of the Earth. Unlike standard geodetic methods using, e.g., total station or a geodetic GNSS-RTK receiver measuring individual significant terrain points, point clouds are typically non-selectively acquired in quasi-regular grids using, e.g., airborne laser scanning (ALS) [1], terrestrial 3D scanning [2], photogrammetry [3], or mobile laser scanning systems mounted on terrestrial vehicles [4]. More recently, unmanned airborne vehicles (UAVs) equipped with cameras [5] and/or lidar systems [6] have become widely used for this purpose and lately, a great improvement in the quality and reliability of measurements has been observed in mobile (typically operator-carried) SLAM (simultaneous localization and mapping) scanners [7] that can be successfully used even in areas without GNSS signal and with limited visibility (such as underground spaces) [8,9].

All point clouds, regardless of their different character (density, absolute or relative accuracy, origin, etc.), inherently include points that do not represent the surface of interest and need to be filtered out using one or more filtering methods. Extraction of points representing ground (ground filtering, vegetation filtering) is a typical example of such filtering. Many methods and algorithms for this process have been developed, typically based on a single attribute allowing us to distinguish between points of interest and other points. The selection of the filtering method depends, among other things, also on the sensor type; the suitability of individual filters for data acquired using different types of sensors is discussed in systematic reviews [10,11,12].

The most common filters are based on slope (e.g., [13,14,15,16]), interpolation (e.g., [17,18,19,20]), morphology (e.g., [21,22,23]), or segmentation (e.g., [24,25]). Other types of filters include, for example, the statistics-based filter [26], cloth simulation filter [27], a combination of cloth simulation filter with s progressive TIN densification [28], or the MDSR filter based on iterative determination of the lowest terrain points taken from multiple perspectives [29]). In principle, all these methods assume that the terrain is relatively level, and the slope (or change of slope) does not exceed certain maximum values. Based on these assumptions, an approximation of the terrain is created (triangular network, square grid, or cloth simulation) and points until a certain distance from it are considered terrain. However, these methods all come with some limitations, such as difficult distinguishing between rocks (ground) and buildings (non-ground).

Recently, however, the trends in new filtering methods have shifted towards the use of machine learning, which can easily be designed as multicriterial and, therefore, more universal. Here, however, a principal problem occurs – to be able to load an irregular point cloud (note that the points within the cloud are not in any regular grid or structure) into a neural network, it is necessary to transform the cloud into a regular structure.

Machine learning-based methods often employ assessment methods based on local features, i.e., characteristics calculated for each point of the cloud from its close vicinity, such as the variance in spheres of various diameters around the point [30], colors of the points [31], or other spatial and spectral features [32,33,34,35]. Researchers also often use approaches from image analysis, transposing the point cloud into 2D, processing it through standard image analysis methods such as convolutional neural networks (CNNs), and, subsequently, reversely transposing it to a 3D point cloud ([36,37]) from one or more perspectives ([38,39,40,41]). The PointNET [42] architecture and its advanced version PointNet++ [43] have brought about another concept for object recognition in the cloud; it is, however, applicable rather to spatially limited clouds and objects.

Voxelization is another possible approach [44] that could be applied even to large clouds. It is based on the transformation of the irregular cloud into a regular structure; the choice of the voxel size is usually the crucial (and limiting) factor defining the resulting resolution.

As mentioned above, most algorithms have been developed for ALS data and are not entirely suitable for dense point clouds from rugged areas with overhanging rocks and rock faces. Hence, in this paper, we propose a novel method that would not depend on the mathematical definition of the surface or on local features from a small vicinity; rather, this method learns the character of the terrain on a data sample, enabling terrain identification also in other areas. The proposed method is based on the detection of the local cloud shape using a neural network (utilizing, therefore, a simple local feature). The input data for learning are voxelized and evaluated within a 9 x 9 x 9 voxel cube. As the Earth’s surface has an unrestricted character, the optimal voxel size (and voxel cube size) cannot be determined beforehand; it is, therefore, necessary to use multiple consequent voxel sizes in a way ensuring that points with the target attribute (such as terrain, rock face, etc.) are represented in each voxel cube. This method allows us to teach the neural network acceptable patterns of the object of interest with various shapes. We propose a novel method, discuss its pros and cons, and, lastly, evaluate the method's performance through its comparison using a standard and freely available Cloth simulation filter (CSF).

2. Materials and Methods

2.1. Method Principle

Figure 1 provides a 2D illustration of the basic principle of the method, showing the profile of a point cloud in a forested area. The point cloud is voxelized into voxels of suitable size (2x2x2m in Figure 1), each voxel is characterized by its position in the grid and the number of points it contains. Voxels below the terrain contain no points (or very few points, e.g., noise). The evaluation of each individual voxel is performed on the basis of its surroundings (here, 9x9x9 voxels, i.e., the central voxel and four voxels in each principal direction), forming a voxel cube (VC) – the principal assessment unit in this algorithm. To avoid issues associated with differences in point cloud density (caused, e.g., by differences in the distance of the 3D scanner from the objects or height of the unmanned aerial vehicle, UAV, above the terrain), the number of points in each voxel is normalized (divided by the total number of points in the entire voxel cube). In this way, each voxel is represented by the percentage of points in the voxel cube. Subsequently, each individual voxel is classified as ground or non-ground by a neural network trained on a part of the investigated point cloud or on a point cloud from a similar area. The points in voxels classified as ground then enter a second pass with a reduced voxel size, providing a finer detail of the terrain. In this way, the voxel size (and detail of the terrain) gradually decreases until the required detail is achieved (Figure 2, Figure 3 and Figure 4).

The initial voxel size (m) and voxel cube size (n of voxels) need to be chosen in a way ensuring that the voxel cube contains a distinguishable part of the real terrain. For this reason, it is safer to initially select a larger voxel size. A bigger voxel cube size (in the sense of the greater number of included voxels) is capable of providing better terrain identification as it gives the neural network better information on the surroundings; at the same time, however, it increases the computational costs; hence, a compromise size of 9x9x9 voxels was used for the verification of this method and its performance. In our experience, the gradual reduction of voxel size should be to approx. 75 % of that in the previous step, which ensures correct detection even in highly rugged terrain.

The above-described process, however, has an inherent flaw. The red arrows in Figure 3 indicate voxels that include only marginal numbers of terrain points, which will lead to their classification as non-ground. This can be, however, easily prevented by evaluating the point cloud in two voxel grids that are mutually shifted by half of the voxel size in each axis (we will refer to them as the “regular” and “shifted”). Next, all points from the cloud are classified as ground or non-ground based on whether or not they lie within a voxel classified as ground at least in one of the grids.

The classification of narrow tall (or deep) terrain features that could be lost during classification using large voxel size poses another possible problem. To rectify this, the algorithm was adjusted to keep (besides the points in voxels classified as ground) also one additional “envelope” layer of voxels for the next step. This envelope, therefore, helps prevent the undesirable removal of features such as narrow rocks that could be erroneously considered trees. Figure 4 shows the gradual action of the filter in 3D, the classification scheme is visualized in Figure 5.

2.2. Deep Neural Network and Its Training

As mentioned above, the classification itself was performed using a deep neural network (DNN). The inputs into the triangular DNN are the normalized numbers of points in individual voxels, i.e., 729 values (9³) for each voxel cube. The first hidden layer has 1,458 neurons, subsequent hidden layers have 729, 364, 182, 91, 45, and 22 neurons, with one output neuron returning binary 0/1 values (1 – ground; 0 – non-ground). The network was created in Python version 3.11.9. Training and use of the neural network were performed using TensorFlow 2.17.0 libraries and Keras interface. In training, two types of regularization were used – the L2 layer regularization in all hidden layers and drop-out regularization, which has been inserted between each two subsequent hidden layers. The definition of the neural network, including regularization coefficients, is shown in Appendix B. As the terrain shape for each voxel/voxel cube size slightly differs, separate training was performed for each step of the algorithm. To optimize the training process, we started the training from the manually classified data at the finest (required) resolution. As the individual steps of the algorithm use voxels of relatively similar size (75%), the patterns of voxel cubes can be rightly expected to be similar and the DNN training was, therefore, performed stepwise, always using the previous network as the initial state.

For training, data manually classified in CloudCompare v 2.13.0 were used. For each voxel size, the training data are voxelized accordingly and if the particular voxel contains at least one point representing terrain, the voxel is considered to be a ground-representing voxel. For training of each voxel size, the point cloud is cropped to a distance of max. 4 voxels from the nearest ground point (which corresponds to the actual classification process as described above). The training was performed with both regular and shifted grids. Moreover, to allow the best possible training of the general shapes, the DNN was trained also on data rotated by 90°, 180°, and 270°. To account for the differences in the numbers of ground and non-ground voxels in the training data, weights were assigned to the voxel classes to compensate for this.

2.3. Training/Testing Data

2.3.1. Data 1

Data 1 describes a steep slope with rocks. The point cloud was acquired photogrammetrically using the structure from motion (SfM) technique from a manually piloted flight with the UAV DJI Phantom 4. The data was processed in Agisoft Metashape ver. 2.0.1 without any filtering to determine as many terrain points as possible. Considering the very dense vegetation, the resolution was set to ultra-high. Both these settings led to the presence of a considerable amount of noise. No further editing was performed for the purposes of testing. The training data (area of 80x40x40 m; Figure 6 (a) contains 34,859,808 points, the test data (54x32x47 m, Figure 6 (b) consists of 17,278,032 points. In both cases, resolution is approx. 1 cm.

2.3.2. Data 2

Data 2 was acquired in a densely forested area with rugged terrain using a lidar system DJI mounted on a UAV DJI Matrice 300 RTK. For this paper, four discrete areas (see Figure 7, Table 1) were cut out from the scanned area of approx. 340 x 250 m – one was used as a training dataset, the remaining three as test datasets. The reason that we did not use the entire area lays in the very laborious manual preparation of the reference ground surface, the location of the individual data in the area is shown in Figure A1. Data 2 – Boulders (Figure 7(c), Figure 7(d)) show a typical area of this site – relatively rugged terrain with large boulders and the predominance of higher vegetation. The Data 2 – Tower (Figure 7(e), Figure 7(f)) area contains built structures, which were not present in the training data. Data 2 – Rugged (Figure 7(g), Figure 7(h)) show the most challenging area as it is covered with lower dense vegetation (shrubs), which even lidar often fails to penetrate, and contains highly rugged terrain features.

2.4. Testing and Evaluation Procedure

Each used dataset was manually classified into two classes – ground and non-ground. One area on each site was used for the training of the DNN, the remaining ones for testing. The results of testing areas produced by DNNs were evaluated (i) visually to be able to assess the character of the potential classification errors and (ii) using standard accuracy characteristics, i.e., balanced accuracy (BA) and F-score (FS). These two characteristics were used as complementary as BA characterizes the average classification of all (in this case, two) classes considering their representation, while FS rather focuses on the success of classification of the element of interest (in our study, ground). Considering points classified as ground by the algorithm positives (P) and points classified as non-ground as negatives (N), the quality of classification is then expressed as true positive rate (% of all ground points in the cloud that were correctly identified) and true negative rate (% of all negatives in the cloud that were identified as negatives). For details, see the calculation formulas in Table 2 [45].

The voxel sizes for the MSVC method were selected according to the principles described above, with an initial voxel size of 6 m. The minimum (final) voxel size was chosen in accordance with the particular dataset – based on the noise level. It was set to 6 cm for Data 1, and to 11 cm for Data 2, respectively. The entire series of steps in meters was 6.00, 4.50, 3.38, 2.53, 1.90, 1.42, 1.07, 0.80, 0.60, 0.45, 0.34, 0.25, 0.19, 0.14, 0.11, 0.08, and 0.06 m.

To compare the classification success with an established method, we used a freely available Cloth Simulation Filter (CSF) implemented as a plugin in CloudCompare ver. 2.13.0. This filter was selected for comparison because it performed best on rocky terrain in our previous paper in which we tested several widely used freely available filters (CSF, SMRF, PMF) [30]. As there is no universal setting of the CSF, multiple settings were tested for each of the test areas and the best result was always used. Namely, the CSF parameters were set as follows: Scene processing – yes; Scene – Steep slope; Cloth resolution – 0.025; 0.5; 0.1 m and threshold of 0.15; 0.20; and 0.25 m.

3. Results

3.1. Data 1 – Rocks

The results of the principal evaluation, i.e., the TPR, TNR, BA, and FS characteristics, are detailed in Table 3. The MSVC outcome is presented not only for the final (minimal) voxel size, but for the purpose of illustration of the process also for the four previous steps (0.08 – 0.19 m).

While the true positive rate (TPR), i.e., the success rate of identification of ground points, of the CSF was approximately 89%, MSVC reached an excellent TPR of 99.94%, i.e., it classified almost all ground points within the cloud as ground. The TNR (True negative rate, i.e., the success rate of identification of non-ground points) in both methods is approx. 76 %. This is caused by the relatively high density of the non-ground points (low vegetation) in the immediate vicinity of the terrain; hence, neither CSF nor MSVC can fully remove them (due to the CSF threshold and MSVC voxel size). BA of the MSVC method is 88.3% and the FS is 98.3%. This indicates that the ground point classification success is high and the method preserves the necessary points. The CSF results are clearly inferior (BA 82.4%, FS 92.6%), as can be seen in Figure 8 (a) and in the detail in Figure 8 (c).

CSF filter, similar to others tested, e.g., in [29,30] is not ready for this type of terrain. MSVC, however, learns the available terrain shapes, which leads to excellent terrain classification. In the top right corner of the detail in Figure 8 (d), a lower part of a tree (bush) that was misclassified by both algorithms is clearly visible. However, this is a very difficult detail to distinguish from the ground – the vegetation was enclosed in a thin crevice; to be fair to the MSVC, we should note that no such feature was present in the training data for the MSVC filter. In a way, therefore, we cannot consider this to be a failure of the MSVC algorithm, as it “saw” no example of this in the training data; had the training data included such an area, the result would likely be even better.

3.2. Data 2

Due to the use of multiple locations within the Data 2 point cloud, only the best achieved results are presented in Table 4 (complete results are shown in the Appendix C, Table A1, Table A2, Table A3). It should be, however, pointed out that many of the accuracy parameters were very close for multiple settings so it was difficult to pick a single setting in some instances.

In all cases, both BA and FS parameters are better for MSVC than for CSF, although the difference is very low and both methods provide good success. However, the differences are more apparent in data visualizations, see Figure 10, Figure 11 and Figure 12. These figures also obviate that MSVC preserves more points in the rugged areas than CSF, thus better “recognizing” terrain shape. The unfiltered vegetation is low (i.e., very close to the terrain) and is approximately the same in both filtering approaches (Figure 10 (a) and (b)). Unlike in Data 1, this dataset contains also terrain features that were not present in the training data (compare, e.g., Figure 7(c), 7(e), 10, and 11 to Figure 7(a), which shows training data, and Figure 9 showing the buildings in the Tower area). Despite this lack of training data, the MSVC algorithm dealt very well with such areas, better than CSF (the data is more complete).

Data 2 – Boulder (Figure 10) shows a similar rate of misclassification of low vegetation for ground for both algorithms (compare Figure 10(b) to 10(a)). However, as seen in Figure 10 (d), MSVC was more successful in identifying ground points in problematic areas than CSF.

Figure 11 shows results from the Data 2 – Tower area. Both filters suffer from noise more than in Data 2 – Boulder (a, b); more importantly, however, CSF incorrectly classified the roof of the building as ground (a) while MSVC correctly recognized it as a non-ground object despite having no training data on built-up areas. The (c) and (d) panels show a similar effect as in the previous case – the MSVC-generated terrain is more complete (with the exception of the precipice just below the building, which was not correctly recognized as terrain by MSVC).

The results of Data 2 – Rugged classification (Figure 12) show a clear difference between the performance of both filters in the area with dense shrubs (a) and (b), where in such areas, CSF misclassifies low vegetation as terrain more often than the MSVC filter (blue ovals). On the other hand, the yellow oval shows an area of low vegetation with the same shape as the terrain underneath that was incorrectly classified as terrain by MSVC (a), (b). Lastly, the green oval shows an area where MSVC identified significantly more terrain points (see Figure 12(d)), it was, however, at the cost of preserving some points belonging to low vegetation.

4. Discussion

In this paper, we have developed and evaluated a novel ground filter based on a combination of point cloud voxelization and deep neural network. The principle of the method lies in the gradual reduction of the voxel size, which gradually, step by step, refines the terrain.

A vast majority of currently used filters are constructed for airborne laser scanning data, where they serve relatively well – it is, however, necessary to point out that applications utilizing ALS data typically need a resolution of just several points per square meter. However, where dense clouds based on UAV data acquisition (be it SfM-based or UAV-borne lidar clouds) are concerned, these filters may fail to provide satisfactory results, particularly in rugged terrain ([30]). Such dense clouds are often used in engineering applications (such as in quarries ([46,47] or in rocky terrain where landslides and overhanging or collapsing rocks may pose a danger to human constructions or lives [48]).

It should be emphasized that we used data extremely difficult for ground filtering to assess the performance of both algorithms. In ordinary practice, the terrain features that are concentrated in our testing areas represent only a small fraction of the clouds, which leads to excellent overall results even when using any other filter. In other words, if such a difficult terrain as shown in our datasets normally represents a very small proportion of the landscape and the rest is correctly classified by any filter, the error caused by the incorrect classification of these features causes only a negligible overall error.

Although both CSF and MSVC filters performed well in identifying terrain points, the MSVC filter consistently outperformed the widely used CSF filter, especially in the most rugged and problematic areas, such as the rock face in Data 1 (best balanced accuracies of 88.3 % vs 82.4 % and F-scores of 98.3 vs 92.6 % for MSVC and CSF, respectively). As obvious from the true positive rates (which were as high as 99.9 % for MSVC and 89.2 % for CSF, respectively), MSVC correctly detects almost all terrain points, which was not at the cost of decreasing the true negative rate (76.6 % and 75.7 % for MSVC and CSF, respectively); contrary, MSVC performs slightly better even in this parameter. Where Data 2 is concerned, MSVC still outperforms CSF, although the difference is not as high. This is given by the fact that the Data 2 point cloud is more “standard” and contains fewer highly problematic spots. The biggest difference between the two algorithms lies in the better preservation of terrain points by MSVC even in the challenging spots (see the points highlighted in green in Figure 10 (d), Figure 11 (d), and Figure 12 (d). In addition, MSVC is less prone to identifying low vegetation as ground (look at the blue ovals in Figure 13–14).

We were surprised by the fact that unlike CSF, the MSVC algorithm managed to successfully remove the buildings in Data 2 – tower, although no such construction was present in the training data. This suggests that the neural network was capable of a certain abstraction, distinguishing between the true terrain and the roof. The likely explanation is that rather than recognizing the roof and the wall as such (which the algorithm has no way of knowing), it only recognized that these structures do not represent the ground as it did not correspond to any terrain feature it “saw” in the training data. This finding makes the algorithm even more promising as it appears that the filter needs to train only on the character of the ground and the knowledge of the exact vegetation type (or other obstacles) might be less important for its correct function. In effect, this opens the door to the possible more universal applicability of the algorithm – if the algorithm is pre-trained to several terrain types, it might be able to correctly identify terrain even without the need of creating reference data for the particular area; rather, selecting just the terrain type (e.g. flat terrain, rocks, urban, etc.) previously trained could be sufficient. This, however, needs to be verified by future research.

The MSVC algorithm brings multiple benefits to the ground filtering of dense point clouds. It does not need any computationally demanding mathematical operations, the algorithm makes its decisions always solely based on the number of points in individual voxels and its surroundings. Last but not least, the algorithm is inherently robust against the presence of outliers and noise. On the other hand, it is not possible to achieve terrain containing only a single layer of ground points due to the cubic character of the voxel, which always contains remnants of low vegetation. This is, however, true for almost any filter as most of them operate with a parameter such as offset or threshold that characterizes the allowed distance of ground points from the approximated terrain surface.

It is very difficult to compare our results to the existing literature as this is a novel filter that has been never employed before. Moreover, few studies only investigate ground filtering on such dense clouds characterizing highly rugged terrain at resolutions comparable to our study (resulting terrain resolution of approx. 1 to 5 cm).

There is ample space for future research on the use of this algorithm besides the above-mentioned construction of neural network that would be universally trained to detect a particular terrain type (which could simplify the use of this algorithm). The simplest direction of future research is increasing the number of voxels in the voxel cube, which could allow faster operation (need fewer steps). At present, this algorithm only uses a single characteristic – the number of points in the voxel; this could be built upon by the addition of other characteristics, such as the spatial variance. Another possible improvement could lie in training – in the present paper, we have trained the filter on terrain rotated in steps of 90 degrees. Reduction of the rotation step could also lead to improvement of training and, thus, to a better terrain detection performance.

5. Conclusions

In this paper, we proposed a novel ground filtering method utilizing the multi-size voxel cube (MSVC) approach combined with deep neural network. We demonstrated its effectiveness in an extremely difficult (rugged) terrain with dense vegetation. Compared to traditional filters, such as CSF, the MSVC method identified terrain points more accurately thanks to the learning feature of the neural network. The results demonstrate that the MSVC filter can be successfully used for digital terrain extraction from dense clouds even in highly demanding environments where other filters fail, such as steep slopes and/or rugged densely vegetated areas. In the present study, we used a part of the study area for learning, necessitating the manual creation of reference terrain, which can be considered a disadvantage of this algorithm; in the future, however, we aim to produce more generally valid neural networks trained for individual types of landscape (such as rocky mountains, urban environment, hills, etc.), which would make the use of this algorithm much more user-friendly.