Multidirectional Shift Rasterization (MDSR) Algorithm for Effective Identiﬁcation of Ground in Dense Point Clouds

: With the ever-increasing popularity of unmanned aerial vehicles and other platforms providing dense point clouds, ﬁlters for the identiﬁcation of ground points in such dense clouds are needed. Many ﬁlters have been proposed and are widely used, usually based on the determination of an original surface approximation and subsequent identiﬁcation of points within a predeﬁned distance from such surface. We presented a new ﬁlter, the multidirectional shift rasterization (MDSR) algorithm, which is based on a different principle, i.e., on the identiﬁcation of just the lowest points in individual grid cells, shifting the grid along both the planar axis and subsequent tilting of the entire grid. The principle was presented in detail and both visually and numerically compared with other commonly used ground ﬁlters (PMF, SMRF, CSF, and ATIN) on three sites with different ruggedness and vegetation density. Visually, the MDSR ﬁlter showed the smoothest and thinnest ground proﬁles, with the ATIN the only ﬁlter comparably performing. The same was conﬁrmed when comparing the ground ﬁltered by other ﬁlters with the MDSR-based surface. The goodness of ﬁt with the original cloud is demonstrated by the root mean square deviations (RMSDs) of the points from the original cloud found below the MDSR-generated surface (ranging, depending on the site, between 0.6 and 2.5 cm). In conclusion, this paper introduced a newly developed MDSR ﬁlter that outstandingly performed at all sites, identifying the ground points with great accuracy while ﬁltering out the maximum of vegetation and above-ground points and outperforming the aforementioned widely used ﬁlters. The ﬁlter dilutes the cloud somewhat; in such dense point clouds, however, this can be perceived as a beneﬁt rather than as a disadvantage.


Introduction
Advances in the current measurement techniques have made it possible to routinely use advanced methods of bulk data collection, which found application not only in the technical and engineering fields. Acquisition of such data is possible using 3D scanners (lidars) or photogrammetric systems (usually based on the Structure from Motion method), mounted on a static platform [1], standard ground-or airborne platforms such as a UAV [2][3][4], or even an airship [5]. Such measurement systems result in point clouds of typically hundreds of millions to billions of points. Such large-scale data cannot be processed without at least semi-automated tools. Since at least some parts of the area of interest are massively oversampled, the initial processing phase usually involves dilution to the necessary (while usable) density, which is followed by the extraction of data needed for further processing. These activities typically include vegetation removal and, thus, extraction of points characterizing the objects of interest-be it the terrain, building structures, or anything else.
Such approaches to capturing the reality are being employed in many fields, from reality modeling for the identification of historical monuments [6], infrastructure mapping [7], mapping of open-pit mines [8], monitoring of changes during mining [9], coastline mapping [10], intertidal reef analysis [11], measurements for the purposes of building reconstruction [12], earthworks volume calculations [13], archaeological documentation [14], etc. Often, the data are processed to produce a general [15] or specialized [16] digital terrain model (DTM).

