MATLAB Virtual Toolbox for Retrospective Rockfall Source Detection and Volume Estimation Using 3D Point Clouds: A Case Study of a Subalpine Molasse Cliff

: The use of 3D point clouds to improve the understanding of natural phenomena is currently applied in natural hazard investigations, including the quantiﬁcation of rockfall activity. However, 3D point cloud treatment is typically accomplished using nondedicated (and not optimal) software. To ﬁll this gap, we present an open-source, speciﬁc rockfall package in an object-oriented toolbox developed in the MATLAB ® environment. The proposed package offers a complete and semiautomatic 3D solution that spans from extraction to identiﬁcation and volume estimations of rockfall sources us-ing state-of-the-art methods and newly implemented algorithms. To illustrate the capabilities of this package, we acquired a series of high-quality point clouds in a pilot study area referred to as the La Cornalle cliff (West Switzerland), obtained robust volume estimations at different volumetric scales, and derived rockfall magnitude–frequency distributions, which assisted in the assessment of rockfall activity and long-term erosion rates. An outcome of the case study shows the inﬂuence of the volume computation on the magnitude–frequency distribution and ensuing erosion process interpretation.


Introduction
3D point clouds are commonly used in different scientific domains and in geosciences. In geosciences, the use of point clouds mainly focuses on capturing outcrop geometry [1], mapping geological layers [2], or investigating surface morphological changes [3,4]. In the topic of natural hazards, point clouds are frequently used for mapping landslides and monitoring topographic changes resulting from rockfall activity or mass movement [5].
In natural hazards, point clouds originate from active sensors, such as laser scanners, also known as LiDAR technologies [6,7], or from passive sensors and workflows, such as classic photogrammetry or structure-from-motion (SFM) [8][9][10][11]. During the last decade, the popularization of point clouds acquired from ground-based LiDAR has allowed the rapid mass collection of slope surface measurements [12,13]. These high-accuracy and highspatial-resolution data have opened new ways of investigating natural hazards [14][15][16][17]. Working with multitemporal high-density point clouds (typically over 100 million points per scan) or large datasets requires a widely applicable and powerful programming language. This has led to the development of open-source software dedicated to processing point clouds, most of which were originally dedicated to civil engineering or robotics, such as CloudCompare [18] or Point Cloud Library [19]. However, several challenges still exist in the management and extraction of this large amount of data collected for natural hazard studies, in particular for detecting and estimating the volume of rockfalls in their source areas. In addition, a key issue is how to gather all tools into a single structure, as has been done for aerial laser scanning in the domains of forestry and land use [20].

1.
Points belonging to topographic changes assumed to result from rockfalls.

2.
Points belonging to unchanged topography assumed to be stable surfaces.
The changes between two point cloud acquisitions are identified by a classical method described in [40]. The detection of rockfall sources uses the difference in distance between two point clouds of different epochs. The computation of the difference is performed using a point-to-surface comparison to obtain the shortest distance (Euclidean distance) from each point to the surface as d i = ∆(P i , S), where (P i ) is, for example, a (i) point in the pre-rockfall event (pre) point cloud (P) of size (n), (S) is the surface built from the post-rockfall event (post) point cloud using the triangulation mesh, and (d i ) is the computed signed distance along the local normal of (S). If no acquisition bias exists in the point cloud, the distribution of distance differences without surface change follows a normal distribution centered on zero. In the locations at which a change in topography occurs, the distance comparison (d i ) must be larger than the standard deviation (σ). According to [5], the standard deviation of the measurements between two epochs can be high and depends on multiple factors, such as the quality of the point cloud datasets, the density of the points, the presence of vegetation, the roughness of the relief, and the quality of the alignment between the point clouds and/or the acquisition locations between the two epochs (LiDAR position or picture positions). According to [41], point cloud points are assumed to be indicative of topographic changes (i.e., here, a rockfall source) without ambiguity when the pointto-surface comparison distances are larger than two times the standard deviation (2σ). As proposed by [41], this threshold can be improved by applying a spatial filter using the mean point-to-surface comparison distance of a point and its 25 nearest neighbors. Thus, the thresholding conditions are defined as follows: Ppost i ⊂ rock f alls source when The points representing the conditions in Equations (1) and (2) are merged for each neighborhood as one set of points belonging to a rockfall source. The result is spatially distributed points that form clusters corresponding to rockfall sources and isolated points associated with noise ( Figure 1).
The changes between two point cloud acquisitions are identified by a classical method described in [40]. The detection of rockfall sources uses the difference in distance between two point clouds of different epochs. The computation of the difference is performed using a point-to-surface comparison to obtain the shortest distance (Euclidean distance) from each point to the surface as di = Δ(Pi, S), where (Pi) is, for example, a (i) point in the pre-rockfall event (pre) point cloud (P) of size (n), (S) is the surface built from the post-rockfall event (post) point cloud using the triangulation mesh, and (di) is the computed signed distance along the local normal of (S). If no acquisition bias exists in the point cloud, the distribution of distance differences without surface change follows a normal distribution centered on zero. In the locations at which a change in topography occurs, the distance comparison (di) must be larger than the standard deviation (σ). According to [5], the standard deviation of the measurements between two epochs can be high and depends on multiple factors, such as the quality of the point cloud datasets, the density of the points, the presence of vegetation, the roughness of the relief, and the quality of the alignment between the point clouds and/or the acquisition locations between the two epochs (LiDAR position or picture positions). According to [41], point cloud points are assumed to be indicative of topographic changes (i.e., here, a rockfall source) without ambiguity when the point-to-surface comparison distances are larger than two times the standard deviation (2σ). As proposed by [41], this threshold can be improved by applying a spatial filter using the mean point-to-surface comparison distance of a point and its 25 nearest neighbors. Thus, the thresholding conditions are defined as follows: The points representing the conditions in Equations (1) and (2) are merged for each neighborhood as one set of points belonging to a rockfall source. The result is spatially distributed points that form clusters corresponding to rockfall sources and isolated points associated with noise ( Figure 1). Figure 1. Sketch of the process of localizing rockfall sources. Where P are the points, S is the reconstructed surface from point cloud at different epochs, and dP is the difference in distance from the point-to-surface comparison. The subscripts pre and post denote pre-and post-rockfall event acquisition, respectively. Sketch of the process of localizing rockfall sources. Where P are the points, S is the reconstructed surface from point cloud at different epochs, and d P is the difference in distance from the point-to-surface comparison. The subscripts pre and post denote pre-and post-rockfall event acquisition, respectively.

