An Iterative Coarse-to-Fine Sub-Sampling Method for Density Reduction of Terrain Point Clouds

: Point clouds obtained from laser scanning techniques are now a standard type of spatial data for characterising terrain surfaces. Some have been shared as open data for free access. A problem with the use of these free point cloud data is that the data density may be more than necessary for a given application, leading to higher computational cost in subsequent data processing and visualisation. In such cases, to make the dense point clouds more manageable, their data density can be reduced. This research proposes a new coarse-to-ﬁne sub-sampling method for reducing point cloud data density, which honours the local surface complexity of a terrain surface. The method proposed is tested using four point clouds representing terrain surfaces with distinct spatial characteristics. The e ﬀ ectiveness of the iterative coarse-to-ﬁne method is evaluated and compared against several benchmarks in the form of typical sub-sampling methods available in open source software for point cloud processing.


Introduction
In recent years, free access to an increasing number of LiDAR (i.e., Light detection and ranging) databases (e.g., the OpenTopography Facility of US and the LiDAR databases from The Environment Agency of the UK) is available to promote the use of LiDAR point clouds. During laser scanning campaigns, the spatial resolution of data acquisition does not often take into account the spatial variation in the surface complexity. In addition, some local areas may be measured multiple times from different stations or flight paths. Consequently, oversampling may occur in some areas and can lead to a large volume of point cloud data, which is more than required for a given application [1][2][3]. This is more likely to be the case for use of open LiDAR data where users have no control on the spatial resolution of data acquisition. In those cases, to make the dataset more manageable and fit for purpose, a thinned point cloud can be sub-sampled from the original point cloud [4,5]. For applications of LiDAR point clouds in Earth science, this may be achieved by selecting a subset of the original point cloud (i.e., the subset selected is the thinned point cloud) or by removing a set of data points from the original dataset (i.e., the remaining data points constitute the thinned point cloud). In either case, the coordinates of individual data points are not altered during the process of thinning from the original point cloud to the reduced one.
The typical sub-sampling methods available in point cloud processing software (e.g., CloudCompare) for data thinning include the random, minimal distance and uniform methods. The random method is the simplest for reducing data density, in which a specified number of data points is selected in a random manner [6,7]. In the minimal distance method, the data point selection is constrained by a minimum distance so that no data point in the selected subset is closer to another data point than the minimum distance specified. In the uniform method, a voxel structure is created and the data point closest to the centre of each voxel is selected. The latter two methods can achieve a more homogeneous spatial distribution of data points in the point cloud of reduced data density [8]. The average data spacing is determined by the minimal distance or the voxel edge length specified. Although these methods are computationally efficient, none of these methods honour the spatial variation of the terrain surface complexity.
From theoretical [9][10][11] and empirical studies [12][13][14][15] on the relation between terrain model accuracy and the density of data points used to build the model, it is widely known that more data points are required to represent local areas where the terrain surface is more complex for a given modelling accuracy requirement (typically, digital elevation model accuracy). To this end, the relative importance of data points needs to be evaluated so that key points are preserved while less important points are removed in the data thinning process [4,5].
A survey of the typical methods for selecting key points was carried out by Heckbert and Garland [16]. Amongst those methods, the point-additive method [5,[17][18][19][20] and the point-subtractive method [18,21,22] are considered suitable and effective methods for scattered point cloud data representing a terrain surface. In the iterative addition method, some initial key points (e.g., local highest or lowest data points [5]) are selected and used to generate a triangulated irregular network (TIN) surface. The deviation of each candidate data point to the TIN surface is calculated globally or locally, and the data point of the largest deviation is classified as a new key point and is added to the key point dataset. The updated key point dataset is used to generate a new TIN surface for selecting the next key point(s). This iterative process continues until a predefined number of data points or a threshold deviation is reached. The point-subtractive method is the reverse of the point-additive method, and starts with a TIN surface based on all candidate data points available. It iteratively removes one or several data points until a predefined number of data points or a threshold error is reached. Although these two types of methods are effective in producing a thinned point cloud with a data density that varies with the terrain surface complexity, such methods often have a high computational cost [5,16,22] because they need to sweep through each and every candidate data point. In addition to the aforementioned methods, sub-sampling with local adaptation to surface roughness can also be achieved using a non-stationary geostatistical approach [23][24][25]. However, this approach requires the production and interpretation of local variograms, which is a time-consuming process and impractical to be applied in practice. A deep learning method for optimised point cloud sampling was also attempted, however it may lead to a thinned point cloud that is not guaranteed to be a subset of the original one [26].
In this research, a spatially non-stationary data sampling method that honours the local surface complexity of a terrain surface was investigated. A new coarse-to-fine resolution method was demonstrated and tested using four LiDAR point clouds that represent bare terrain surfaces of different characteristics. The typical sub-sampling methods were also implemented as benchmarks and compared to the proposed method.

The Coarse-to-Fine Method
The coarse-to-fine method involves an iterative process for sub-sampling from a point cloud, and it includes the following steps (the key steps are also shown in Figure 1). In Step 1, the whole area is partitioned into a set of local sub-areas. In Step 2, regular grid locations G xy over the whole area are defined and allocated to each individual sub-area. In other words, a total number of m k grid locations (i.e., G xy,k ) surrounded by the kth local sub-area are allocated to the sub-area A k . In Step 3, the original point cloud is used to predict the elevation values (z o ) at all the grid locations G xy . In Step 4, a subset point cloud is produced using homogeneous voxel cubes (of edge length S i ) where the data point closest to the centre of each voxel is selected. The method starts with a large edge length S 1 and hence a coarse-resolution subset point cloud. In Step 5, the subset point cloud is used to predict the elevation values (z s ) at all the grid locations G xy . The elevation differences e = z s − z o are calculated and can be readily allocated to each individual sub-area (i.e., based on the allocation of the grid locations in Step 2). These lead to local elevation differences e k,i = z s,k,i − z o,k,i , where z s,k,i and z o,k,i represent the predicted elevations of the ith (1 ≤ i ≤ m k ) grid location within the kth sub-area A k using the subset point cloud and the original point cloud, respectively. The values e k,i are then used to calculate the local metric RMSE k (i.e., root mean square error) for each individual sub-area using the following equation.
where m k represents the number of grid locations in the kth sub-area, e k,i represents the difference between the predicted elevations of the ith grid location within the kth sub-area using the subset point cloud and the original point cloud, respectively.
Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 12 of the grid locations in Step 2). These lead to local elevation differences , = , , − , , , where , , and , , represent the predicted elevations of the i th (1 ) grid location within the k th subarea using the subset point cloud and the original point cloud, respectively. The values , are then used to calculate the local metric RMSE (i.e., root mean square error) for each individual subarea using the following equation.