Current Ground Filtering Approaches
Vegetation removal or, more generally, identification of cloud points representing the terrain (ground) is a specific activity for which a considerable number of filters (algorithms and their implementations in software) have been developed. They can be classified into slope-based, interpolation-based, morphology-based, segmentation-based, statistics-based, and hybrid filters; machine learning; and other unclassified filters.
Slope-based filters work with the assumption that the terrain does not abruptly change; these methods determine the slopes in the vicinity of a point and compare the results with a maximum preset slope. These methods work fine in a flat terrain, but they are less effective in a rugged or forested landscape (e.g., [17][18][19]). Interpolation-based filters first select the "seed" ground points from the point cloud using a window size slightly exceeding the size of the biggest object in the area. Then, an approximation of the surface among such points is created and, using a threshold comparing the distance of the remaining points in the cloud from this surface, additional points are added (densify the cloud) (e.g., [20][21][22]). Morphology-based filters are based on the set theory and use operator definition (opening, closing, dilation, and erosion). Mostly, they compare the heights of the neighboring points, working with acceptable elevation differences depending on the horizontal distances of the points. Filtering quality depends on the filter cell size; to improve the filtering quality, morphological filters with gradually increasing cell size and threshold have been designed [23][24][25].
Segmentation-based filters divide the point cloud into segments, and each segment is labeled as either terrain or object; success depends on the quality of the segmentation. They work with point height, intensity, height difference, or geometric attributes (area, perimeter) [26][27][28][29]. Statistic-based filters are based on the fact that terrain points have a normal distribution in terms of height in a certain window, but an object or vegetation points violate this distribution. Based on this assumption, skewness and kurtosis parameters are calculated and used to judge whether a point is on the terrain or on an object [30,31]. These are methods in which no threshold is specified; for example, the method proposed by Bartels (skewness balancing [30]) takes into account the skewness value of the surrounding point cloud and deletes the highest point from the cloud if the skewness is greater than 0. The remaining group of unclassified filters can include, for example, a hybrid filter [32], cloth simulation filter [33], multiclass classification filter using a convolutional network [34,35], or a combination of a cloth simulation filter with progressive TIN densification [36]. Algorithms based on deep learning (neural networks) have also been investigated (e.g., [37][38][39][40][41][42]), but these need to be trained on specific data and are, therefore, probably not yet widely used.
As none of these approaches is universally applicable, the evaluation of the suitability of individual filters for specific purposes has recently been the subject of many scientific papers [42,43] and so has their sensitivity analysis [44].
Most known filtering techniques have been primarily developed for point clouds derived from airborne laser scanning (ALS), the nature of which is significantly different from the data acquired by systems mounted on small aerial vehicles (UAVs) and mobile or static ground platforms (regardless of whether they were mounted with laser scanners or collect imagery for the photogrammetric derivation of the cloud). The latter types of scanning systems typically provide high coverage density and extremely detailed and dense point clouds. Such data have been, in the past, typically predominantly used for relatively flat areas covered with buildings or vegetation. At present, remote sensing is increasingly used for measurements of relatively small areas with rugged rocks, quarries, or similar objects where vegetation removal is highly complicated even for a human operator. Thus, special filters and algorithms that go beyond the "standard" techniques and are capable of utilizing the dense clouds provided, for example, by UAV systems are being developed even for such specific situations [45][46][47].
Based on the above, it can be concluded that the general widely available ground filtering approaches are not suitable for many engineering applications, especially where rugged areas of smaller extent are concerned. Therefore, in this paper, we presented a newly developed filtering algorithm particularly designed for dense point clouds originating from UAV photogrammetry, UAV lidar, or terrestrial laser scanning (with typical densities of hundreds to thousands of points per square meter) in nature-close environments-the MDSR (multidirectional shift rasterization). This algorithm allows the acquisition of a terrain point cloud appropriately describing even the most rugged areas, which does not suffer from problems and artifacts typically found in geometric filters. In this algorithm, the terrain points are individually selected from the cloud as the lowest points of cells of a square grid, which is being gradually shifted in small steps and tilted about all three coordinate axes. In this paper, we aimed to (i) present the principle of the novel filtering algorithm, (ii) test its performance on three areas differing in the degree of ruggedness and vegetation cover, and (iii) compare its performance with that of selected widely used filters.

Materials and Methods
The proposed algorithm differs from the commonly used ones mainly in that it assumes a massive number of redundant points-a situation common in, for example, terrestrial, mobile, or UAV scanning, i.e., scanning with high detail. This assumption contrasts with the assumptions used in the majority of current algorithms (which were created for point clouds acquired by ALS). The algorithm does not aim to select all ground points existing in the point cloud, but rather it aims to identify a sufficient number of such points for accurate characterization of the course of the terrain with emphasis on the reliability of those points that were identified as the terrain.