Step 2: Clustering Rockfall Sources
The second step consists of the identification of each cluster representing a single source and the removal of outliers. Currently, the segmentation of point clouds into clusters uses the density-based clustering algorithm DBSCAN [42], as applied in [24]. The principle behind DBSCAN [42] consists of scanning a point cloud with a selected neighboring radius (ε) assuming the availability of the minimum number of points (k) required to create a cluster. The choice of the minimum number of points (k) and the neighboring radius (ε) depends on the point cloud point density and the volume to be detected. This point is not intuitive and is a limitation of DBSCAN as it is necessary to first identify a reasonable measure of similarity for the dataset, before selecting the optimal (ε) It can be determined by a trial-and-error procedure.
A solution to these shortcomings is to use the OPTICS algorithm (ordering points to identify the clustering structure) [43] as a variant of DBSCAN. When using the OPTICS algorithm, only the minimum number of points (k) considered as a cluster is needed. Then, the value of (ε) can be calculated from the chosen (k). Based on the reachability plot structure, an optimal neighborhood radius (ε) containing a predefined (k) points is given by [43] where (n) is the number of points in the dataset, (m) is the dimensionality of the experimental space, (Γ) is the gamma function, and (V) is the volume of the space formed by m points. The OPTICS algorithm sorts the data within a point cloud according to their distance and core distance and categorizes the points into one of three categories: core points, border points, or outliers ( Figure 2). A point (i) of a point cloud is defined as follows: • A core point, if the neighborhood of radius (ε), has at least k-points (reachable points); • A border point possesses at least one core point within a radius (ε); • An outlier is a point with no point or no core point within its radius (ε).

Step 2: Clustering Rockfall Sources
The second step consists of the identification of each cluster representing a single source and the removal of outliers. Currently, the segmentation of point clouds into clusters uses the density-based clustering algorithm DBSCAN [42], as applied in [24]. The principle behind DBSCAN [42] consists of scanning a point cloud with a selected neighboring radius (ε) assuming the availability of the minimum number of points (k) required to create a cluster. The choice of the minimum number of points (k) and the neighboring radius (ε) depends on the point cloud point density and the volume to be detected. This point is not intuitive and is a limitation of DBSCAN as it is necessary to first identify a reasonable measure of similarity for the dataset, before selecting the optimal (ε) It can be determined by a trial-and-error procedure.
A solution to these shortcomings is to use the OPTICS algorithm (ordering points to identify the clustering structure) [43] as a variant of DBSCAN. When using the OPTICS algorithm, only the minimum number of points (k) considered as a cluster is needed. Then, the value of (ε) can be calculated from the chosen (k). Based on the reachability plot structure, an optimal neighborhood radius (ε) containing a predefined (k) points is given by [43] where (n) is the number of points in the dataset, (m) is the dimensionality of the experimental space, (Γ) is the gamma function, and (V) is the volume of the space formed by m points. The OPTICS algorithm sorts the data within a point cloud according to their distance and core distance and categorizes the points into one of three categories: core points, border points, or outliers ( Figure 2). A point (i) of a point cloud is defined as follows: • A core point, if the neighborhood of radius (ε), has at least k-points (reachable points); • A border point possesses at least one core point within a radius (ε); • An outlier is a point with no point or no core point within its radius (ε). Now, if a point (i) is a core point, then it forms a cluster together with all points (core or border) that are reachable from it. They are mutually densely connected. Each cluster contains at least one core point. Now, if a point (i) is a core point, then it forms a cluster together with all points (core or border) that are reachable from it. They are mutually densely connected. Each cluster contains at least one core point.
According to [43], the OPTICS algorithm is a well-suited clustering method for datasets with a large number of points. In addition, OPTICS was successfully used in chemometrics to reveal clusters of arbitrary shapes with differing densities [44].
After processing, the point cloud is divided into a series of point cloud sets corresponding to its clusters and outliers. Then, at the end of step 2, a cluster is considered a rockfall source when it is a compound of points from different epochs, pre-and post-rockfall sources. The single epoch clusters and outliers are removed.

