Adaptive feature-conserving compression for large scale point clouds

Abstract

In this work we introduce a practical method for reducing big point clouds of buildings and infrastructure. The proposed method introduces bilateral filtering with a tailored set of evaluation functions, that will conserve as much information as possible. To determin the actual statistical parameters to perform this filtering and reason about our development, we investigate different point properties on a comprehensive dataset. The dataset contains artificial, photogrammetric and laser scanned point clouds and was made publicly available. We showcase our filtering method by preserving more information than voxel grid or density filters challenging even sparser photogrammetric datasets. Finally, we discuss some encoding strategies as well as the sweet spot between size and resolution.

1. Introduction

2. Related Work

Practitioners from the civil engineering industry mainly have two options when capturing as-is environments as point clouds - laser scanning and photogrammetric reconstruction. While both have advantages and disadvantages (Golparvar-Fard et al. 2011; Baltsavias 1999) , the properties of the resulting point cloud can differ significantly. Laser scanners provide a fixed amount of measures per arc which leads to a radially decreasing point density. Most of the point clouds consist of multiple aligned scans (Boehler et al. 2003) creating a combined but locally inconsistent density: all points are created equal. In contrast, photogrammetric reconstructions are influenced by the number of images, parameter and colour gradients on the surfaces. Some surface materials with restless colouring are reconstructed with higher density and precision, which provides a natural grouping of points in regions of image diversity. Besides this distribution difference, the number of dimensions created by the recording device stands out the most. Original laser scan data consists of the three spatial properties - sometimes with an attached intensity, that encodes the strength of the signal that was the reflected by the hull. Some laser scanners attach colour information with a build-in camera device as a followup. The properties of the photogrammetric reconstruction can differ from one toolchain to another. Usually normals are estimated in the progress of dense reconstruction (Zheng et al. 2014; Schönberger et al. 2016) creating depth and normal maps. All usual properties and their common types are listed in Table 1. Considering these properties, the memory consumption is respectively 28 bytes per point stored by a laser scanner. The same size is taken by several properties of the photogrammetric reconstruction. The size increases during processing due to the optimal Single Instruction Multiple Data (SIMD) memory alignment. Reducing the memory footprint has been in interest of research and work on reducing the consumption of the structure (Elseberg et al. 2013; Schnabel & Klein 2006) . For large point clouds re-encoding the static part in a spatial structure results in huge memory saving but results in an increase of processing time. Another way to reduce the memory consumption is filtering the point cloud. Traditionally the filtering methods are generalized methods, that will be applied to the recording point stream of scanners or simultaneous localization and mapping (SLAM) methods. Most methods such as bilateral or PDE based filtering (Han et al. 2017; Moorfield et al. 2015) try to optimize the cloud with respect to meshing performance. For Airborne Laser Scanning (ALS) a study was undertaken to compare different algorithms (Sithole & Vosselman 2005). The main goal of the examined filter types was not to reduce data, but to increase the signal-to-noise ratio (SNR). Filtering can be archieved by changing single properties of an outlier or by culling (removing) the point in question. Most of the laser scanners, as well as photogrammetric toolchains, employ multiple stages of filtering in the process of creating the merged point cloud. If the original signal is not known, evaluating the SNR can be a challenge. The signal is a set of multidimensional measures along a surface with seemingly non-continuous behaviour. Divergent approaches have been taken to evaluate point clouds: While one approach is the evaluation of suitability in use cases (Rebolj et al. 2017), a second tries to evaluate a certain set of properties and their relations. Example for such analysis are the volumetric density or the number of geometric features (Haala et al. 2013; Angel Alfredo Martell 2017; Dyer 2001) . A recorded set of points can be compared to ground truth data that was generated by a more accurate procedure, such as rapid manufacturing or laser scanning (Tóth et al. 2013). A third approach benchmarks the toolchains and information density (Eickeler et al. 2018) by combining and trading between multiple criteria. As an extension, the cluster density (avg. nearest-neighbour-distance) was measured on multiple spatial cells. This lead to the introduction of density fields and is used to pre-align the clouds prior to a Iterative Closest Point Algorithm (Chen et al. 2017b) . In a follow up the same authors use cluster density as criteria, while examining the deformation and the impact of density on a dataset containing bridges. Leading them to the conclusion that working on dense data rapidly increases the cost of recording and processing, but only adds a small plus at the detection rate (Chen et al. 2018) . Culling measures in an infrastructural context (and therefore addressing the same problem as this work) as well as redefining points by evaluating segments was analysed by Chen (Chen et al. 2017a) and therefore addressing the same problem as this work. In a first partitioning process segments are created and grown. Segments group points with the same direction (eigen values) together and extended by adding points in the vicinity of the point of origin (seed) and completed if a certain threshold is exceeded (Sampath & Shan 2010). For each segment two different key point indicators - smoothness and curvature - are evaluated and a local compression value is calculated. In a consecutive process, a voxel grid is used to reduce the number of measurements in each voxelated segment by this calculated factor. Diverging from this approach, we present a point-wise evaluation of informational gain and reduce the point cloud by grouping and culling points based on our evaluation model. Clusters are formed by similarity and frequency rather than eigen values. 3.