Algorithm Description
The basic idea of the algorithm is that the lowest point from a certain area is likely to be a terrain point (or, to be more general, a point of the surface of interest covered by vegetation). Obviously, if such an area is square, this is similar to rasterization (but in rasterization, the lowest point detected in the cell is typically not left in its original location but used as a proxy for the whole square and cell and placed in its center). Without any doubt, the likelihood of the lowest point being a ground point increases with the increasing cell size. On the other hand, however, the larger the selected cell, the fewer the points left to characterize the terrain. The basic premise for the choice of the cell size is that it must be larger than the largest object that may lie on the terrain (i.e., a hole in the terrain data caused, for example, by occlusion by vegetation under which there is no point on the terrain, a building, etc.). In the real world, such a cell must necessarily cover an area of at least several meters. Thus, if the point cloud is processed by its division into square cells of sufficient size (in the horizontal plane), the lowest point in each of these cells is very likely to represent the terrain; the number of such points will, however, be very small, and many characteristic points of the terrain (e.g., the highest points of a ridge or pile) will be missing.
The presented method uses two techniques for point cloud densification and capturing the ground points even in the maxima (i.e., points with the highest ground elevation): 1.
Moving the raster as a whole (i.e., shifting the cells in which the lowest point is identified) in steps much smaller than the cell size in both horizontal axes (X, Y); 2.
Tilting (or, rather, rotation) of the raster around all 3 axes (X, Y, Z) changes the perspective for the evaluation of the elevations of individual points in the given projection. After such tilting, the raster is again gradually shifted as in the previous step, and points with the lowest elevation in each cell after each displacement are selected. Figure 1 shows the basic ideas of the algorithm-the idea of shifting and rotating the grid. The figure shows only a cross-section; note, however, that the algorithm itself works in 3D, and, thus, if the cell size is large enough, the algorithm is able to bypass areas that are completely free of ground points (e.g., under a large tree; in such a case, the ground on the edge of the cell is selected, and, in the area without ground points, none are falsely selected).
First, the cloud is divided into cells using a preset cell size ( Figure 1, gray lines), and the lowest point within each cell is detected (Figure 1a, red dots). Subsequently, the grid moves (is shifted) by a preset distance and determines the lowest points in the cells that are newly created in this way, which densifies the cloud (Figure 1b,c). After that, the true innovation of this algorithm is applied-the point cloud is rotated, and the new projection of the grid on the terrain is used to detect the lowest points in thus newly created cells (Figure 1d, blue dots). Then, the grid is again shifted multiple times, and the lowest points after each such shift are determined to further densify the cloud (Figure 1e,f). Figure 1g shows the filtering result. (d-f) rotated raster in 3 steps (i.e., one rotation shifted 3 times)-blue points; (g) resulting terrain points identified in the respective color.
The algorithm in pseudocode is shown in Appendix A and the flowchart in Appendix B. The algorithm is designed to assign each point into a grid cell; if the assigned point is the lowest one in the cell at the time, it is recorded and the previous lowest point in the cell is deleted. In this way, all points are assigned and processed during a single pass for one position (shift and tilt) of the raster, which greatly reduces the computational demand compared with a situation in which all points within a cell would be first identified and the lowest of them subsequently selected. The algorithm is also suitable for easy parallelization; in our tests, 12 processors were simultaneously used.
The computational complexity, therefore, linearly increases with the number of cloud points and quadratically with the number of grid shifts, and, of course, it also increases with the increasing number of preset cloud rotations. The size of the grid (cell) determines the size of the objects that the method is able to bypass and, together with the grid shift size, also the detail of the resulting cloud.
It must be mentioned that if the edge of the cloud contains above-ground points, this thin margin will not be automatically corrected due to the oblique evaluation of the cloud, and it must be removed during further processing (manually or simply by cropping the entire filtered point cloud, which is preferable).

