Surface Reconstruction for Three-Dimensional Rockfall Volumetric Analysis

Abstract

Laser scanning is routinely being used for the characterization and management of rockfall hazards. A key component of many studies is the ability to use the high-resolution topographic datasets for detailed volume estimates. 2.5-Dimensional (2.5D) approaches exist to estimate the volume of rockfall events; however these approaches require rasterization of the point cloud. These 2.5D volume estimates are therefore sensitive to picking an appropriate cell size to preserve resolution while minimizing interpolation, especially for lower volume rockfall events. To overcome the limitations of working with 2.5D raster datasets, surface reconstruction methods originating from the field of computational geometry can be implemented to assess the volume of rockfalls in 3D. In this technical note, the authors address the methods and implications of how the surface of 3D rockfall objects, derived from sequential terrestrial laser scans (TLS), are reconstructed for volumetric analysis. The Power Crust, Convex Hull and Alpha-shape algorithms are implemented to reconstruct a synthetic rockfall object generated in Houdini, a procedural modeling and animation software package. The reconstruction algorithms are also implemented for a selection of three rockfall cases studies which occurred in the White Canyon, British Columbia, Canada. The authors find that there is a trade-off between accurate surface topology reconstruction and ensuring the mesh is watertight manifold; which is required for accurate volumetric estimates. Power Crust is shown to be the most robust algorithm, however, the iterative Alpha-shape approach introduced in the study is also shown to find a balance between hole-filling and loss of detail.

1. Introduction

Over the past decade, terrestrial laser scanning (TLS) has proven to be an efficient and effective tool for slope stability assessment and monitoring. TLS systems have been used to investigate fundamental rockfall processes. Multi-temporal TLS scans have been used to identify and characterize patterns of precursor rockfall activity [1,2,3]; and to detect and analyze millimeter-scale pre-failure deformation in rock slopes [4,5]. TLS systems have also been used for rock bridge characterization [6], assessing cliff-talus interactions [7], and most commonly; the development of remotely sensed rockfall databases [8,9,10,11,12,13].

Many of these studies have used detailed multi-temporal TLS point clouds or surface models to derive a frequency-magnitude relationship for the study site. All of these studies are hindered by the temporal longevity that is required to generate realistic frequency-magnitude relationships. Work by Williams et al. [12] demonstrates the effects of temporal censoring by TLS acquisition rate and its implications on frequency-magnitude relationships. Regardless of the time period the slope is being monitored and the data acquisition rate, these analyses rely on confidence in the volumetric calculations to derive frequency-magnitude relationships.

Recent work by Olsen et al. [10] raises a key issue with many of these analyses. At many study sites, complex slope geometry results in occlusions (i.e., holes) in the point cloud captured by the TLS system. Additionally, the divergence and attenuation of the laser beam with increasing range can result in non-uniform point density across the study slope. Both factors can present challenges for change detection analyses and subsequently, volumetric analysis. To compute the volume of detected topographical change between sequential scans, the areas of change are commonly converted from a point cloud to a triangulated surface mesh. However, many surface reconstruction techniques were not developed for natural environments. Some algorithms assume uniform point spacing and smooth surface geometry which might not present in the natural environment [14]. As a result, the conversion from a point cloud to a mesh can result in topological artifacts which inhibit accurate volumetric calculations. For accurate volume calculations, the triangulated surface mesh must be watertight, free of intersecting triangles and have consistent normal vector orientation. As noted by Olsen et al. [10], ensuring that each mesh is topologically correct can require significant manual editing which is time-consuming and subjective. Therefore there is a need to evaluate automated surface reconstruction algorithms and their influence on resulting volumetric analyses, in order to improve methods used to monitor and quantify geomorphological processes.

In this work, the authors address the 3-Dimensional (3D) surface reconstruction of rockfall point clouds, here on denoted as rockfall objects, to assess the implications for volumetric analysis. Commonly used surface reconstruction algorithms are implemented and compared using a synthetic rockfall object and three natural rockfall events that occurred in the White Canyon, British Columbia, Canada.