Step 3: Rockfall Source Volume Estimation
The last step consists of calculating the volumes of the different sets of clustered points identified as rockfall sources. The volume computation utilizes the full 3D point cloud using the α-shape hull algorithm [45,46] to take into account the natural geometric complexities encountered for rockfall source shapes. This algorithm, originally developed in the domain of computer graphics to recreate complex surfaces, is applied to geometrical studies in biosciences, [47], and its uses include visualization and volume estimations. The α shape is a mathematical expression used for the representation of complex shapes (convex or concave) by the linear approximation of the original shape from a set of points [48]. According to [45], the α shape is the generalized form of the concave hull in which it is a complex of Delaunay triangulation (DT) in 2D or 3D. For a given value of α, the α-complex includes all the simplexes in the DT that have an empty circumscribing sphere with a squared radius equal to or smaller than α ( Figure 3).
According to [43], the OPTICS algorithm is a well-suited clustering method for datasets with a large number of points. In addition, OPTICS was successfully used in chemometrics to reveal clusters of arbitrary shapes with differing densities [44].
After processing, the point cloud is divided into a series of point cloud sets corresponding to its clusters and outliers. Then, at the end of step 2, a cluster is considered a rockfall source when it is a compound of points from different epochs, pre-and postrockfall sources. The single epoch clusters and outliers are removed.

Step 3: Rockfall Source Volume Estimation
The last step consists of calculating the volumes of the different sets of clustered points identified as rockfall sources. The volume computation utilizes the full 3D point cloud using the α-shape hull algorithm [45,46] to take into account the natural geometric complexities encountered for rockfall source shapes. This algorithm, originally developed in the domain of computer graphics to recreate complex surfaces, is applied to geometrical studies in biosciences, [47], and its uses include visualization and volume estimations. The α shape is a mathematical expression used for the representation of complex shapes (convex or concave) by the linear approximation of the original shape from a set of points [48]. According to [45], the α shape is the generalized form of the concave hull in which it is a complex of Delaunay triangulation (DT) in 2D or 3D. For a given value of α, the α-complex includes all the simplexes in the DT that have an empty circumscribing sphere with a squared radius equal to or smaller than α ( Figure 3). To define the closeness to a real shape or the envelope (Sα) formed by a set of points, the value used is an α value corresponding to a research radius in the point cloud ranging from 0 to ∞ and follows the subsequent conditions: If 0 < α < ∞, Sα will be the largest polyhedron or shape connecting m points of the point cloud.
In addition, we look for an efficient determination of the volume of this specific shape (Sα). As Sα is a complex of DT, we use the DT to decompose volumes defined by the envelope of S into tetrahedrons. The volume of a i tetrahedron with a triangular base of a given area, A, and a height, h, are given as follows: The volume of the source of a rockfall is equal to the sum of all n tetrahedron volumes, Vi, inside an α-shaped Sα defined as follows: To define the closeness to a real shape or the envelope (S α ) formed by a set of points, the value used is an α value corresponding to a research radius in the point cloud ranging from 0 to ∞ and follows the subsequent conditions: If 0 < α < ∞, S α will be the largest polyhedron or shape connecting m points of the point cloud.
In addition, we look for an efficient determination of the volume of this specific shape (S α ). As S α is a complex of DT, we use the DT to decompose volumes defined by the envelope of S α into tetrahedrons. The volume of a i tetrahedron with a triangular base of a given area, A, and a height, h, are given as follows: The volume of the source of a rockfall is equal to the sum of all n tetrahedron volumes, V i , inside an α-shaped S α defined as follows: This process is applied to all rockfall sources. The precision of the volume calculation depends on the α value. It is possible to have an intuitive perception of the real volume by plotting the volume versus the α value, and the best α value is then located at the first asymptotic behavior instance ( Figure 4).
This process is applied to all rockfall sources. The precision of the volume calculation depends on the α value. It is possible to have an intuitive perception of the real volume by plotting the volume versus the α value, and the best α value is then located at the first asymptotic behavior instance ( Figure 4).