Illustration of the Algorithm Principle
The performance of the algorithm will be demonstrated on a simple example of data-capturing a rugged and heavily forested hill. The point cloud covering an area of 16,700 m 2 was scanned by a UAV-lidar system UAV-DJI L1 (further information can be found in Appendix C); after dilution to a 0.07 m minimum distance between points, it contains 2,995,078 points. The site can be described as a ridge on a flat area (Figure 2), and, besides ruggedness, the data are characterized also by large holes in the terrain data caused by vegetation cover and a small stream running through the area (water absorbs the lidar beam).  As the area size is approximately 140 m × 130 m, the number of ground points will be very small if a 10 × 10 m cell is used. Given the cell size, these points are very likely to be indeed on the ground (and not on vegetation), but in terms of describing the landform, their density is certainly insufficient. Higher densities would be achievable by reducing the cell size, but in such a case, the cell would be likely, in many cases, smaller than the areas without ground data and filtering would be ineffective, identifying vegetation as ground points. Shifting the entire raster can lead to the potential identification of additional ground points. Figure 3b presents the situation after calculations using 5 × 5 displacements in the XxY direction (the grid was, therefore, shifted 5 times by 1/5 of the cell size in each direction). Here, we can observe the desired densification, but the terrain points on the vertices (or on the ridge as here) cannot be, on principle, correctly identified-in the vicinity of the ridge, any 10 × 10 m cell is bound to contain points with elevations lower than the ridge. This problem can be solved by changing the direction of the terrain view, i.e., by rotating the data (see Figure 3c-e-showing data filtered using individual rotations). Rotation also changes the projection of the cell on the terrain and, thus, changes the effective grid size, which also densifies the resulting cloud. Rotations should always be selected in a way corresponding to the estimated slopes of the surface (see below in the Methods). Figure 3f shows the filtering result of the combined data from 18 rotations (a combination of −50 gon, 0 gon, and +50 gon for rotations about the X-and Y-axes; 0 and 50 gon about the Z-axis), each of which was analyzed using 5 × 5 raster positions. Of the original approximately 3 million points, 29,000 remain, with an average density of 1.5 points/m 2 . Figure 3 also demonstrates the presence of gaps in the data; this is most obvious along the ridge (e.g., Figure 3b,e). This gap position changes with individual rotations according to the apparent positions of the ridge relative to the remaining points from the virtual bird's eye view. It is worth noting that the grid shift was also used as a part of the algorithm in [48], but without tilt (rotation), it is not possible to successfully filter areas close to the horizontal and especially the tops (ridges) of landforms or overhangs. The result of filtering without using rotation can be seen in Figure 3b with the already described gaps.

Data for Testing
The algorithm has been tested on point clouds acquired by mobile and terrestrial scanning systems. We present the results from terrains of varying complexity (different from the illustration dataset), from a virtually flat area through an extremely rugged and sloping rock face to a quarry area with a densely forested part of the surface, causing large gaps in terrain points. Lastly, as an example of an area with tall anthropomorphic objects, filtering of a rugged terrain with vegetation and bridge pillars is presented.
Site 1 (Figure 4a)-a virtually flat terrain with low vegetation; an area of 20,000 m 2 ; diluted to 1,408,072 points with a minimum distance between neighboring points of 0.07 m; acquired using a Leica Pegasus mobile scanner mounted on a car (further information can be found in Appendix C). This dataset was selected as a basic one because it can be assumed that virtually any filter should be able to satisfactorily deal with its filtering.
Site 2 ( Figure 4b)-a rural road surrounded by forest on steep slopes; a relatively rugged area of approximately 25,000 m 2 ; 10,016,451 points; diluted to 0.07 m; acquired with a Leica Pegasus mobile scanner mounted on a car (further information can be found in Appendix C). Site 3 ( Figure 4c)-a virtually vertical cliff covered with vegetation, nets, and structures; an area of approximately 1240 m 2 ; diluted to 0.05 m resolution; 2,785,912 points in total; acquired with a Trimble X7 terrestrial scanner (further information can be found in Appendix C). The data from this site contain high-as well as medium-and low-dense vegetation.
Site 4 ( Figure 4d)-an area of approximately 1207 m 2 with rugged terrain, low and tall vegetation, and tall bridge pillars combined with a curved bridge construction; diluted to a resolution of 0.07 m; 830,836 points in total; acquired by the combination of a terrestrial 3D scanner Leica P40 (mainly the construction) and UAV photogrammetry using a DJI Phantom 4 RTK (more information can be found in Appendix C).