RMSE = ∑
, where represents the number of grid locations in the k th sub-area, , represents the difference between the predicted elevations of the i th grid location within the k th sub-area using the subset point cloud and the original point cloud, respectively. In Step 6, the following condition is checked: RMSE a predefined threshold . For the sub-areas where the condition is met, the data points (of the subset point cloud) surrounded by those sub-areas are kept. The other sub-areas where the condition is not satisfied are left empty.
Step 6 leads to a sub-sampled point cloud with a combination of filled and empty sub-areas. To fill all empty sub-areas with data points, Steps 4-6 are iterated. In each subsequent iteration, a smaller edge length for the voxel cubes is used, leading to a subset of higher data density. However, the same threshold RMSE is used in all iterations.
The logic underpinning the aforementioned steps is described in the following. The local RMSE is an indicator for the complexity of local spatial variation. If a local RMSE value exceeds the pre- In Step 6, the following condition is checked: RMSE k ≤ a predefined threshold τ RMSE . For the sub-areas where the condition is met, the data points (of the subset point cloud) surrounded by those sub-areas are kept. The other sub-areas where the condition is not satisfied are left empty.
Step 6 leads to a sub-sampled point cloud with a combination of filled and empty sub-areas. To fill all empty sub-areas with data points, Steps 4-6 are iterated. In each subsequent iteration, a smaller edge length S i for the voxel cubes is used, leading to a subset of higher data density. However, the same threshold RMSE is used in all iterations. The logic underpinning the aforementioned steps is described in the following. The local RMSE k is an indicator for the complexity of local spatial variation. If a local RMSE k value exceeds the pre-defined threshold, it indicates that the spatial variation in the local area is more complex than that the threshold RMSE represents. Therefore, more data points are required to preserve the local surface characteristics. This is achieved by reducing the voxel size.