Case Study
To illustrate the application of the abovementioned methodology and its implementation in MATLAB, we carried out a series of point cloud acquisitions with a terrestrial LiDAR device at a study site located a few kilometers eastward of Lausanne (Vaud, Switzerland) in the Lavaux region ( Figure 5A) within the Molasse Swiss Plateau [49]. The region is affected by several critical slow-moving landslides and numerous rockfall activities [50].

Case Study
To illustrate the application of the abovementioned methodology and its implementation in MATLAB, we carried out a series of point cloud acquisitions with a terrestrial LiDAR device at a study site located a few kilometers eastward of Lausanne (Vaud, Switzerland) in the Lavaux region ( Figure 5A) within the Molasse Swiss Plateau [49]. The region is affected by several critical slow-moving landslides and numerous rockfall activities [50].
The pilot study area, referred to as the La Cornalle cliff, is part of the lateral scarp of a slow-moving landslide ( Figure 5B) that has been well documented since the 18th century and has been active for 10,000 years, resulting in several hazardous events [50]. The cliff, which is over 35 m high and 110 m long ( Figure 5C), is an interesting actual-scale laboratory used to study rockfall processes, including the triggering conditions needed for the erosion of molassic rock, which leads to estimation of the erosion rate induced by rockfall activity. The cliff is hardly accessible for direct measurement.
The La Cornalle cliff is composed of subhorizontal alternations of decimetric to metric beds (0.1-3 m) of sandstone and marls ( Figure 5D) of Chattian age from the Lower Freshwater Molasse [51]. Molassic rocks are generally weak rocks that are highly affected by erosional processes. The bedding is monocline and dips in the opposite direction of the slope. Four sets of discontinuities are present: stratification S 0 125/14 • , joint system J 1 234/86 • subparallel to the cliff face, and joint systems J 2 150/75 • and J 3 325/80 • nearly perpendicular to the slope. Discontinuities are clearly present on sandstone beds. The spacing varies from 0.2 to 1.0 m for J 1 and from 0.15 to 2 m for J 2 and J 3 , respectively. The pilot study area, referred to as the La Cornalle cliff, is part of the lateral scarp of a slow-moving landslide ( Figure 5B) that has been well documented since the 18th century and has been active for 10,000 years, resulting in several hazardous events [50]. The cliff, which is over 35 m high and 110 m long ( Figure 5C), is an interesting actual-scale laboratory used to study rockfall processes, including the triggering conditions needed for the erosion of molassic rock, which leads to estimation of the erosion rate induced by rockfall activity. The cliff is hardly accessible for direct measurement.
The La Cornalle cliff is composed of subhorizontal alternations of decimetric to metric beds (0.1-3 m) of sandstone and marls ( Figure 5D) of Chattian age from the Lower Freshwater Molasse [51]. Molassic rocks are generally weak rocks that are highly affected by erosional processes. The bedding is monocline and dips in the opposite direction of the slope. Four sets of discontinuities are present: stratification S0 125/14°, joint system J1 234/86° subparallel to the cliff face, and joint systems J2 150/75° and J3 325/80° nearly perpendicular to the slope. Discontinuities are clearly present on sandstone beds. The spacing varies from 0.2 to 1.0 m for J1 and from 0.15 to 2 m for J2 and J3, respectively.
The significant rockfall sources located in the sandstone beds are caused by the differentiated erosion rate between marls and sandstone. This phenomenon induces the overhanging and toppling failure mechanisms of sandstone beds. Numerous rockfalls accumulate material on the top of the slow-moving landslide [52], which promotes movement. Based on structural analysis, the expected block volumes could vary from 0.003 m 3 (0.1 × 0.2 × 0.15 m) for the minimal volume to 6.0 m 3 (3.0 × 1.0 × 2.0 m) for the maximal volume.
We acquired a first point cloud of the cliff from a ground-based LiDAR survey in June 2010. A second acquisition was completed in September 2012; afterward, the cliff was The significant rockfall sources located in the sandstone beds are caused by the differentiated erosion rate between marls and sandstone. This phenomenon induces the overhanging and toppling failure mechanisms of sandstone beds. Numerous rockfalls accumulate material on the top of the slow-moving landslide [52], which promotes movement. Based on structural analysis, the expected block volumes could vary from 0.003 m 3 (0.1 × 0.2 × 0.15 m) for the minimal volume to 6.0 m 3 (3.0 × 1.0 × 2.0 m) for the maximal volume.
We acquired a first point cloud of the cliff from a ground-based LiDAR survey in June 2010. A second acquisition was completed in September 2012; afterward, the cliff was monitored by LiDAR four times per year (seasonally) up to May 2015. The data were collected with an Optech/Teledyne ILIRS 3D ER with theoretical accuracy of 7 mm at 100 m and a standard deviation of 10 mm [53]. The footprint at the cliff range (~300 m) was approximately 50 mm. The effective point surface density was of 860 pts/m 2 .