Data Processing and Evaluation
The filter presented in this paper was programmed in the Scilab environment interpreted language (www.scilab.org (accessed on 6 May 2022), open-source software for numerical computation; version 6.1.1 was used). Since the latest versions of the environment do not allow parallel executions, the task was pre-divided into the required number of X parts (usually 12 on a 16-processor system) to speed up the computation.
As the data differ in character, various settings were used (detailed in Table 1). Partial rotations about individual axes were used in all combinations. As mentioned above, the estimated slopes of the terrain serve as a basis for the selection of the rotation angles-the angles of rotation should (i) roughly correspond to and (ii) slightly exceed the angles between all normals of the surface and the vertical direction (see the high number of angles needed for site 3 with extreme ruggedness). Note that if the steps are fine in one or two directions, the step in the remaining direction(s) can be somewhat greater as the combination of all rotations already ascertains sufficient detail. As the dataset for site 3 is highly rugged and the surface is very steep, the range of rotations about the X-axis differs from those in other datasets. In site 4, the gaps in ground data (represented by the base of the pillars) are notably bigger than those in other sites; to overcome this, a greater raster size was chosen (7.5 m, i.e., bigger than the gaps), and a greater number of shifts (25 shifts) were used to achieve an appropriate detail.
The evaluation of the quality of MDSR ground filtering using standard methods is difficult in the used complicated datasets. In view of their density and ruggedness, the creation of a reference terrain model that would be used for determining types I and II errors is practically impossible as no perfectly filtered data exist, and even manual assessment would be highly dependent on individual judgment, which is thus associated with high uncertainty. Moreover, the inherent dilution of the cloud would also bias such an evaluation. Another widely used approach using independent control points determined by the GNSS or total station is not suitable here either, as this approach is, in this type of terrain, burdened by the differences between methods (their accuracy, coverage, etc.).
As this study aimed to evaluate ground-filtering quality, this can be only reliably performed using datasets created from the same original data.
For these reasons, the filtering quality was evaluated by comparing with the original cloud and with the results of filtering by widely used conventional filters. First, a visual comparison was performed, and, subsequently, differences between the surface created from MDSR filtering results and those obtained by other filters (and with the lowest points from the original cloud) were calculated.
Choosing the optimal settings for those filters is quite complicated in the rugged areas (as shown, e.g., in [45]). For this reason, multiple settings were operator-tested to identify (for each site and filter) the settings preserving the highest number of terrain points while removing the maximum amount of vegetation points. The selected settings for individual datasets and filters are shown in Table 2.

Visual Evaluation
Visual evaluation of the results of MDSR and those of other filters was performed both on the level of entire surfaces and of profiles from selected potentially problematic areas.

Evaluation of the Differences from the MDSR-Based Surface
This evaluation was performed as follows: MDSR-filtered data were used for the creation of TIN-approximating ground (MDSR-TIN surface), which was cropped at the edges to remove artifacts described above. The signed distance (+ above and − below) of all points from the original cloud as well as from filtered clouds obtained by other filters was calculated, and root mean square deviations of individual points left below and above the filtering result were calculated.
For the comparison with the original (unfiltered) point cloud, only points below the MDSR-TIN surface are of importance (representing points that were omitted in the detection of the ground), whereas points above the MDSR-TIN surface obviously contain all vegetation.
Subsequently, the results of ground filtering by other filters were evaluated in the same way; in these, however, the RMSD for points above the MDSR-TIN surface was also calculated to evaluate the quality of above-ground point removal. The RMSD was calculated using Equation (1): where d is the distance from the created surface approximation (TIN), and n is the number of points.

Results
The above-described filtering techniques were applied on the point clouds; Table 3 shows the number of cloud points before and after filtering using all algorithms.