1.1. Surface Reconstruction Background

Surface reconstruction aims to construct a complete surface model from a discrete point cloud input. Triangular meshes are often used as surface models because they are the simplest of polygon meshes, and therefore reduce the amount of computation used for their various computer graphical and 3D modeling applications. A key requirement for many applications, including accurate volume calculations, is that the constructed triangular mesh surface needs to be a watertight manifold surface with correct topology. A mesh is watertight if the mesh does not have any missing triangles (i.e., surface holes). A mesh is a manifold if: (1) each edge is incident to only one or two faces and; (2) the faces incident to a vertex form a closed or an open fan (Figure 1e,f).

A 2-manifold mesh must satisfy the following conditions: every point has a neighborhood homeomorphic to the open unit disk [15]. This neighbourhood can be continuously deformed to an open disk. Finally, every edge of a 2-manifold mesh is a manifold edge, with only two attached triangles (Figure 1b). If a mesh does not satisfy the aforementioned conditions (i.e. it is non-manifold), the mesh must be healed prior to computations. Healing the mesh involves hole-filling routines and manifold extraction.

Point cloud artifacts challenge the ability of surface reconstruction algorithms to produce watertight manifold meshes that approximate the surface topology of the object. Berger et al. [14] provide a review of the state of the art and challenges in surface reconstruction algorithms. In their review, they address common point cloud artifacts that impact on reconstruction algorithms. These artifacts include sampling density, noise, outliers, misalignment, and missing data (Figure 2).

1.2. Rockfall Objects

In this work, the authors make use of four rockfall objects as input for surface reconstruction. These objects include a synthetic rockfall and three natural rockfalls derived from sequential TLS scans in the White Canyon, British Columbia, Canada.

1.2.1. Houdini—Synthetic Rockfall Object

A synthetic rockfall object with a known volume was created to test the volumetric accuracy of the surface reconstruction methods. Houdini [16], a procedural modeling and animation software package, was used to create the rockfall object. The rockfall object was modeled from an initial ellipsoid with dimensions of 1 by 1 by 2 meters. A portion of the ellipsoid was intersected with a rectangular prism and then deleted to give the synthetic rockfall a hypothetical failure plane (Figure 3a–c). Surface topology was added to the ellipsoid and intersecting prism to give the synthetic object surface detail, which would challenge the surface reconstruction methods while representing realistic rockfall geometry. A notable detail is the lobe protruding on the right side of the object (Figure 3d).

The rockfall object has a volume of 4.68 m³. The volume was calculated by summing the signed volumes of tetrahedra with a common apex based at individual surface triangles. The surface mesh contains 642,573 faces, 632,010 vertices and is watertight manifold. The vertices were subsampled within the open-source software package CloudCompare [17] to minimum point spacings of 1, 2, 5 and 10 cm (Figure 3d).

1.2.2. Natural Rockfalls—White Canyon Rockfall Cases

In this study, the authors utilize three rockfall case studies that occurred in the White Canyon, British Columbia, Canada (50.266261°, −121.538943°). These events are a subset of the remotely sensed rockfall database from Bonneau et al. [13]. The database is derived from sequential TLS scans taken with an Optech Illris 3D-ER time-of-flight system. The three events correspond to events large enough to be of potential impact on the railway infrastructure and were interpreted to be the result of discrete individual events, based on fairly well constrained shape, relative to rockmass structure present at the rockfall source location (Figure 4). For additional details on the rockfall database and its construction, readers are referred to Bonneau et al. [13].

2. Methods

2.1. Surface Reconstruction Background

This section outlines the surface reconstruction approaches utilized within the study. The methods include Convex Hull, Alpha-shape and Power Crust.

2.1.1. Convex Hull

A Convex Hull defines the smallest possible polygon that contains every point within a set. More formally, the Convex Hull of a set of points S, in n-dimensions, is the intersection of all convex sets containing S [15]. Figure 5 displays an 2D example of a convex hull for a given set of points. In this work, the authors implement the Convex Hull algorithm in MATLAB [18].

