Processing outcrop point clouds to 3D rock structure using open source software

This study proposes a processing procedure to extract the 3D rock structure directly from point clouds using open source software. The basic processing involves: (1) estimating the Hough’s normal of each point; (2) converting the normal to dip direction and dip of the corresponding plane; (3) coloring each point in HSV color space according to its normal; (4) decoding the sets number using the multivariate kernel density estimators; (5) extracting and visualizing the set-based points; and (6) estimating the set-based geometric parameters and conducting stereographic projection. The result is an actual discrete fracture network aggregated with the set-based point clouds having HSV colors. From the initial point cloud to the completion of processing, we manage all data in one single file. The case studies show that the processing procedure can identify, extract, and quantify the fracture sets that have less exposed areas, which facilitates the evaluation of main risks.


Introduction
The establishment of overall geometrical configuration of geological fractures in hard rocks is the fundament of a prudent engineering design. The degree of detail required for characterization varies considerably from applications and project stage. For an assessment of the overall rock mass quality, for example, the techniques of scanline and scan window on two-dimensional rock exposures can sample the data such as orientation, spacing and trace length. The estimated fracture parameters, together with intact rock properties, groundwater and stress condition, are the basis for rock mass classifications. In order to evaluate the main risks, it is necessary to consider the three-dimensional spatial organization of fracture sets. For example, in order to maintain the desired excavation geometry during a drill and blast tunneling it is desirable to identify the boundary fractures of removable blocks. The assessment of rockfall hazard in the Alps is an example where the actual fracture pattern of the rockfall source area is indispensable. However, the attempts to investigate the location-dependent rock structure are often confronted with difficult and dangerous accessibility to the rock outcrops, e.g. in an operational environment of tunneling or on rock cliffs situated a few hundred meters above valley floor.
Consequently, the remote collection of the spatial data using geomatics techniques such as LiDAR and digital photogrammetry has expanded rapidly in recent years, which can capture outcrops as point clouds with unprecedented resolution and accuracy. Application examples show that quickly gaining meaningful 3D rock structure data from the clouds including millions of points is still a challenge, especially in daily practices of rock engineering. In this paper, we first summarize the known methods of using point cloud to study the engineering rock structure. Then, we focus on how to use open source

Methods
The method presented here can process any generic point cloud, independent of the data acquisition technology used. If the original point cloud has an RGB code, it will be helpful for the spatial interpretation of the rock structure, but the processing itself does not require a colored point cloud. The entire workflow encompasses the following steps.
First, the point normals of a cloud are calculated using the Hough method in CloudCompare. The normal components (Nx, Ny, Nz) of a point enables the estimation of dip and dip direction of its corresponding plane. The advantage of the Hough algorithm for computing point normal is described in IOP Conf. Series: Earth and Environmental Science 833 (2021) 012054 IOP Publishing doi:10.1088/1755-1315/833/1/012054 3 [9]. As the second step, assign each point an HSV color according to its normal [13]. The HSV colored point cloud can help in visualizing the spatial pattern of fracture sets and identifying the interaction between fractures and the morphology [4]. So far, the processing has extended the original point cloud having only three columns of (x, y, z) coordinate into a big data, in which the newly added columns are the normal (Nx, Ny, Nz) of each point, its HSV color, and dip/dip direction of the corresponding plane. This big data is the input for further processing.
In order to identify the existing set number, we compute the joint distribution of dip direction and dip using the multivariate kernel density estimators bkde2D() of the package 'KernSmooth' and kde2d() of the package 'MASS' in R. We focus here on identifying those sets that have relatively low density but are decisive for rock mass fragmenting. bkde2D() is extremely powerful and can handle millions of points simultaneously. Riquelme [8] used the MATLAB version of the kde2d for recognizing rock joints from 3D point clouds. The design of kde2d() is model-based. Although the treatable number of points is limited, kde2d() can identify the orientation set that has extremely low density. The R algorithm rgl.surface() of the package 'rgl' visualizes the joint distribution as the density field (z) composed by of dip direction (x) and dip (y), which makes the clusters having extremely low density visible. At this point, the step 3 is completed.
The fourth step is to extract the set-based points using the R package 'skmeans' [16][17][18]. This algorithm uses one minus the cosine angle between the vectors as the dissimilarity measure to partition the entire cloud into subclasses, which is what we want to achieve by clustering sets using the normal orientations. The classes number to be entered into skmeans() should be the number of the sets determined in step 3 plus an integer from 1 to 3. This specific value depends on the surface characteristics of an outcrop because in addition to the points of fracture surfaces, there are also slope face, plants or surface deposits in the original point cloud. The silhouette width s(i) of the R package 'cluster' helps in verifying which set of a point belongs to. Points with a large s(i) are very well clustered to its own set, a small s(i) around 0 means that the point lies between two sets, and points with a negative s(i) are probably placed in the wrong set. The new information about the set membership of each point is added to the big data again. In this way, our data management is always concentrated in one R file. For the extraction of a set-based cloud, we only use the points having positive silhouette width. The aggregation of the set-based point clouds results in the in-situ DFN having HSV colors.
The calculation of the set-based orientation parameters follows Scheidegger's eigenvector approach [19]. The used R algorithm is eigen(). The estimated parameters include dip and dip direction of the best mean axis, Fisher's k, and variability as well confidence cones limited with the probability of 0.68, 0.95, and 0.99 for each set. The stereographic projection of fractures as poles can help in visualizing the spatial relationship of the different sets, but using all the millions of points in the projection at one time will make everything unclear. Lyman [20] indicated that if a sample of size 50 is taken from a population of 300 points, the probability that the orientation difference between the sample and the population exceeds 2° is negligibly small. In this paper, the smallest set has about 40000 points, from which we randomly took 30 samples with the sample size of 5000 for calculating the mean orientation and variability. The differences between the orientation estimated using the population (i.e. all 40000 points) and that using the samples did not exceed 1°. Therefore, we randomly selected 5000 points from each extracted set to carry out stereographic projection using the R package 'RFOC'.
Riquelme [21] suggested a method to calculate the spacing of non-persistent discontinuity using 3D point clouds. If the points of a set have been extracted and the set-based mean orientation estimated as introduced above, Stariha [22] shows that the most effective way to measure the set-based spacing is to use the Compass plugin [23], which has been implemented in CloudCompare since 2015. Stariha's approach uses the tool "Measure two-point thickness", with which at least 30 normal spacing values can be directly measured and automatically processed within one minute.
Persistence is the areal extent of a rock discontinuity [24]. It is difficult to quantify the persistence of a fracture in 3D because apart from opaque of rocks, parts of a fracture can also appear in incipient form [25]. Although our quantification on persistence is based on the acquirable area of fractures by IOP Conf. Series: Earth and Environmental Science 833 (2021) 012054 IOP Publishing doi:10.1088/1755-1315/833/1/012054 4 remote sensors, the compilation of the set membership of each point and its HSV rendering has facilitated the establishment of the in-situ DFN.