Visual Evaluation
The first region is practically flat; in principle, any filter should satisfactorily work here. All filters yielded similar results; for this reason, only the result of the MDSR filter is shown ( Figure 5).  Since the quality of ground filtering cannot be distinguished when displaying the entire cloud in a single image, it will be presented on selected profiles showing both the original data and data obtained by the filters under comparison (including the MDSR filter). The profile shown in Figure 6 shows the success of both the geometrical filters and MDSR filters in basic filtering. It is worth noting that in parts of the profile with low vegetation (most likely grass), the MDSR profile is significantly smoother (the profile curve aggregates points from a 0.5 m wide strip of the point cloud, and, hence, it cannot be perfectly thin and smooth).
Clearly, the ATIN and MDSR have the best and very similar results, with MDSR results being slightly more complete. Filtering results of all filters are shown in Appendix D.
The second area captured by data 2 is significantly more rugged, covered with dense, both tall and low, vegetation. Ground filtering in this area is no longer a simple matter, the slopes are quite steep, and the point density rapidly decreases with the distance from the vehicle path. The ground surface data are greatly affected by tree cover. An overall view of the CSF and MDSR filtering results is shown in Figure 7. Results of all filters are shown in Appendix D. After comparing the results, we can conclude practically the same as in the previous case ( Figure 8). The MDSR appears to correctly select ground points even in the areas with very sparse data (left and right edges of the image), with the ATIN performing the second best.
The numbers of points of individual filtered point clouds greatly differ (Table 2). While the original cloud consists of 10 mil. points and the CSF a point cloud of 1.25 mil. points, the MDSR identified only 0.28 mil. ground points. However, as can be obviously seen from the figures, even this is sufficient for terrain description. On the road (center), the point density is approximately 45 points/m 2 ; on the adjacent slopes (outside the covered-up areas), there are about 30 points/m 2 .
It can be, therefore, concluded that the MDSR algorithm in this case successfully filtered the ground and, in addition, suitably and appropriately diluted the point cloud.
Such an approximately 1:5 reduction in the number of points in the cloud will bring a significant benefit for the simplification of further processing while not compromising the terrain description (tens of points per square meter are more than sufficient for a high-quality description of the terrain).
The area where data 3 was acquired is a rugged steep cliff, practically vertical in the upper part; it is, therefore, of high complexity and highly problematic for both capture and subsequent evaluation. In such an area, the filtering results differ a lot. For this site, results from all filters are shown, as in this area, the differences in the filtering success are best visible. The data processed by the CSF are visibly full of holes, especially where the terrain slope is higher or under an overhang (center of the left part in the figure). To capture the ruggedness, shape, and steepness of the area, the grid tilting (explained in Methods) used multiple angles, even as high as 120 deg, which allowed us to capture even the area under the overhang (Figure 9f). The profiles show similar results as for previous sites, but in this case, the differences are even greater (Figure 10). In site 4, the surface is somewhat less rugged than in site 3; however, the tall bridge pillars and supports, the contact of which with the terrain is not just perpendicular, represent the dominant filtering problem in this area. The results of CSF and MDSR filtering are shown in Figure 11a,b, respectively. The results of the remaining filters are similar to those of the CSF and are presented in Appendix D ( Figure A5). The profiles in Figure 12 show that even in this anthropomorphic site, the MDSR filter offers the smoothest and thinnest profile of all tested algorithms. Note that in places where the pillars are in contact with the ground, all the filters identify a part of the pillar surfaces as the ground. Again, the MDSR filter performed the best of all used filters (i.e., identifying the smallest area of the pillar surfaces as the ground-only on a single pillar surface).