2.1.2. Alpha-shapes

An Alpha-shape is a generalization of the convex hull of a point set [19]. For a given point set S, a family of Alpha-shapes can be defined. To define the Alpha-shape, an Alpha-complex must be defined, which is a subcomplex of the Delaunay triangulation approach. Therefore, for a given Alpha-radius, the Alpha-complex includes all the simplices in the Delaunay triangulation which have an empty circumscribing sphere with squared radius equal or smaller than the Alpha-radius. Here “empty” means that the open sphere does not include any points of S. The Alpha-shape is then simply the domain covered by the simplices of the Alpha-complex (Figure 6). Readers are referred to Edelsbrunner and Mücke [20] for more details on 3D Alpha-shape formulation. In this work, the authors implement the Alpha-shape algorithm in MATLAB [18].

An iterative Alpha-shape approach is developed in this study. The approach iterates through the list of possible Alpha-radii that produce unique shapes and finds the smallest Alpha-radius that produces a watertight manifold surface mesh. This process is accomplished by assessing if each edge in the mesh is attached with only two triangles, and therefore meets the definition of a 2-manifold mesh.

2.1.3. Power Crust

The Power Crust algorithm, originally developed by Amenta et al. [21] theoretically guarantees that for a given set of input points, the output surface is the watertight boundary of a three-dimensional polyhedral solid described by the approximate medial axis transform (MAT). A 2D illustrative example of the algorithm is displayed in Figure 7.

The algorithm relies on the concept that a given surface can be expressed as an infinite union of spheres centered on the inner medial axis by the MAT. The MAT is approximated by a subset of the Voronoi vertices of the input point cloud. These vertices are referred to as poles and are located near the medial axis. Polar balls are spheres surrounding the poles and touch the nearest input point samples. The radius of each polar ball determines the weighting applied to each pole. A power diagram of the weighted poles is used to determine the approximated inverse transform. The inverse transform divides the space into polyhedral cells. A subset of the polyhedral cells is labeled as the interior of the object. The Power Crust forms the surface that separates the inner cells from the outer cells. Further details on the Power Crust algorithm can be found in Amenta et al. [21]. In this work, the authors implemented the Power Crust algorithm in C++. The code is available from: https://github.com/ialhashim/powercrust.

As noted by Amenta et al. [21], not all Power Crust faces are necessarily triangles (Figure 8b). As a result, the authors implemented a triangulation scheme which checks each of the Power Crust faces and triangulates polygonal faces (Figure 8c).

2.2. Point to Mesh Comparisons

Point to mesh distance calculations were conducted between the input TLS point clouds of the rockfall objects and the triangulated surface meshes. All distance computations were completed in CloudCompare.

The distance from the point of interest to a triangle on a surface mesh in 3D, is defined as the point to the plane of the triangle if its projection onto the plane is inside the triangle. Otherwise, it is defined as the distance between the point and the closest edge of the triangle [22]. Figure 9 displays an example, where if P is the point of interest and ABC is a triangle included in the plane α, H is defined as the orthogonal projection of P onto the plane α. Therefore, there are two cases for the distance calculation:

• H is within the interior of the triangle ABC then the distance is:

  $$d\left(P, ABC\right)= \left|\right|\overrightarrow{HP}\left|\right|$$

  (1)
• H is within the interior of the triangle ABC then the distance is:

  $$d\left(P, ABC\right)= \left|\right|\overrightarrow{KP}\left|\right|$$

  (2)
  where, the point K is defined as the point on triangle ABC which minimizes the distance to H. If the triangle ABC has an associated normal vector, as an additional product of the calculation, the calculated distance can be signed (e.g., loss or gain). The normal vector of the triangle provides information if the point is within the interior or exterior of the surface represented locally by the triangle.

3. Results

3.1. Houdini—Synthetic Rockfall Object

Each of the sub-sampled point clouds were meshed using each of the surface reconstruction algorithms. The number of faces, vertices, and volume are displayed within