Application examples
The first example concerns the processing of a limestone outcrop near Moosen in Sankt Johann of Tyrol Austria. We photographed this outcrop with a simple digital SLR camera and reconstructed the 3D scene using open source SfM software [26]. The 3D point cloud has 3.7 million points in a section of 40 m in length and 15 m in height ( figure 1(a)). The processing as described in step1 and 2 shows that there are totally seven orientation groups ( figure 1(b)). The results of clustering, set-based parameter computing, stereographical projection, and examples of the actual DFN as well as spacing are in figure 1(c) to (l). Apart from the bedding (B), there are three joint sets (J1, J2 and J3) ( figure 1(c) -(d)). Figure 1(e) is a zooming-in of the x-shaped intersection of J2 with J1 and J3 that orient in E-ESE. At the same time, the beddings are the top and bottom of the rhombic blocks ( figure 1(f)). Figure 1(g) shows bedding spacing measured using Stariha's method. Figure 1(h) shows the spacing histogram and density distribution of the bedding. The skmeans() divided the slope into lying and hanging parts (figure 1(i) -(j)). The deposits (GS) on the ground and slope surface are a class of its own (figure 1(k) -(l)).
The second example is to show the applicability of our method for processing the point cloud of complex rock structures. The object is the slope face of a hydropower plant in Dachstein limestone, in which the exposed fractures are weathered and eroded by water. Figure 2(a) shows part of this slope. We scanned the whole slope using a Riegl LMS-Z620 scanner in combination with Leica RTK-GPS-Rover for geo-referencing. As an example, we use the point cloud in a section of 60 m in length and 20 m in height merged with the clouds of four scan-positions. At least seven classes are clearly visible in the joint distribution of dip direction and dip ( figure 2(b)). Figure 2(c) shows the HSV-coloring of the cloud. The R algorithms successfully extracted all fracture sets, especially those that cut into the slope and have relatively low point density (figure 2(d)-(e)) by comparison with the huge amount of dense points representing the slope surface (figure 2(f)-(g)). This processing greatly facilitates our judgement on the persistence continuity of the sets having less exposed areas. The algorithms employed also effectively differentiate the areas having strong dissolution (figure 2(h)-(i)) that mainly occurs where slope face and the joints meet. This comprehensive and precise characterization of fracture system provided a reliable basis for us to identify the locally scour-prone blocks.

Conclusion
Geomatics techniques such as LiDAR and digital photogrammetry can acquire engineering rock outcrops as point clouds with unprecedented resolution and accuracy. However, the poorly exposed fractures, although which may be very important for engineering, have a low density in the point cloud. This study proposes a workflow to identify, extract, and quantify such fracture sets that have less exposed areas using open source software. On the outcrops of the two examples in this paper, we also carried out direct contact measurements according to the traditional scanline method. Since the first outcrop is small and not steep, the representative scanlines could be arranged according to our needs. The results of our field scanline mapping confirmed the correctness of the digital approach [26]. However, the digital reconstruction proposed here has a significant advantage over traditional contactbased measurements because of its effectiveness in the comparable time. The outcrop of the second example is a nearly erective slope with a height of more than 20 meters, on which the traditional scanline could only be arranged within a range of 1.5 meters above the ground surface. The orientation sets based on the compass measurements [27] are consistent with those of the algorithm-based clustering introduced in this paper. However, the digital dividing of the orientation sets in this paper is based on the three-dimensional surface measurements of the entire outcrop, which greatly facilitates the spatial extraction of the set-based point clouds. Moreover, the aggregation of the set-based point clouds is the actual discrete fracture network, which is extremely useful for evaluating the engineering behavior of blocky rock masses. Further use and improvement of the introduced procedure may help engineers identify risks at the right time and make more optimal decisions.