Comparison of Point Clouds with an MDSR-Based Surface
The overall RMSDs of points of the original (i.e., unfiltered) point cloud found below the MDSR-generated surface in individual study sites are 0.7 cm, 2.5 cm, 0.6 cm, and 1.7 cm for sites 1, 2, 3, and 4, respectively. This indicates that the MDSR filter indeed selected the lowest points in the cloud. This is comparable with the results of the remaining tested filters, all of which had similar RMSDs of points below the MDSR-generated cloud (Table 4). In other words, all filters were capable of very good detection of the lowest points. However, when comparing the RMSDs of the filtering results taking into account only the points found above the MDSR-generated surface, the RMSDs of the other filters compared with the MDSR ranged from 4.5 cm to 95 cm, indicating a much higher number of unremoved vegetation points. This can be also seen from the thicker profiles of groundfiltered points yielded by the traditional algorithms compared with the new MDSR filter (Figures 6, 8, 10 and 12). Therefore, of all the five tested filters, the MDSR shows the lowest amount of above-ground points (especially unremoved vegetation) that remained in the clouds after filtering.
On the flat surface (site 1), all filters yielded satisfactory results, whereas at sites 2 and 3, only the ATIN filter performed similarly well as the MDSR filter. In site 4, the situation is different due to the large gaps in ground data caused by the pillars-the bigger pillar is 10 m wide and about 5.5 m long at the point of contact with the ground. To successfully perform, the filters including the MDSR need different settings than in the previous scenarios. This particularly affected the ATIN filter; where the "wilderness" setting that was successfully used in the previous sites (with a default step of 3 m) provided better vegetation filtering but above the gaps, parts of unfiltered pillars (up to the height of 14 m) were filtered as the ground. For this reason, the setting had to be changed to "nature" (default step size of 5 m).
The RMSD above the MDSR surface for the original (unfiltered) point cloud basically shows just the average vegetation height.
The results of the ATIN filter implemented in lasground_new software are the most comparable with those obtained by MDSR filtering; this also corresponds to the profiles shown in the previous chapter.

Discussion
Many papers aiming at testing filters and selecting those best suited to the particular purpose have been published. For example, Moudry et al. [42] tested the ground filtering of low-density ALS data capturing a slopy terrain of a spoil heap largely employing the same filters as this paper and comparing the filtering results with control points measured using a total station. Similar to our study, the ATIN also performed very well in their work (root mean square error (RMSE) 0.15 m), followed by the SMRF (0.17 m), PMF (0.18 m), and CSF (0.19 m); nevertheless, the differences between the results of individual filters were relatively small. Štular et al. [49] tested a variety of filters for the determination of the ground from low-density ALS data for archeological purposes. The terrain analyzed in their study was much more rugged with many projections. Again, the ATIN was among the best-performing filters. However, data on the ground filtering of point clouds of a density as high as in our paper capturing rugged areas are practically nonexistent in the literature.
The presented MDSR algorithm works on an entirely different principle compared with not only "classical" geometric filters that create an approximation of the surface using the selected points and then use a threshold to select the terrain points but also other alternative approaches such as that previously presented by our group [45]. Unlike all of these, this algorithm selects the lowest points from the cloud even if only few such points are present, and, as a result of this scarcity of ground points, they fail to form an algorithmically discernible surface necessary for filtering using the above-mentioned approaches.
The main features of our algorithm, therefore, are as follows:

1.
No approximations, simplifications, and assumptions of the terrain are made; 2.
The filtering step also dilutes the point cloud; 3.
Compared with the common geometric filters, this one can be used for much more complex terrains and, therefore, is much more versatile; 4.
The computational demands approximately quadratically increase with increasing required detail (due to the number of raster shifts and raster size); 5.
Similar to all filtering methods, this one also needs verification by an operator; here, typically, a thin layer of filtering artifacts on the edges of the area needs to be removed either manually or simply by cropping by (typically) several decimeters to a few meters; 6.
Where a dense point cloud is needed, the MDSR method could be also used as the first step of an advanced multistep algorithm for the acquisition of the first terrain approximation, after which the remaining points can be identified based on a threshold as in standard filters.
Although this algorithm could work for ALS scanning data, it is primarily designed to work with high-density data of a local character with a high surface density in a complex terrain measured at a short distance (UAV lidar, mobile lidar, terrestrial 3D scanner, SfM method). We should, however, mention a potential pitfall of this algorithm-outliers below the real terrain that the algorithm would detect as terrain points. This is not much of a problem for lidar data used in our paper as during the initial lidar data processing, such outliers are typically removed by the manufacturer-provided software, but this problem might occur when working with SfM data. In such a case, the outliers should be removed by prefiltering using, for example, noise-removal filters.
The algorithm is primarily designed for dense point clouds of nature-close environments (as urban landscapes, typically formed by buildings and planar areas, do not typically need particularly strong vegetation-removal algorithms). However, as illustrated in site 4, which is highly anthropomorphic, the algorithm can work well even in such areas providing that the cell size is set to be larger than the largest gap in the ground points. This, however, has to be compensated by the use of a greater number of stepwise shifts to achieve the terrain model of sufficiently fine detail.
Due to the high computational complexity with a large number of shifts of the base raster, the data can be with advantage gradually reduced. In other words, a smaller window and a smaller number of shifts could be used in the first pass, thus reducing the total number of points (especially above-ground points), and subsequently, a second pass with a larger window and a larger number of shifts could be employed on thus diluted data. Such specific combinations of parameters for stepwise reduction are currently investigated and will be the subject of future papers.