Table 1. Comparison of the surface reconstruction results for the four different subsampled point clouds of the synthetic rockfall object.

Point Spacing	Surface Reconstruction Method	Faces	Vertices	Volume (m3)	Volumetric Error (%)	Alpha-Radius (m)
1-cm	Default Alpha-shape	43,110	20,888	1.69	−63.90	0.3841
	Iterative Alpha-shape	21,046	10,525	5.18	10.65	0.62619
	Power Crust	954,648	464,855	4.68	−0.02	N/A
	Convex Hull	1182	593	5.87	25.39	Inf.
2-cm	Default Alpha-shape	16,782	8063	2.08	−55.63	0.4307
	Iterative Alpha-shape	10,500	5244	4.54	−2.93	0.56824
	Power Crust	338,322	166,774	4.68	−0.01	N/A
	Convex Hull	782	393	5.83	24.57	Inf.
5-cm	Default Alpha-shape	5940	2866	2.86	−38.82	0.493
	Iterative Alpha-shape	4392	2198	5.06	8.09	0.59193
	Power Crust	71,822	35,646	4.68	−0.09	N/A
	Convex Hull	514	259	5.74	22.69	Inf.
10-cm	Default Alpha-shape	2354	1157	3.62	−22.73	0.5377
	Iterative Alpha-shape	2070	1035	4.75	1.40	0.58245
	Power Crust	18,644	9261	4.66	−0.40	N/A
	Convex Hull	364	184	5.62	20.00	Inf.

. The volumetric error reflects the comparison of the volume from Houdini to the calculated volume of each of the triangulated surface meshes. The volume of each of the meshes was calculated using the divergence theorem [23].

Figure 10 displays example triangulated surface meshes of the synthetic rockfall object. The Convex Hull approach is shown to simplify and overestimate the surface geometry, while Power Crust is able to capture many of the details of the surface. The default Alpha-shape seems to capture a fair amount of detail, however produces boundary facets on the interior of the shape (Figure 10b), resulting in a drastic underestimation of volume. The iterative Alpha-shape method captures some of the surface details, while producing a hole-free shape.

3.2. Natural Rockfalls—White Canyon Rockfall Cases

Triangulated surface meshes were generated for each of the three rockfall case events using the methods described in Section 2. The face count, number of vertices, holes, and volume are shown for all three objects, using the four reconstruction methods, in

Table 2. Surface reconstruction results for Rockfall 1.

Surface Reconstruction Method	Faces	Vertices	Holes	Volume (m3)
Default Alpha-shape	6746	2937	30	0.25
Iterative Alpha-shape	2922	1464	0	1.34
Power Crust	37,412	18,597	0	1.45
Convex Hull	310	157	0	1.82

,

Table 3. Surface reconstruction results for Rockfall 2.

Surface Reconstruction Method	Faces	Vertices	Holes	Volume (m3)
Default Alpha-shape	2848	1337	11	0.94
Iterative Alpha-shape	1668	836	0	1.63
Power Crust	28,704	14,076	0	1.39
Convex Hull	222	113	0	2.44

and

Table 4. Surface reconstruction results for Rockfall 3.

Surface Reconstruction Method	Faces	Vertices	Holes	Volume (m3)
Default Alpha-shape	21,640	10,492	46	57.20
Iterative Alpha-shape	13,368	6688	0	129.63
Power Crust	183,044	90,359	0	123.81
Convex Hull	226	115	0	211.77

.

The default Alpha-shape in all three cases produced surface reconstructions with holes in the triangulated mesh. These holes result in smaller volumes being calculated for each of the rockfall objects. The Convex Hull, by definition, encompasses all the input points and therefore overestimates the volume of each of the rockfall objects. For all three cases, the Power Crust approach produced an order of magnitude more faces than the other reconstruction methods. It should be noted that although Power Crust interpolates the input point cloud, not all input points are vertices and not all vertices are input points. In comparison to the Alpha-shape or Convex Hull approaches, the input points form the vertices of the mesh.