Other Typical Methods for Comparison
Several other typical methods (i.e., the random method, the uniform method and the minimal distance method) are also considered for comparison. Detailed descriptions of these methods can be found in the Introduction.
To evaluate the effectiveness of the methods considered, the thinned point clouds obtained using each method are compared to the original point cloud by interpolation to a total number of N grid locations distributed regularly over the whole area, which leads to the error e j representing the difference between the interpolated values derived using the thinned and the original point clouds at the jth grid location. The values (e j ) at all the grid locations are used to calculate several global metrics, including RMSE, mean error (ME = 1 and maximum deviation (i.e., the maximum value of e j ). Amongst these statistical metrics, the global RMSE is considered the main metric for evaluating the effectiveness of each method and is widely adopted in the relevant literature [4,5].

Study Data
The study datasets considered in this research comprise four airborne LiDAR point clouds. These datasets represent a mix of bare terrain surfaces of distinct characteristics: a hilly and relatively rough surface (Figure 2a defined threshold, it indicates that the spatial variation in the local area is more complex than that the threshold RMSE represents. Therefore, more data points are required to preserve the local surface characteristics. This is achieved by reducing the voxel size.

Other Typical Methods for Comparison
Several other typical methods (i.e., the random method, the uniform method and the minimal distance method) are also considered for comparison. Detailed descriptions of these methods can be found in the Introduction.
To evaluate the effectiveness of the methods considered, the thinned point clouds obtained using each method are compared to the original point cloud by interpolation to a total number of N grid locations distributed regularly over the whole area, which leads to the error representing the difference between the interpolated values derived using the thinned and the original point clouds at the j th grid location. The values ( ) at all the grid locations are used to calculate several global metrics, including RMSE, mean error (ME = ∑ ), standard error (SE = ∑ − ME ), and maximum deviation (i.e., the maximum value of ). Amongst these statistical metrics, the global RMSE is considered the main metric for evaluating the effectiveness of each method and is widely adopted in the relevant literature [4,5].

Study Data
The study datasets considered in this research comprise four airborne LiDAR point clouds. These datasets represent a mix of bare terrain surfaces of distinct characteristics: a hilly and relatively rough surface (Figure 2a

An Example for the Coarse-to-Fine Method
The data processing was carried out in MATLAB. To implement the coarse-to-fine method, an example is presented where specific values of the parameters involved were used. A simple approach for area partition was used, in which the whole area was divided into square sub-areas (e.g., 20 × 20 = 400 square sub-areas) of the same edge length over the whole area. A grid resolution of 1 m was used for the calculation of e k,i . The grid resolution should be fine enough (compared to the size of individual sub-areas) so that an adequate number of e k,i is available for calculating the local RMSE k . For prediction to the grid locations, the TIN with linear interpolation was used. An initial edge length (e.g., S i = 8 m) for the voxel edge length was used. In the following iterations, the edge length S i was gradually reduced at an equal decreasing interval (e.g., 0.2 m). The ending edge length depends on a user-specified threshold RMSE.
To show a sub-sampled point cloud using the coarse-to-fine method, an example is presented using the point cloud representing the hilly rough surface and a user-specified threshold RMSE of 0.085 m. The values of the parameters applied were those reported in the previous paragraph. The original and the sub-sampled point clouds are shown in Figure 3a,b, respectively. The sub-sampled point cloud consists of 67,457 data points (approximately 22.67% of the original dataset), and its spatial variation in data density is obvious. For the MATLAB (2017a version) codes used, the time taken for the subsampling was approximately 58 seconds using an average office computer (4 cores of i5-3470 CPU @ 3.2 GHz and 8 GB RAM).

An Example for the Coarse-to-Fine Method
The data processing was carried out in MATLAB. To implement the coarse-to-fine method, an example is presented where specific values of the parameters involved were used. A simple approach for area partition was used, in which the whole area was divided into square sub-areas (e.g., 20 × 20 = 400 square sub-areas) of the same edge length over the whole area. A grid resolution of 1 m was used for the calculation of , . The grid resolution should be fine enough (compared to the size of individual sub-areas) so that an adequate number of , is available for calculating the local RMSE . For prediction to the grid locations, the TIN with linear interpolation was used. An initial edge length (e.g., Si = 8 m) for the voxel edge length was used. In the following iterations, the edge length Si was gradually reduced at an equal decreasing interval (e.g., 0.2 m). The ending edge length depends on a user-specified threshold RMSE.
To show a sub-sampled point cloud using the coarse-to-fine method, an example is presented using the point cloud representing the hilly rough surface and a user-specified threshold RMSE of 0.085 m. The values of the parameters applied were those reported in the previous paragraph. The original and the sub-sampled point clouds are shown in Figure 3a and Figure 3b, respectively. The sub-sampled point cloud consists of 67,457 data points (approximately 22.67% of the original dataset), and its spatial variation in data density is obvious. For the MATLAB (2017a version) codes used, the time taken for the subsampling was approximately 58 seconds using an average office computer (4 cores of i5-3470 CPU @ 3.2 GHz and 8 GB RAM).
It is expected that more data points should be selected in sub-areas where the surface is more complex. To confirm this, a local surface roughness map was produced using the local root mean square height (RMSH) method, which is a roughness metric applied widely in Earth Science [27,28]. The window size used for calculating the RMSH was the same as that for the individual sub-areas for data sampling. The local surface roughness (shown in Figure 4b) has an evident correlation with the local data density (in Figure 4a).  It is expected that more data points should be selected in sub-areas where the surface is more complex. To confirm this, a local surface roughness map was produced using the local root mean square height (RMSH) method, which is a roughness metric applied widely in Earth Science [27,28]. The window size used for calculating the RMSH was the same as that for the individual sub-areas for data sampling. The local surface roughness (shown in Figure 4b) has an evident correlation with the local data density (in Figure 4a).

Comparison on Effectiveness
The commonly used methods and the coarse-to-fine method were implemented for the four LiDAR datasets considered. The original point clouds were reduced to lower data densities of various levels (represented by the number of data points selected). The level of reduction ranges from approximately 4% to 80% of the original point clouds.
The thinned point clouds were compared to the original point cloud at a grid resolution of 1 m. The statistical metrics introduced in Section 2.2 were calculated. For a more direct comparison, the parameters for the minimal distance method and the uniform method were determined using a trialand-error approach so that the number of data points selected using these two methods were approximately the same as those selected using the coarse-to-fine method. For the random method, multiple rounds (i.e., 20 in this study) of random data selection were performed, and the results were averaged over all rounds to produce the mean values of the metrics considered.
The results in terms of RMSE are shown in Figure 5. It is expected that the RMSE decreased with increasing data density for the sub-sampling methods and the terrain surfaces investigated. The random method produced the lowest accuracy for all the datasets considered. Similar accuracy (higher than that of the random method) was observed for the minimal distance method and the uniform method, suggesting that evenly distributed data points lead to smaller differences from the original point cloud than a random distribution. The coarse-to-fine method performs more effectively than the other methods that do not honour the spatial variation of the terrain surface characteristics.

Comparison on Effectiveness
The commonly used methods and the coarse-to-fine method were implemented for the four LiDAR datasets considered. The original point clouds were reduced to lower data densities of various levels (represented by the number of data points selected). The level of reduction ranges from approximately 4% to 80% of the original point clouds.
The thinned point clouds were compared to the original point cloud at a grid resolution of 1 m. The statistical metrics introduced in Section 2.2 were calculated. For a more direct comparison, the parameters for the minimal distance method and the uniform method were determined using a trial-and-error approach so that the number of data points selected using these two methods were approximately the same as those selected using the coarse-to-fine method. For the random method, multiple rounds (i.e., 20 in this study) of random data selection were performed, and the results were averaged over all rounds to produce the mean values of the metrics considered.
The results in terms of RMSE are shown in Figure 5. It is expected that the RMSE decreased with increasing data density for the sub-sampling methods and the terrain surfaces investigated. The random method produced the lowest accuracy for all the datasets considered. Similar accuracy (higher than that of the random method) was observed for the minimal distance method and the uniform method, suggesting that evenly distributed data points lead to smaller differences from the original point cloud than a random distribution. The coarse-to-fine method performs more effectively than the other methods that do not honour the spatial variation of the terrain surface characteristics. Figure 6 shows very small mean errors as compared to the corresponding RMSE, suggesting that systematic (or biased) deviations of the thinned point clouds from the original point clouds were minimal for all the sub-sampling methods considered. It was also observed that a more complex surface was likely to result in a larger mean error (e.g., Hilly Rough versus Flat Smooth). The standard errors are shown in Figure 7, which illustrates the methods had similar results to those based on the RMSE metric. This is expected as the mean errors were all comparatively small. Figure 8 shows the maximum elevation difference between the thinned point cloud and the original point at the pre-defined grid locations. Smaller values were observed for the coarse-to-fine method for the four different types of terrain surfaces considered. The maximum deviation decreased with increasing data density, but those changes were not in a monotonic fashion (especially for the minimal method and the uniform method).  Figure 6 shows very small mean errors as compared to the corresponding RMSE, suggesting that systematic (or biased) deviations of the thinned point clouds from the original point clouds were minimal for all the sub-sampling methods considered. It was also observed that a more complex surface was likely to result in a larger mean error (e.g., Hilly Rough versus Flat Smooth). The standard errors are shown in Figure 7, which illustrates the methods had similar results to those based on the RMSE metric. This is expected as the mean errors were all comparatively small. Figure 8 shows the maximum elevation difference between the thinned point cloud and the original point at the predefined grid locations. Smaller values were observed for the coarse-to-fine method for the four different types of terrain surfaces considered. The maximum deviation decreased with increasing data density, but those changes were not in a monotonic fashion (especially for the minimal method and the uniform method).

Discussion
The minimal distance method and the uniform method are algorithmically different, leading to different sub-sampled point clouds. However, these two methods tend to produce a sub-sample of

Discussion
The minimal distance method and the uniform method are algorithmically different, leading to different sub-sampled point clouds. However, these two methods tend to produce a sub-sample of largely homogenous distribution. This is likely the reason for similar RMSE values shown in Figure 4. Therefore, either method can be considered if one wishes to obtain an evenly distributed sub-sample.
The coarse-to-fine method is developed for sub-sampling point clouds that represent terrain surfaces. Although it was tested using airborne LiDAR-derived terrain point clouds, it is expected that it can readily be applied to terrain point clouds obtained using other surveying techniques. The coarse-to-fine method is able to produce a thinned point cloud that is adapted to local terrain surface roughness and can lead to more uniform local RMSE values over the whole area of interest.
For a given threshold RMSE in the coarse-to-fine method, the set of data points selected is affected by (1) the interval (e.g., 0.2 m used in Section 4.1) of the voxel edge length decreasing from a large value (i.e., coarse-resolution) to a small value (i.e., fine-resolution), and (2) the size of individual sub-areas. When a smaller interval is used, more varied data density bands are available. This will lead to a sub-sampled point cloud that has a smaller difference to the original point cloud, which is confirmed by the analysis results shown in Figure 9a. When a smaller sub-area size (i.e., more sub-areas) is used, the local surface roughness or complexity is better represented by the finer resolution, which also leads to a more favourable sub-sampled point cloud (Figure 9b). However, the increase in accuracy is less significant when the number of sub-areas used is large (e.g., from 400 to 900 sub-areas). The coarse-to-fine method is area (or block) based. Within a single local area, the relative importance of data points is not evaluated. The data density is adapted to only the between-block variation in surface roughness (often calculated over a particular scale), which can lead to an abrupt change in data density at the edges of the blocks. This issue may be mitigated by using a smaller block size and a smaller decrement internal for density reduction. In this study, the TIN with linear interpolation is used for prediction. However, the coarse-to-fine method can be readily implemented using other interpolation methods.
As shown in Figure 4, the spatial variation of data density of the thinned point cloud has a large correlation with that of the surface roughness map, and as such it may be used to infer the local roughness of the terrain surface.
A simple area partition method (i.e., square blocks) was used in this study. As a suggestion for further work, it would be interesting to investigate the performance of the coarse-to-fine method when more complex area partition methods (e.g., over-segmentation approaches such as voxel cloud connectivity segmentation and point cloud local variation [29][30][31][32][33]) were used. Such methods may take into account the terrain surface complexity for the area partition. The coarse-to-fine method is area (or block) based. Within a single local area, the relative importance of data points is not evaluated. The data density is adapted to only the between-block variation in surface roughness (often calculated over a particular scale), which can lead to an abrupt change in data density at the edges of the blocks. This issue may be mitigated by using a smaller block size and a smaller decrement internal for density reduction. In this study, the TIN with linear interpolation is used for prediction. However, the coarse-to-fine method can be readily implemented using other interpolation methods.

Conclusions
As shown in Figure 4, the spatial variation of data density of the thinned point cloud has a large correlation with that of the surface roughness map, and as such it may be used to infer the local roughness of the terrain surface.
A simple area partition method (i.e., square blocks) was used in this study. As a suggestion for further work, it would be interesting to investigate the performance of the coarse-to-fine method when more complex area partition methods (e.g., over-segmentation approaches such as voxel cloud connectivity segmentation and point cloud local variation [29][30][31][32][33]) were used. Such methods may take into account the terrain surface complexity for the area partition.

Conclusions
In this article, a simple non-stationary data sampling method for reducing point cloud density was investigated. The application of this method to the LiDAR datasets confirms its ability to produce a thinned point cloud with a data density (in sub-areas) which is adapted to local terrain surface roughness. For the different types of terrain surfaces, it was confirmed that the coarse-to-fine method performed more effectively than several other sub-sampling methods considered.