Conclusions
The novel MDSR ground-filtering method primarily designed for nature-close areas presented in this paper was proven to be highly effective. Its testing on four point clouds obtained using various methods at sites with varying ruggedness and vegetation cover and comparison of the performance with other widely used filters (ATIN-Lastools; PMF-PDAL; SMRF-PDAL; CSF-CloudCompare) showed it to be superior to those filters (especially in the terrain with the greatest ruggedness and vegetation cover), with the ATIN performing the best of the rest. This was shown both on the visual presentation of the filtered point clouds and by the evaluation of the distance of the MDSR-based surface from the original (unfiltered) point cloud and from the ground data acquired using the other filters.
The selection of the lowest points in various positions and angles leads, besides the identification of the lowest points in the cloud, also to a dilution of the point cloud. Thus, we can assume that the presented filter can be also used in combined approaches with this filter acting as the first step for creating an approximation of the terrain that can be subsequently densified by the standard threshold-based approach employed in other filters.
However, it should be noted that as the filter is not designed for urban scenes with large roofs, it is not universally applicable, and for such instances, other filters with a changeable window size, such as the SMRF, might be more suitable. Funding: This research was funded by the Grant Agency of CTU in Prague-grant number SGS22/046/ OHK1/1T/11, "Optimization of acquisition and processing of 3D data for purpose of engineering surveying, geodesy in underground spaces and 3D scanning", and by the Technology Agency of the Czech Republic-grant number CK03000168, "Intelligent methods of digital data acquisition and analysis for bridge inspections".

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Appendix A. MDSR Algorithm Code
Below, the algorithm pseudocode is shown. The structure of the data is very simple-it is a matrix where details of every point are in a single line, with the first three numbers always representing the point coordinates (X, Y, Z); additional point properties (such as RGB color components, intensity, etc.) are preserved but not used during filtering.

5.
Saving the identified data-this is solely a routine programming issue of data management function result = SaveLowestPoints(LPointsID, Points) //here, the storage pathways and form are to be programmed. End;

DJI Zenmuse L1 UAV Scanner
For the acquisition of the illustration dataset, the UAV scanner DJI Zenmuse L1 was used. It was mounted on a DJI Matrice 300 quadrocopter, and the basic characteristics of both devices are given below in Tables A1 and A2. Further information can be found on the manufacturer's website (https://www.dji.com, accessed on 6 May 2022).

Leica Pegasus Mobile Scanner
For the acquisition of the point cloud on sites 1 and 2, the mobile system Leica Pegasus Two was used; the basic characteristics of the device are given below in Table A3. Further information can be found on the manufacturer's website (https://leica-geosystems.com/, accessed on 6 May 2022).

Trimble X7 Terrestrial Scanner
For the acquisition of the point cloud on site 3, the 3D terrestrial scanner Trimble X7 was used; the basic characteristics of the device are given below in Table A4. Further information can be found on the manufacturer's website (https://geospatial.trimble.com, accessed on 6 May 2022).

Leica P40 Terrestrial Scanner
For the acquisition of the point cloud on site 4, the 3D terrestrial scanner Leica P40 was used; the basic characteristics of the device are given below in Table A5. Further information can be found on the manufacturer's website (https://leica-geosystems.com, accessed on 6 May 2022).

DJI Phantom 4 RTK
For the acquisition of the point cloud on site 4 also, the DJI Phantom 4 RTK UAV was used; the basic characteristics of the device are given below in Table A6. Further information can be found on the manufacturer's website (https://dji.com, accessed on 6 May 2022).