A point to mesh distance computation was conducted for all three of the objects. The input point clouds were compared to the triangulated surface meshes generated using each of the four surface reconstruction implementations. The results of these distance computations for each of the three rockfall objects are depicted in Figure 11, Figure 12 and Figure 13. The statistics for the distances are shown in Figure 14, Figure 15 and Figure 16 and in

Table 5. Descriptive statistics for the distance computations for each surface reconstruction method for Rockfall 3.

Surface Reconstruction Method	Mean (m)	Standard Deviation (m)	Variance (m)	Skewness
Default Alpha-shape	0.0017	0.0055	3.02 * 10−5	4.22
Iterative Alpha-shape	0.011	0.015	2.12 * 10−4	2.36
Power Crust	4.8 * 10−5	0.0047	2.22 * 10−5	0.10
Convex Hull	0.038	0.029	8.62 * 10−4	1.39

,

Table 6. Descriptive statistics for the distance computations for each surface reconstruction method for Rockfall 2.

Surface Reconstruction Method	Mean (m)	Standard Deviation (m)	Variance (m)	Skewness
Default Alpha-shape	0.013	0.023	5.3 * 10−4	2.32
Iterative Alpha-shape	0.025	0.036	0.0013	2.12
Power Crust	−5.1 * 10−5	0.0060	3.6 * 10−5	−2.35
Convex Hull	0.078	0.072	0.0052	1.17

and

Table 7. Descriptive statistics for the distance computations for each surface reconstruction method for Rockfall 3.

Surface Reconstruction Method	Mean (m)	Standard Deviation (m)	Variance (m)	Skewness
Default Alpha-shape	0.024	0.048	0.0023	2.97
Iterative Alpha-shape	0.051	0.088	0.0078	2.74
Power Crust	0.0062	0.040	0.0016	6.0056
Convex Hull	0.35	0.19	0.036	0.52

.

For all three cases, each of the surface reconstruction approaches over- or under-interpolates the input point cloud. The default Alpha-shape approach not only interpolates in comparison to the input point cloud, but also produces holes in the surface mesh. The holes corresponded to locations of lower point density in comparison to the overall surface. The Power Crust and iterative Alpha-shape approaches produced watertight meshes. Interestingly, the Power Crust and iterative Alpha-shape approach produced similar volumes that differed between −17% to 8% of the Power Crust volume. However, the iterative Alpha-shape approach had the second highest degree of interpolation of the input point clouds in the reconstruction. The highest degree of interpolation was the with the Convex Hull approach, which can be seen in the results of the distance computations (Figure 11, Figure 12 and Figure 13). In all three cases, Power Crust was the only algorithm that displayed under-interpolation (i.e., negative change) when compared to the input point cloud in comparison to the other reconstruction methods.

Figure 14, Figure 15 and Figure 16 display histograms of the distance computations displayed in Figure 11, Figure 12 and Figure 13 for each rockfall object and surface reconstruction method. Table 4, Table 5 and Table 6 display the mean, standard deviation, variance and skewness of the distance computation for each of the surface reconstruction methods.

4. Discussion and Conclusions

Surface reconstruction remains a challenge for objects in natural environments. At present, there is a trade-off between accurate surface topology reconstruction and ensuring the mesh is watertight manifold. In natural environments point cloud artifacts will be present in many situations, however, optimal survey design and post-processing techniques can help reduce these artifacts.

Previous studies that make use of Alpha-shape approach, such as van Veen et al. [11], Carrea et al. [9] and others, do not address nor specify how they selected an appropriate Alpha-radius for surface reconstruction. Selecting a constant Alpha-radius for a study slope is not appropriate in most situations without significant manual post-processing. Manual post-processing is subjective and not practical at highly active sites where many rockfall objects must be analyzed. As TLS acquisition and processing workflows are becoming automated, practitioners require methods that are reliable and ensure that the volumes being calculated are as accurate as possible. As demonstrated with the rockfall objects analyzed in this study, picking too small of an Alpha-radius results in holes in the mesh. Holes in the surface mesh result in underestimating the volume of the rockfall being analyzed. To overcome the aforementioned issues, an optimized iterative Alpha-shape approach has been presented that ensures watertight, manifold triangulated surface meshes that can be used to calculate the volume of rockfall events. However, as demonstrated with all the rockfall objects in this study, the iterative Alpha-shape approach interpolates the reconstruction. The volumetric error, for the synthetic shape, as an example, ranged from −2.9% to 11%. However, the approach is robust, automated and will ensure the reconstruction is watertight.