Results
The results are illustrated to show some output from typical datasets from 3D point clouds. Table 1 reports the input parameters used in the illustrated results, which were obtained automatically by the algorithm and manually by the user. When available, the automatic parameters proposed by the software were applied. The user-defined parameters were as follows: With a minimum of 34 points, it is expected to detect volumes larger than 0.5 L according to the described point surface density. The α value was determined for each point rockfall source based on the procedure described in Section 3.3. As an example, Figure 6A shows the observed and identified major rockfall events from sandstone beds between September 2012 and March 2013 during which the northern part of the cliff experienced important rockfall activity. Thirty-one rockfall events were identified via photographic comparison, mostly the important volumes. Figure 6B shows the identified rockfall sources after computation using the developed toolbox and applying the presented methods on the LiDAR data collected between September 2012 and March 2013. Each rockfall source belongs to a different-colored point cloud. In total, 49 individualized rockfall sources extracted with the presented procedure were identified. When comparing the two results, we noted that the difference in the identified number of rockfall events comes from the small undetected rockfall sources in the picture analysis. We analyzed the sensitivity of the computed volume for a restrained series of identified rockfall source varying the α values for the whole series. We increased the α values by increment up to infinity, hence increasing the trend toward a convex shape. Figure 7 shows the magnitude-frequency relationship of the series of computed volumes from the identified rockfall sources with different complex shapes, between Fall 2012 and Spring 2013, illustrated in Figure 6B. Figure 7 outlines the influence of the α value on each series of computed volumes from rockfall sources using power-law regression on each magnitude-frequency distribution. The magnitude-frequency representation is used here to To illustrate the output results of the presented procedure, and particularly the computed rockfall volumes, we used the classic rockfall magnitude-frequency representa-tion [22,31,[33][34][35][36][37]. According to [37], volumes from rockfall source areas are assumed to follow a power-law regression: where V is the rockfall volume, n is the number of rockfalls larger than V occurring in a rock wall during an investigation period, and the constant, a, represents the number of rockfalls whose volume is greater than 1 m 3 . It depends on the size of the cliff, the length of the investigation period, and the geological and geomorphological context. The exponent, b, depends on the geological and geomorphological context only. Thereby, recent studies showed that a and b are correlated to the GSI (Geological Strength Index) of the rock cliff [37]. In addition, the value of b indicates the proportion of small events as compared to larger events. Therefore, this is important in the context of a power-law distribution, where small events in sum could contribute significantly to overall volume loss. The powerlaw regression has some limits as described in [22], as we have potential temporal and/or spatial resolution bias. We analyzed the sensitivity of the computed volume for a restrained series of identified rockfall source varying the α values for the whole series. We increased the α values by increment up to infinity, hence increasing the trend toward a convex shape. Figure 7 shows the magnitude-frequency relationship of the series of computed volumes from the identified rockfall sources with different complex shapes, between Fall 2012 and Spring 2013, illustrated in Figure 6B. Figure 7 outlines the influence of the α value on each series of computed volumes from rockfall sources using power-law regression on each magnitudefrequency distribution. The magnitude-frequency representation is used here to show the trend. The regressions are made from volumes ≥0.01 m 3 , which is the expected smallest volume identified without ambiguity from the point cloud. We observed an increase in the number of volumes greater than 1 m 3 , or an increasing parameter a for an increase of the α value. In contrast, we observed a decrease in the b exponent. For values α < 0.1, which are a too low research radius, we observed that most of the reconstructed volumes are not filled or contained inner holes, as in Figure 4B. For a value 0.1 ≤ α ≤ 1.25, we can reconstruct volumes conserving a close-to-reality geometry depending on block shape complexity. For a value α > 1.25, we start connecting farther points, increasing shape convexity for all rockfall volumes, and moving away from real geometry up to an infinite α value, which indicates the maximal convex shape and the maximal rockfall volume (i.e., Figure 4D).
We applied the presented methodology to the overall monitoring period from 2010 to 2015 data. Thus, we identified 394 rockfall sources. Based on observation made on Figure 7, we applied on each of these rockfall sources a detailed assessment of the volume in order to compute close-to-reality volume, applying the method illustrated in Section 3.3 and Figure 4. As final result, we obtained a total volume of 105.2 m 3 , ranging from 0.0015 to 7.63 m 3 with a mean volume of 0.2 m 3 (see Figure 8). Less than 1% of the rockfall sources were smaller than 0.003 m 3 (0.1 × 0.2 × 0.15 m), which is the smallest volume expected from structural analysis. Less than 25% of the rockfall sources are smaller than 0.01 m 3 (0.2 × 0.2 × 0.2 m), which is the expected smallest volume identified without ambiguity from the point cloud. In Figure 8, the regression shown in red represents the rockfall sources with volumes higher than 0.01 m 3 , and the other regression shown in blue represents all volumes, including the smallest volumes. Green triangles were not included in the regression analysis because they represent multiple rockfall sources and were too scarce during the period of monitoring; thus, they were not considered to be representative [22,54]. We observed that the choice of the volume interval considered in the regression affects the a and b parameters. With distribution containing all rockfall sources, a increases and b decreases compare to distribution with only rockfall sources >0.01 m 3 . We applied the presented methodology to the overall monitoring period from 2010 to 2015 data. Thus, we identified 394 rockfall sources. Based on observation made on Figure 7, we applied on each of these rockfall sources a detailed assessment of the volume in order to compute close-to-reality volume, applying the method illustrated in Section 3.3 and Figure 4. As final result, we obtained a total volume of 105.2 m 3 , ranging from 0.0015 to 7.63 m 3 with a mean volume of 0.2 m 3 (see Figure 8). Less than 1% of the rockfall sources were smaller than 0.003 m 3 (0.1 × 0.2 × 0.15 m), which is the smallest volume expected from structural analysis. Less than 25% of the rockfall sources are smaller than 0.01 m 3 (0.2 × 0.2 × 0.2 m), which is the expected smallest volume identified without ambiguity from the point cloud. In Figure 8, the regression shown in red represents the rockfall sources with volumes higher than 0.01 m 3 , and the other regression shown in blue represents all volumes, including the smallest volumes. Green triangles were not included in the regression analysis because they represent multiple rockfall sources and were too scarce during the period of monitoring; thus, they were not considered to be representative [22,54]. We observed that the choice of the volume interval considered in the regression affects the a and b parameters. With distribution containing all rockfall sources, a increases and b decreases compare to distribution with only rockfall sources >0.01 m 3 .  using the optimal α value according to method described in Section 3.3 for each rockfall source. One power law is fitted on overall rockfall sources and a second on the rockfall sources with volume >0.01 m 3 . The green triangles were clearly identified as multiple rockfall sources and are not considered for the regression fit.

Discussion
From a processing and methodological standpoint, we note that applying the presented package for the detection and volume calculation of rockfall sources based on point clouds provides the benefit of being less time intensive compared to procedures that use manual methods and multiple software programs. This advantage is even more relevant when monitoring is performed on large outcrops based on several time series. The presented approach provides automatic thresholding, making it a robust method. Based on our experience, the presented results illustrate that segmentation permits efficient extraction of rockfall sources. However, some limitations appear in the point-to-surface method in terms of precision. Improvements could be made using finer point-to-surface comparisons in the first step by, for example, using the M3C2 algorithm comparison developed Figure 8. Magnitude-frequency plot of the rockfall sources from the full survey period (2010-2015) using the optimal α value according to method described in Section 3.3 for each rockfall source. One power law is fitted on overall rockfall sources and a second on the rockfall sources with volume >0.01 m 3 . The green triangles were clearly identified as multiple rockfall sources and are not considered for the regression fit.

Discussion
From a processing and methodological standpoint, we note that applying the presented package for the detection and volume calculation of rockfall sources based on point clouds provides the benefit of being less time intensive compared to procedures that use manual methods and multiple software programs. This advantage is even more relevant when monitoring is performed on large outcrops based on several time series. The presented approach provides automatic thresholding, making it a robust method. Based on our experience, the presented results illustrate that segmentation permits efficient extraction of rockfall sources. However, some limitations appear in the point-to-surface method in terms of precision. Improvements could be made using finer point-to-surface comparisons in the first step by, for example, using the M3C2 algorithm comparison developed by [55] and investigated in [56]. The DBSCAN algorithm used in clustering datasets requires a homogenously spaced point cloud since the minimum number of points and ε cannot then be chosen appropriately to detect all clusters at once. If the data and scale are not well analyzed before applying this method, the results can be biased. This limitation is eclipsed by the benefits of the OPTICS algorithm, as it offers an alternative to requiring an optimal ε, even if the user does not know the ε value. The advantages of using the OPTICS algorithm are as follows: (a) This algorithm does not require prior knowledge of the number of clusters in the data; (b) this variation of the DBSCAN algorithm enables the identification arbitrarily shaped clusters and clusters completely surrounded by other clusters; and (c) using OPTICS classification, we are able to remove outliers. However, the proposed algorithm also has some disadvantages, the OPTICS algorithm was not entirely deterministic; depending on the order in which the data are processed, border points were reachable from more than one cluster. This situation could occur if two or more rockfall sources were very close.
For the volume computation, it is possible to say that using the α shape takes full advantage of the 3D surface (with overhanging, etc.) and does not reduce it to 2.5D [30]. We see that for shapes that were more complex, our estimates approach the real volume and avoid overestimation (see Figure 4D). The α shape can compute a very complex shape, as often observed in natural rockfall geometry related to a discontinuous layout. The main problem with α shapes is the determination of the best α values for surface reconstruction when the point cloud does not present a uniform density of points [32]. Some studies have shown improvement of the surface and volume computation when using the α shape as was done in [57]. Furthermore, some issues with holes from shadowing during acquisition were encountered, as was also observed by [25] with the potential to negatively impact the accuracy of the results. This indicates that acquiring the point cloud from multiple positions can help overcome these errors.
Concerning the results from the case study, a good correlation can be seen between the rockfall sources identified during fieldwork (green in Figure 6A) and sources detected with our package (shown in Figure 6B). The segmentation of rockfall sources worked sufficiently for the efficient identification of volumes at different scales, e.g., from 0.0025 to 10 m 3 with a point spacing at cliff range of 3.4 cm.
When looking at the volume-frequency plot in Figure 7, we observed a decrease in b exponents, when increasing α value, which is in agreement with previous studies pointing out that the decrease in the value of b indicates a rise in the portion of larger events as compared to smaller events [33]. This is logically correlated to increasing α value, which generates a more convex and hence larger volume. Therefore, we observed variations in the volume ranging from 5% to 90% when computing concave or convex shapes depending on the selected α value and the complex geometry of the rockfall sources. The variation is particularly significant for large volumes with an increase trend for larger volumes. The structural analysis of the cliff showed four main joint sets that produced a parallelepiped unitary rockfall geometry. Thus, we can expect that larger rockfall events are a combination of several unitary parallelepipeds inducing a more complex geometry and thereby producing some concave geometries. In contrast, we also observed that volume computation influences a as computing less-complex geometries increases the overall volume size, which leads to an overestimation of volume size during the monitoring period.
This highlights the importance of measuring the complexity of the geometry of large rock falls and the use of an α shape to obtain a volume and a magnitude-frequency analysis that are closer to reality.
When looking at the volume-frequency plot in Figure 8 an underestimation of the number of rockfalls with small volumes (<0.05 m 3 ) is suspected. The same observation is made on Figure 7. This may have been caused by the limited spatial resolution (one point every 3.4 cm), as described in [35], and the minimal number of points used for cluster detection (see Section 3.2). Some large rockfalls were related to multiple events located in the same area (i.e., triangle symbol in Figure 8) but computed as a single event due to insufficient temporal resolution [35] or too close sources. This presence of multiple unresolved events can be resolved with a higher-frequency data acquisition, such as monthly or permanent acquisition [22], instead of seasonal acquisition. Thus, errors linked to coalescence and superposition of events can be reduced with enough temporal sampling [33].
The volume estimations provided the statistics of all volumes that provided information about the erosional processes or cliff retreat [58]. Based on the volume-frequency plot, the estimated exponent value (b) was between 0.48 and 0.78, and the factors were similar to those found in the literature by [34] for seaside sandstone cliffs. This can be correlated locally high fracturation in some sandstone layers or the poor GSI of the rock cliff but suggests that the volume computation and the regression line boundaries also influenced the erosional process interpretation [37,59]. The power law shows some limits in this study site as we have temporal resolution bias and effect related to variation of thickness in sandstone layers in a too small area leading to a too large effect on sampling. Nevertheless, using the power law, it is possible to extrapolate to large volumes that were not reached during the period of monitoring [58,60,61]. The interpretation made on the rockfall activity from the presented results allowed us to state that, even though the failures were clearly episodic, the mean cliff retreat rate by rockfall was estimated to be approximately 10 mm/year for a surface of 2240 m 2 in the zone of interest, assuming a mean detachment of material of 25 m 3 per year between October 2011 and October 2013 and of~21 m 3 per year over the full monitoring period from 2010 to 2015 on sandstone beds. In addition to these results, the erosional volume of marls was not considered in this study.

Conclusions
The proposed toolbox tries to fill an existing gap in the open-source tools available for the treatment of point clouds oriented toward the quantification of rockfall activity. This package is a compound of multiple state-of-the-art algorithms (DBSCAN, OPTICS, α shape) that are able to identify rockfall sources and their volumes by considering its complex shapes. Perspectives for the RockfallQuantification packages include an additional tool for describing the shape of each rockfall event in terms of the principal axis, which could bring interesting inputs to rockfall propagation software that takes into account the shape of the blocks, such as RocFall [62] or Rockyfor3D [63].
The application of the RockfallQuantification package to the La Cornalle cliff demonstrated a decrease in processing time and provided a robust procedure to quantify rockfall activity and determine erosion rate. We were able to compute volume with concavities observed in the field and to finally calculate the erosion rate and subsequent magnitudefrequency distribution of the local Molassic cliff. In addition, the results highlighted the importance of correctly estimating the volume of rockfall events when magnitude-frequency plots are used.
The RockfallQuantification package as well as the 3DPointCloudToolBox can also be applied to other high-resolution surface data techniques, such as aerial laser scanning or photogrammetric techniques, including Structure-from-Motion. In addition, the Rockfal-lQuantification package can be used for other applications in geosciences to extract volume information from other natural events (e.g., lava flow, sediment transport, etc.). Finally, as a free-code package, it allows users to constantly improve the package for use in their specific needs. The RockfallQuantification package contains a series of functions encapsulated in the main script ( Figure A2). This inner structure allows either extracting intermediary results at each step or running the full package to obtain the volume for each rockfall event as output. Figure A2 describe all functions used in the RockfallQuantification package that enable the extraction, individualization, and computation of rockfall volumes. As input data, the package needs two point clouds acquired at different epochs. Other necessary input parameters include the minimum number of points (k) to be considered to create a cluster in step two (see Section 3.2) and the research radius (r) needed to compute the shape of the rockfall point cloud at step three (see Section 3.3). It is possible to define all automatic research parameters (such as σ and ε) to decrease processing time or affine the results. It defines the threshold to determine σ at step one (see Section 3.1) or the neighborhood radius (ε) for step two (see Section 3.2). This solution is easier if the user has an idea of the research radius parameter, but this will increase the computation time. To compute the point-to-plane comparison, RockfallDetect needs the PointCloudComparison package. The provided final output is a database with labeled rockfall event volumes and associated point cloud objects defined at the second stage.
research parameters (such as σ and ε) to decrease processing time or affine the results. It defines the threshold to determine σ at step one (see Section 3.1) or the neighborhood radius (ε) for step two (see Section 3.2). This solution is easier if the user has an idea of the research radius parameter, but this will increase the computation time. To compute the point-to-plane comparison, RockfallDetect needs the PointCloudComparison package. The provided final output is a database with labeled rockfall event volumes and associated point cloud objects defined at the second stage. Figure A1. The structure of the 3DPointCloudToolBox contains folders labeled as Documents for tutorials, Examples of scripts to start up, Beta for functions in development, and Data to store raw data that can be used in the 3DPointCloudToolBox. The main folder Toolbox contains subfolders for point cloud data treatments (MainLibrary, PointCloudAlignment, PointCloudComparison, SyntheticPointCloud) and a specific subfolder landslide for landslide monitoring, with one for rockfall monitoring and another for landslide surface displacement monitoring (RockfallQuantification, LandslideTracking). Figure A1. The structure of the 3DPointCloudToolBox contains folders labeled as Documents for tutorials, Examples of scripts to start up, Beta for functions in development, and Data to store raw data that can be used in the 3DPointCloudToolBox. The main folder Toolbox contains subfolders for point cloud data treatments (MainLibrary, PointCloudAlignment, PointCloudComparison, Syn-theticPointCloud) and a specific subfolder landslide for landslide monitoring, with one for rockfall monitoring and another for landslide surface displacement monitoring (RockfallQuantification, LandslideTracking).

RockfallQuantification Functions
Step 1: RockfallExtract Extract point belonging to surface change from two PointCloud objects Step 2: RockfallSegment Individualize single rockfall event by clustering index on PointCloud

RockfallQuantification Functions
Step 1: RockfallExtract Extract point belonging to surface change from two PointCloud objects Step 2: RockfallSegment Individualize single rockfall event by clustering index on PointCloud dbscan_optics Density-Based Spatial Clustering of Applications with Noise [42] and OPTICS improvement [43] dist Compute Euclidean distance between points in the cloud epsilon Compute optimal epsilon radius according to gamma function approximation (Daszykowski et al., 2002) Step

Add
Add the content of a given point cloud to this one addNoise Add simulated noise to the true point positions with following possibility: Gaussian position smearing Outliers to simulate completely wrong position Drop out some points by replacing points position by NaNs

ComputeBoundaries
Compute the Boundary points

ComputeCurvature
Compute the curvatures at each point using: Estimation of the curvature based on [64] Variation of the surface from correlation of point clouds based on [65] ComputeDelaunayTriangulation Compute a 3D Delaunay triangulation using built-in MATLAB ® function

PointCloud
Constructor of the object PointCloud and related methods

AffinTransform
Apply an affine transformation to object PointCloud AlphaBoundary Determine the convex hull of the object PointCloud using [45] EuclDist Compute the Euclidean distance between two vectors of 3D points.

ComparePoint2Point
Compute comparison using the shortest point to point distance. Calculation is made using Euclidean distance between a given point in PointCloud A to the closest point in PointCloud B with output as absolute differences.

ComparePoint2Surface
Compute comparison using the shortest point to surface distance. Calculation is made between a given point in PointCloud A and the distance parallel to the normal to the closest point in PointCloud B with output as signed differences.