Methods

The basic concept of our method is to evaluate each measurement on a point-by-point basis. We take multiple properties such as point normals, colours, and meta data into account. Based on our preceding analysis of our test dataset, we determin a weighting function and the statistical parameters for each property. These weighting is executed on each point of the cloud and based on their score the point is being kept or removed. For each measurement the score is revaluated as the filtering of the dataset progresses. This ensures that the removed point is always the lowest scoring one.

Taking this concept one step further, we cluster points of similar information content and try to increase the information density. This optimization will drop multiple points if their information content is marginal and add a new important point in between. Our algorithm merges certain points if their weights fall under a certain threshold. The information of the grouping is conserved by encoding the information in the length of the normal vector, without changing the memory layout of the point. In a similar way we compare the weights of a seeded set to their superset and calculate the roughness of the partitions. If the roughness is similar the information of the subset is encoded in the normal.

Datasets. Our test data is constructed of a set of synthetic and real captures. The synthetic consists of geometric primitive and simple construction and point clouds were extracted from multiple meshes. Per square unit we distributed measures with densities of log10 randomly on the containing triangles (Turk 1990) creating 10, 100, …, 1 000 000 [ ⁄ ]. The points spread is random (Mersenne Twister). For the real captures 3 captures were used: A photogrammetric reconstruction of a clock tower at the TU Munich, a photogrammetric reconstruction of a (advanced) construction site at [48°08'50.6"N 11°31'33.4"E] and laser scan of the same construction site that was taken in parallel. The photogrammetric reconstructions were created using the program colmap (Schönberger et al. 2016). The space correlation is established by formula 2. We have chosen a polynomial of 3rd order for distance evaluation and use the euclidean distance as parameter. The use of this 3rd order emphasises sharp angles as well as far distance. The first term, the heaviside function enables thresholding to remove duplicates. All are individual constants for dampening the influence and should be chosen according to the specifications needed (∑ = 1).

The weight evaluates the angle between the normal and the power of 2, we accentuate vivid change. Normal change has a great impact on the surface smoothness and the quality of the Concept. We use bilateral filtering to estimate the information content ( ) for every point.

∈ ( ) = ( , ) + 1 | | ∈ ( , ) Where is the point of interest and part of the set nearest neighbours. The w is a weighting function and ( , ), the eucledian distance between those measures. The last part consists of constant and the sum of the weighted influence of properties of and . The properties are evaluated by their dedicated function . The constant was added since a minimal density is required by multiple guidelines and algorithms (GSA BIM Guide 2019). We propose following weighting functions: ( , ) = ( ( , ) − ℎ ) 1 ( , ) 3 − 2 ( , ) 2 + 3 ( , )

+ 4 1 ( , ) = 2 −1 ⋅

2 | | | |

1 ( , ) = 2

, ( , ) = ( − ) edges. Small normal changes may be due rough materials, greater changes hint interesting geometric features. In theory angles greater that 180 degrees are not covered, but as these angles cannot exists in photogrammetric reconstructions and no patch growing algorithm can generate such normals, this fact can be neglected. The colour weight, is the normed distance in the colour-space. The factor 2 is needed since the colour range is not fixed and determined by the input point cloud. As alternative the intensity can be evaluated, since this measurement must be taken with care as multiple scans can target the same feature but record different intensities. This term is also used for all other properties that were encoded to the cloud. After calculating and removing the point based on the score, the algorithm reevaluate the points in the vicinity. Programmatically this is quite complicated since the lower scores might change drastically, and order can be cached beforehand. For optimizing the calculation, we can use a tree structure and start in multiple cells in parallel, before rearranging the point list. The number of filter cycles should not exceed 50% of the median for dense point clouds and even less for sparse point clouds.

4. Parameter extraction

This framework of formulas was defined by analyzing a broad set of point clouds. For our use we need to populate all statistical parameters and fit the model to our context. We used our full set dataset for the analysis, but only few examples are discussed. Two primitive sets and the high-res construction site are shown to mark the influence of density on the minimal, maximal and average distance (see Figure 1a-c). While on sparse data the reduction of the minimal distance is noticeable, the impact diminishes while using the datasets of medium or high density. This behavior does not differ using NN-2 or NN-10 (data online). For all further evaluation NN-5 was chosen to express a certain robustness and locality. The photogrammetric dataset shows equal properties, which is surprising since the toolchain does not guarantee an even distribution. The density was measured on the roof of the building (see Figure 2c). With the determined polynomial function one can estimate the density and parameters of photogrammetric toolchain before starting with the actual calculation. For further analysis we create histograms of the existing point clouds. For all histograms a double peak can be observed (see Figure 2e-j). The fitted hull curve is based on a double gaussian distribution, and perfectly frames our data. The shape of the distribution does not change with higher density: The first peak can be explained as the influence of the angles on the density, close to the edge the density increasing because of the angular change. If the euclidian distance between the points increases, this influence gets diminished and the histogram of the synthetic examples are uniform. While this seems counter intuitive, these are properties of the nearest neighbour search combined with the influence of planar and almost planar surfaces. In the photogrammetric reconstruction, the peaks melt together and many of the measurements are at lower end. The reason for the faster, rising edge roots in the used reconstruction toolchain, where the maximal resolution is limited by the input image. The flat response of the second gaussian is the outcome of the applied photo consistency function during the dense point cloud creation. As shown in earlier research, the lack of features on rather smooth walls such as concrete tend to have lower density values (Eickeler et al. 2018) . This effect flattens the response and shifts the peak to the right. However interestingly the median and the mean still have a similar distance. For determining good parameters for c1 to c4 we need also to take the contextual properties of our filtering into account. With increasing distance, points become more valuable since we want to guarantee a certain deterministic minimal space distribution. We also want to balance the first gaussian peak as an indicator for geometric turbulence and therefore higher information density. Based on our chosen dataset of this work, we determin: c1 = 8.76, 2 = 9.82, 3 = 3.09, 4 = 0.02. The dampening parameter: = 0.1, = 0.7, = 0.2, = 0.2.

5. Implementation & Evaluation

The filtering was implemented in C++ featuring the Point Cloud Library (PCL) and tinyply. While removing measurements from the cloud the information content needs to be recalculated. Because the needed nearest neighbour search is very costly, a dynamic tree structure is optimal. We were able to achieve similar results and good performance by recalculating a static tree structure after removing a random set of spatial unrelated, lower scoring points. This might not lead to the highest scoring point cloud, but we achieved over 98% of the total score on our test data with lower density. For our real-world data, some colours, intensities and normal were not available and their information context was set to 0. Since we used our own quality criteria to filter the cloud, these same criteria cannot be used to evaluate the performance of our filtering. As an alternative benchmark, we discuss multiple properties of the filtered point clouds to ensure the effectivity of our approach.

In Figure 3a we are visualizing the effect of the filtering to our dense cube. The dense cube was filtered by 80% of the orignial density leaving 120 k measures. The filtering has two properties: (1) The surfaces are reduced to equidistant density. (2) At the edge of the cube the information content is greatest and more points are preserved.

As a showcase of colour filtering we are visualizing the information content of the clock tower data (Figure 3b). This photogrammetric reconstruction has high density and features intresting colouring on the old copper plating. The points are again reduced by 80%, reducing the binary filesize from 260 MB to 60 MB. The figure is split by a fine line, inidicating the crossing of the orignal and the filtered cloud. In the clocktower cloud the texture quality was conserved by the increase in the local density. Overall the original data, the standard deviation of the density was increased by 15%. On features that have both, a normal and a colour change the information content is high. This can be observed at the contact points of the roof plating. The original score of was 12 975 which was improved to 16 778. The filtering of the clocktower without parallelization or optimizations took around 28 min (~10 million points). Further analysis was undertaken to determine the information gain compared to a lower valued reconstruction. We compared the construction site data with two settings of our dense reconstruction pipeline. The underlying geometry and the camera positions were kept constant while the resolution of the depth and normal maps were increased (1000px to 3000px). Additionally, we filtered the high-density dataset to contain the same number of points as the original low-resolution point cloud. The filtering increased the contrast by a vast amount. One example of these results is shown in Figures (14-17).

Grouping: The second part of the filtering is the grouping of points of similar properties. The amount of grouping heavily depends on the grouping parameters provided. By using the weights for grouping the elements below a certain threshold, every parameter that delivers a constant first derivative can be filter points on planes easily. On the point cloud with a higher density the results are basically the same as culling with different parameters. The algorithm behaves as if the eucledian distance in ( , ) to 1.

For the roughness, the distance to the best fitting plane is evaluated. If the superset estimates the same roughness as the subset the points are grouped, and the roughness is saved. We achieved good results by creating a subset from the median. The grouping and the roughness are encoded as the length of the normal vector. We used the first bit to keep the information steady and then encoded the weight by big endian notation in the first 2 bytes and the roughness in the following. While further processing such a point cloud, we faced some problems with existing software. Some of our day-to-day software normalized the vectors before evaluation, others when writing results to a file. In both cases the normals were replaced and the information was lost. The altered point cloud may show different behaviour in a certain set of evaluation algorithms and should not be used for further analysis.

6. Conclusion

We previously showed that we can preserve geometry and colour by selective density variations. Our proposed method shows multiple improvements over the classical filtering approaches of large-scale datasets. Nonetheless employing the filtering increases the complexity of the tool chain and the actual speed up needs to be researched on a case by case basis. While the filtering of small point clouds is relatively fast, its complexity is of ( ) and therefore can get computational expensive. This is the same for most of the algorithms used in the analysis which may be more expensive or less parallelizable. As grouping can be achieved with less extra effort, we would love to recommend it over culling the limitations of the encoding as normal length is not feasible. The determination of the roughness contains the same limitations.

Another impressing fact is that the increase of the minimal resolution is only marginal with the addition of points on surfaces. Adding more points will result in a reduction of error rather than rising the minimal spatial resolution (see figure 1a-c). Practitioners should use historical data or available data sets to evaluate their sweet spot between resolution and effort. As with most recording, if the properties are not completely known, resolution is king. As shown in this work reducing the amount of data in a second step can be done effectively with filtering. However, this does not apply to captures of bigger datasets as the effort is increasing manifold. Before recording, users should ask the question: “what is the statistical error that I can accept” rather than “what are the smaller features that I want to capture”. Naturally the trade of for increasing the resolution of the capture gets worse. The most important part of our filtering model is, that there is no filtering by property. As speciality of bilateral filtering we can apply our method without influencing the latter analysis. This is in stark contrast to the family of PCA based methods, as they will reinforce the later PCA used Recognitions algorithms and change the parameters practitioners are using on smaller, unfiltered samples.

7. Outlook

The focus of the discussion is mainly on the evaluation of speed increase and data reduction opposed to the quality loss. Most operations on points are rather cheap if they are performed on single points but putting them to a bigger scale soon reaches computational limits. One of the limitations is the data throughput and splitting strategies for big data sets. With modern technologies there are certain improvements that should be discussed: First and most importantly, the encoding of the point cloud can be improved. The choice of floats and their non-linear accuracy over their range maybe suboptimal. Another improvement can also be made by encoding the normal vector with radians leaving 4 bytes of information for local functions. All improvements need to be adapted to the availability of extended SIMD instructions which opens the possibility to further enhance the throughput of points. Following up this discussion, a smart process of annotation of points, facets and volumetric elements needs to be developed. The encoding of semantic data on grouped or single points should be researched thoughtfully. Many formats have domain specific solution, such as geodetic information encoded in *.las files, that could be generalized.

Haala, N, Cramer, M & Rothermel, M 2013, 'QUALITY OF 3D POINT CLOUDS FROM HIGHLY OVERLAPPING UAV IMAGERY', ISPRS - International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XL-1/W2, pp. 183–188.

Han, X-F, Jin, JS, Wang, M-J, Jiang, W, Gao, L & Xiao, L 2017, 'A review of algorithms for filtering the 3D point cloud', Signal Processing: Image Communication, vol. 57, pp. 103–112.

Moorfield, B, Haeusler, R & Klette, R 2015, 'Bilateral Filtering of 3D Point Clouds for Refined 3D Roadside Reconstructions' in Computer analysis of images and patterns. 16th international conference, CAIP 2015, Valletta, Malta, September 2 - 4, 2015; proceedings, eds G Azzopardi & N Petkov, Springer, Cham, pp. 394– 402.

Rebolj, D, Pučko, Z, Babič, NČ, Bizjak, M & Mongus, D 2017, 'Point cloud quality requirements for ScanvsBIM based automated construction progress monitoring', Automation in Construction, vol. 84, pp. 323–334. Sampath, A & Shan, J 2010, 'Segmentation and Reconstruction of Polyhedral Building Roofs From Aerial Lidar Point Clouds', IEEE Transactions on Geoscience and Remote Sensing, vol. 48, no. 3, pp. 1554–1567.

Zheng, E, Dunn, E, Jojic, V & Frahm, J-M (eds.) 2014, Patchmatch based joint view selection and depthmap estimation.

References

Angel Alfredo Martell 2017 , Benchmarking structure from motion algorithms with video footage taken from a drone against laser-scanner generated 3D models . Available from: https://robotik.informatik.uniwuerzburg.de/telematics/download/MSc_Angel_Martell. pdf [22 February 2018 ].

Baltsavias , EP 1999 , 'A comparison between photogrammetry and laser scanning' , ISPRS Journal of Photogrammetry and Remote Sensing , vol. 54 , 2 - 3 , pp. 83 - 94 .

Boehler , W , Vicent, MB & Marbs , A 2003 , 'Investigating laser scanner accuracy' , The International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences , vol. 34 , Part 5, pp. 696 - 701 .

Chen , J , Tang, P & Xiong, W 2018 , 'Characterizing Point Cloud Data Density for Spatial Change-Based Maintenance Planning of Civil Infrastructure Systems' . Construction Research Congress 2018 . Construction Information Technology, eds C Wang, C Harper , Y Lee , R Harris & C Berryman , American Society of Civil Engineers, Reston, pp. 776 - 785 .

Chen , J , Zhang, C & Tang, P 2017a , ' Geometry-Based Optimized Point Cloud Compression Methodology for Construction and Infrastructure Management ' in Computing in Civil Engineering 2017 , pp. 377 - 385 .

Chen , M , Wang , S , Wang, M , Wan, Y & He, P 2017b , ' Entropy-Based Registration of Point Clouds Using Terrestrial Laser Scanning and Smartphone GPS' , Sensors (Basel, Switzerland), vol. 17 , no. 1.

Dyer , CR 2001 , 'Volumetric Scene Reconstruction from Multiple Views' in Foundations of Image Understanding , ed LS Davis, Springer US , Boston, MA, pp. 469 - 489 .

Eickeler , F 2019 , Dataset EG -ICE 2019 , Lehrstuhl für Computergestützte Modellierung und Simulation .

Available from: http://vision.cms.bgu.tum.de/eg_ice2019.

Eickeler , F , Borrmann, A & Mistre , A 2018 , 'Defining quality metrics for photogrammetric construction site monitoring' in Proc. of the 12th European Conference on Product and Process Modelling , Copenhagen, Denmark.

Elseberg , J , Borrmann , D & Nüchter , A 2013 , 'One billion points in the cloud-an octree for efficient processing of 3D laser scans' , ISPRS Journal of Photogrammetry and Remote Sensing , vol. 76 , pp. 76 - 88 .

Golparvar-Fard , M , Bohn, J , Teizer , J , Savarese , S & Peña-Mora , F 2011 , 'Evaluation of image-based modeling and laser scanning accuracy for emerging automated performance monitoring techniques', Automation in Construction , vol. 20 , no. 8 , pp. 1143 - 1155 .