Power Crust was shown to be a robust approach that eliminates the need for the hole-filling, or post-processing steps for manifold extraction. In addition, in comparison to Alpha-shape approaches, a unique Alpha-radius does not need to be defined. However, with Power Crust, the results of reconstruction produce an order of magnitude more faces in comparison to other methods. These faces and vertices could create microscale topology which is not present in the input point cloud. In addition, there is no way to verify if these subtle topological artifacts exist since the algorithm is interpolating between the sample inputs. Therefore, if the back scar of a rockfall is being reconstructed, there may be topological artifacts induced as a result of the reconstruction. However, mesh simplification schemes exist as mentioned in Amenta et al. [21], to reduce the number of faces while still preserving the detail of the input point cloud.

All the surface reconstruction methods presented in the study will operate on any point cloud input. The point clouds used in this study were based on TLS scans, however, Structure-from-Motion Multi-View-Stereo photogrammetry (SfM-MVS) and aerial laser scanning (ALS) point clouds could also be used. Each data acquisition method (e.g., TLS, ALS) has its respective challenges with generating the most complete representation of the slope surface [24,25]. Data acquisition aside, the output of the surface reconstruction is a direct function of the point cloud input. If there are significant point cloud artifacts, such as occlusions in the input point cloud, each of the surface reconstruction algorithms will require significant interpolation to fill in the areas where there are no points present. Further work is required to test the application of other surface reconstruction algorithms (e.g., Berger et al. [14]) and their applicability to generating triangulated surface meshes of rockfall objects.

In conclusion, surface reconstruction aims to build a triangulated surface mesh representation based on an input point cloud. Point clouds may contain artifacts which present challenges for surface reconstruction algorithms. In this work, four surface reconstruction algorithms were used to assess the volumetric implications of surface reconstruction for rockfall back-analysis. From this analysis the authors can draw the following general conclusions:

• The Convex Hull approach in all cases over-interpolates the input point cloud. The method will generate the minimum smallest polygon that encapsulates all points. Therefore, as demonstrated in all rockfall cases within this study, concave features on all objects are not captured. The volume, using the Convex Hull method, in almost all cases will be an over-estimate.
• The Alpha-shape approach is sensitive to point density, especially when it is non-uniform. As a result, it is difficult and sometimes impossible to choose an Alpha-radius to balance hole-filling against loss of detail.
• The iterative Alpha-shape approach tries to find a balance between hole-filling and loss of detail. It is a robust approach that can be integrated into automated TLS processing workflows. In all cases, analyzed in this study the iterative Alpha-shape over-interpolated the surface in comparison to the input point cloud.
• Power Crust appears to the optimal surface reconstruction algorithm that was tested in the present study. The approach theoretically guarantees that the output surface will be watertight and manifold. The staggering amount of faces in the output surface mesh may present problems for some applications however simplification approaches exist.
• Each approach has its drawbacks, however, it is critical to ensure that the mesh is watertight and manifold for accurate volumetric calculations. The volume of Rockfall 3 in this study rose from approximately 57 m³ (e.g., default Alpha-shape) to around 125 m³ (e.g., iterative Alpha-shape & Power Crust) once the object was watertight (e.g., holes were filled). That represents a gain in 68 m³ of volume attributed to the rockfall event. This case is an excellent example to show the implications of method selection on frequency-magnitude analyses and using sequential TLS datasets to help design mitigation measures. As more monitoring approaches become automated, the importance of having quality control on the inputs to databases becomes paramount.