An Open-Source Algorithm for 3D ROck Slope Kinematic Analysis (ROKA)

: The Markland test is one of the most diffused and adopted methods of kinematic analysis for the identiﬁcation of critical intersections of rock discontinuities that could generate rock failures. Traditionally, the kinematic analysis is based on the use of a stereographic approach that is able to identify the critical combination between the orientations of discontinuities and the rock wall. The recent improvements in the use of Digital Outcrop Models (DOMs) created the conditions for the development of a new automatized approach. We present ROck Slope Kinematic Analysis (ROKA) which is an open-source algorithm aimed at performing the Kinematic Analysis using the discontinuity measures collected onto a 3D DOM. The presented algorithm is able to make a local identiﬁcation of the possible critical combination between the identiﬁed discontinuities and the orientation of the slope. Using this approach, the algorithm is able to identify on the slope the presence of critical combinations according to the traditional kinematic analysis of planar failure, ﬂexural toppling, wedge failure, and direct toppling modes of failures and then visualize them on DOMs. In this way, the traditional approach is more effective and can be adopted for a more detailed analysis of large and complex areas.


Introduction
Discontinuities have a great influence on rock slope stability because, along their surfaces, the cohesion of rock can be very small or null, and therefore, they can induce different Modes of Failure (MOFs). Rock slope stability can be analyzed using different approaches, from simple geometrical methods (e.g., kinematic analysis) to complex numerical methods (e.g., Limit Equilibrium (LE); Finite Element Method (FEM); Block Element Method (BEM); Discrete Element Method (DEM; FEM-DEM; etc.)).
The Kinematic Analysis (KA) is certainly the simplest and quickest technique because it is purely geometric. It evaluates the possible MOFs, analyzing the angular relationships between slope surface and discontinuities [1] without considering the forces responsible for the MOF [2]. In any case, before carrying out more sophisticated analyses by numerical methods, it is often useful to proceed with an evaluation of the potential mechanisms of kinematic instability affecting a slope. Several studies have been realized to determine the possible MOFs using stereographic projection-based techniques [3][4][5][6][7][8][9][10]. Among them, the Markland test [4] and its refinement proposed by Hocking [6] is one of the most used methods.
The main pitfalls based on the principle of the traditional Markland test are as follows [11]: (a) the assumption of a uniform slope angle and slope azimuth, (b) the possible presence of non-tightly clustered discontinuity data, and (c) the possible non-uniform spatial distribution of discontinuity orientations and sizes.
The assumption of a uniform slope angle and slope azimuth could negatively influence the KA of a natural slope with more than one oriented surface (multi-faced slope) [12]. As emphasized by [11], the traditional KA considers only the mean and most representative orientation of the identified discontinuity sets, due to the large number measures. This assumption makes the method qualitative and, anyhow, it is applicable only to tight discontinuity clusters. Moreover, the authors emphasize that the discontinuities orientation and size could change onto the rock slope surface. Therefore, some discontinuities with a specific orientation or dimension could be present only onto limited portions of the multi-faced rock slope. [13] show that the size of the discontinuities, not considered in the KA, has a crucial role in the possible formation of unstable intersections.
These limitations are principally related to a stereographic approach, where the KA is performed using only the mean attitudes of the rock slope and of the discontinuity sets.
Today, the computer technology permits overcoming these limits: the problem of the consideration of uniform slope angle, slope azimuth, and friction angle could be easily resolved performing the KA with the different orientations of the rock slope, using the method proposed by [12] or, as suggested by several authors (e.g., [11,[14][15][16][17]) performing a sensitivity analysis using different KAs calculated for each possible slope orientation; the problem of the presence of non-tightly clustered discontinuity data could be easily resolved considering all the orientations of the discontinuities (not only the representative) in order to consider all the orientation variability [11].
Notwithstanding, the spatial distribution of discontinuity orientations and sizes still remains a limitation because no proposed method considers the 3D location of the discontinuity (a certain set of discontinuities can be present only in a specific portion of the rock slope) and its dimension (small discontinuities normally have less influence on stability than larger ones).
This limit is principally related to techniques of data acquisition: with a traditional field survey, it is not possible to map all the discontinuities affecting a rock slope due to the inaccessibility of some portions of it, and it can be really difficult to define exactly the 3D position of the discontinuity measured by hand in the field.
The principal advantage of DOM-based rock slope stability analysis is the possibility of performing digital measurements, mapping all the visible discontinuities of the rock, and extracting their information, such as position, orientation, and size, with high precision and accuracy. The knowledge of the exact orientation and dimension of the measured discontinuities, but, in particular, of their location onto the rock slope, allows an easier verification if each discontinuity is critical for the portion of the slope where it outcrops.
In this paper, we present an open-source algorithm, written in MATLAB language, aimed at performing the KA using the discontinuity measures collected onto a 3D Digital Outcrop Model. This code performs the KA of four MOFs (planar failure, wedge failure, flexural toppling, and direct toppling) automatically for all the different portions of a rock slope, highlighting only the real possible MOFs caused by 3D mapped discontinuities (with a defined position, orientation, and size) that really intersect the calculated different portions of the slope. The main advantages of ROck Slope Kinematic Analysis (ROKA) are (i) the identification of the discontinuities responsible for the real possible MOFs in the different parts of the slope, (ii) the evaluation of their criticality and the visualization of the critical discontinuities, and (iii) the critical portions of the rock slope on the digital outcrop model in a 3D environment. In this way, KA can be easily performed even in complex environments to evaluate the local stability conditions, allowing focusing an eventual detailed in situ analysis only on limited portions of the slope.

The ROck Slope Kinematic Analysis (ROKA) Algorithm
The stability of the rock slopes is frequently controlled by the discontinuities affecting the rock mass. In such cases, the geometrical relationships between the structures and the slope's orientation determine the kinematic stability of the slope. Kinematic Analysis (KA) refers to the geometrically possible motion of blocks of rock without considering the forces involved but evaluating the allowed and constrained directions of movement.
Thanks to its simplicity and speed, the KA of a rock slope is generally used to assess whether masses of rock may move along geologic structures and slide down a slope. In particular, it can be used to evaluate four principal modes of potential sliding, which are also called Modes of Failure (MOFs): plane failure, wedge failure, flexural toppling, and direct toppling.
As already described, the traditional KA has some limitations, such as the impossibility of considering different slope orientations and taking into account the dimensions and positions of the rock discontinuities.
Here, we present an algorithm called ROck slope Kinematic Analysis (ROKA) and written in MATLAB language that is able to overcome these limits performing an ad hoc Kinematic Analysis (KA) of the possible Modes of Failure (MOFs) of each differently oriented part of a rock slope, considering only the discontinuities really measured onto the DOM, even if it is characterized by a complex geometry (i.e., considering all the different local orientations of the slope).
This code calculates the criticality of the discontinuities in a deterministic way, correlating their measured position, attitude, and size with the slope's local orientation and determining for every fracture the potential MOFs (planar failure, wedge failure, flexural toppling, and direct toppling).
ROKA has been organized into four main stages ( Figure 1):

1.
Preparation and definition of the input data (described in Section 2.2), which essentially consists of the point cloud representing the slope, the geometric information of the discontinuities, and the calculation of the variables (slope-discontinuity intersection cutoff, scan-radius, friction angle, and lateral limits); 2.
Preprocessing calculation to obtain the geometrical information of the real intersections between the discontinuities, described in Section 2.3. This step is fundamental for the KA of the wedge failure and direct toppling MOFs; 3.
Ad hoc KA described in Section 2.4, which represents the main core of the algorithm. After the two previous phases, the KA is performed in each position of the rock slope considering its local orientation (calculated using the scan-radius) and the discontinuities that really intersect the rock slope in the analyzed sector; 4.
Exportation of the results (potentially unstable discontinuities and portions of the rock slope) and their visualization in the open-source software Cloud Compare (Section 2.5); Cloud Compare is a 3D rendering environment where the 3D point cloud and meshes representing the rock slope can be visualized, and "virtual" discontinuity surfaces or traces can be measured [23].
x FOR PEER REVIEW 4 of 29

Assumptions of the Analysis
ROKA has been developed considering three main basic assumptions. The first and principal assumption is that the discontinuities' surface is considered planar and circular, following the Baecher's disc model [33]. This assumption is widely used in material mechanics [34][35][36][37], and it considers discontinuities as discs; therefore, a discontinuity can be described by the coordinates of the disc center (cx, cy, cz), disc radius (r), and orientation of the 3D plane in which the disc lies ( Figure 2). The orientation of the disc can be defined by its attitude (e.g., dip/dip direction) or by its normal N.

Assumptions of the Analysis
ROKA has been developed considering three main basic assumptions. The first and principal assumption is that the discontinuities' surface is considered planar and circular, following the Baecher's disc model [33]. This assumption is widely used in material mechanics [34][35][36][37], and it considers discontinuities as discs; therefore, a discontinuity can be described by the coordinates of the disc center (cx, cy, cz), disc radius (r), and orientation of the 3D plane in which the disc lies ( Figure 2). The orientation of the disc can be defined by its attitude (e.g., dip/dip direction) or by its normal N.

Assumptions of the Analysis
ROKA has been developed considering three main basic assumptions. The first and principal assumption is that the discontinuities' surface is considered planar and circular, following the Baecher's disc model [33]. This assumption is widely used in material mechanics [34][35][36][37], and it considers discontinuities as discs; therefore, a discontinuity can be described by the coordinates of the disc center (cx, cy, cz), disc radius (r), and orientation of the 3D plane in which the disc lies ( Figure 2). The orientation of the disc can be defined by its attitude (e.g., dip/dip direction) or by its normal N. The relationship between N, its components Nx, Ny, and Nz, and the dip/dip direction can be defined calculating the direction cosines [38]: The relationship between N, its components Nx, Ny, and Nz, and the dip/dip direction can be defined calculating the direction cosines [38]: The dip angle can be defined as: whereas the dip direction can be defined • if cos α > 0 and cos β > 0: • if cos α > 0 and cos β < 0 or cos α < 0 and cos β < 0: • if cos α < 0 and cos β > 0: The second important assumption is that all the considered discontinuities are only those visible and mapped onto the rock slope, and their surfaces are hypothesized as circular with a diameter equal to their maximum visible extension.
The third assumption is that the discontinuities or their intersection can be potentially unstable if they intersect the slope surface and satisfy these geometric relationships with the local orientation of the slope:

•
A discontinuity is potentially critical for planar failure if (6a) it is a dip slope discontinuity and (6b) its dip angle is higher than the friction angle φ and (6c) is lower than the apparent angle of the slope along the discontinuity dip direction; (6) • Two discontinuities are critical for wedge failure if (7a) they form a dip slope intersection and (7b) the dip angle of the intersection is higher than the friction angle φ and (7c) is lower than the apparent angle of the slope along the discontinuity dip direction; A discontinuity is critical for flexural toppling if (a) it is an anti-dip slope discontinuity and (b) its dip < 90 • − slope dip + friction angle φ; (a) anti dipdir slope − 90 • < dipdir discontinuity < anti dipdir slope • Two discontinuities are critical for direct toppling if (a) they form anti-dip slope intersection and (b) the dip of the intersection < 90 • − slope angle.
(a) anti dipdir slope − 90 • < trend intersection < anti dipdir slope This last assumption can partially change when the lateral limits are considered. In this case, assuming a general value of 20 • [5], the planar failure and flexural toppling MOFs can be activated by a discontinuity if they satisfy the previous assumptions and have a direction within ±20 • from the slope direction, whereas the direct toppling can be activated if two discontinuities form an intersection that satisfies the previous assumption and have a trend within ±20 • from the opposite dip direction of the slope.

Algorithm Input Data
The input required for the KA performed by the ROKA algorithm are as follows: 1.
The 3D point cloud representing the study rock slope; 2.
The geometric features (center coordinates, orientation, dimension/radius) and the set of attribution of the discontinuity; 3.
The dimension of the scan-sphere that will be centered on every point of the point cloud to select the surrounding points used to calculate the "average" local orientation of the slope; and the dimension of the circular scan-plane that will be centered on every point of the point cloud with the local orientation of the 3D model, which is used to calculate the intersections between the discontinuities and the slope surface; 4.
The cutoff value of the intersection between discontinuities and the moving scan-plane; 5.
The KA parameters, in particular the friction angle and the lateral limits.
The input point cloud must be a text file in which all the metric coordinates of the 3D points that represent the rock slope (x, y and z) and their normal vector components (Nx, Ny and Nz) are reported. To perform a correct analysis, the point cloud must have a near constant surface point density because it is required to ROKA calculations, especially for the definition of critical value of the discontinuity planes and intersection that could activate MOFs (described in Section 2.5).
The geometric information about the discontinuities must be stored in an Excel file (.xlsx) where the rows correspond to the discontinuities and the columns to their geometrical information. The XLSX file must report the header and the data as the example shown in Table 1. The set to which the discontinuity belongs must be stored in another Excel file (.xlsx) where dip, dip direction and set must be reported in three column for each discontinuity.
The dimension of the scan-radius can be decided considering several parameters, as the density of the point cloud representing the slope, the complexity of the slope geometry and the dimension of the discontinuities. In general, for each point of the 3D model the scan-radius defines ( Figure 3): 1.
the dimension of the scan-sphere (Figure 3b) that is used to calculate the 'average' local orientation of the slope;  The choice of the scan-radius size must be done considering also its influences on results, because modifying the scan-radius the different slope orientations and consequently the number of critical discontinuities and intersections can undergo variations.
The cutoff value of the length of the intersections between the discontinuities and the scan-plane must be defined in order to avoid an overestimation of the critical discontinuities. The cutoff must be expressed in meter and the default value selected (already considered by the algorithm) is equal to 90% of the scan-plane diameter.
The parameters required to perform the kinematic analysis are the friction angle and the lateral limits. The default values are 30° and 20°, respectively, because, as reported by [5], they are the most common values for many types of rocks slope. These values can be changed by the users.

Preprocessing, Real Discontinuity Intersections Calculation
Preliminarily, before the Kinematic Analysis (KA) process, the real discontinuity intersections are calculated taking into account the position and dimension of the discs representing the discontinuities ( Figure 4) and using the function 'IntersectionCalculator.m'. In particular, it consists of the following: 1. calculation of the intersection line between the two infinite planes that fit two discontinuities using the Hessian normal form [39,40]; 2. exclusion of the intersection lines that cannot be formed by the calculated circular discontinuities that have a finite dimension; 3. calculation and exportation of the dimension, orientation, and position of the possible "real" discontinuity intersections ( Figure 4). The choice of the scan-radius size must be done considering also its influences on results, because modifying the scan-radius the different slope orientations and consequently the number of critical discontinuities and intersections can undergo variations.
The cutoff value of the length of the intersections between the discontinuities and the scan-plane must be defined in order to avoid an overestimation of the critical discontinuities. The cutoff must be expressed in meter and the default value selected (already considered by the algorithm) is equal to 90% of the scan-plane diameter.
The parameters required to perform the kinematic analysis are the friction angle and the lateral limits. The default values are 30 • and 20 • , respectively, because, as reported by [5], they are the most common values for many types of rocks slope. These values can be changed by the users.

Preprocessing, Real Discontinuity Intersections Calculation
Preliminarily, before the Kinematic Analysis (KA) process, the real discontinuity intersections are calculated taking into account the position and dimension of the discs representing the discontinuities ( Figure 4) and using the function 'IntersectionCalculator.m'. In particular, it consists of the following: 1.
calculation of the intersection line between the two infinite planes that fit two discontinuities using the Hessian normal form [39,40]; calculation and exportation of the dimension, orientation, and position of the possible "real" discontinuity intersections ( Figure 4).

Ad Hoc Kinematic Analysis Process
After defining the input data and the preprocessing calculation, the algorithm performs the ad hoc Kinematic Analysis (KA) for each point of the point cloud representing the rock slope ( Figure 1).
In particular, for each point of the cloud, the algorithm performs the following: 1. Selects all the points of the model contained in the moving scan-sphere, which has a scan-radius defined by the user (see . If the scan plane has an overhanging attitude, the algorithm considers its attitude as vertical (dip equal to 90°) and facing toward the outside the rock slope.

ROKA Results and Their Visualization
The outputs of ROKA algorithm consist of (a) the identification of the discontinuities and intersections that could induce failures in all the different positions of the slope, (b) their classification based on their criticality, and (c) a stability classification of every point of the slope. (a) To indicate the discontinuity planes and intersections that could activate an MOF (the critical discontinuities), the ROKA algorithm exports them as DXF files into specific directories representing all the different MOFs that could be activated, allowing a further simple 3D visualization of these results. (b) To classify the critical discontinuities, ROKA considers the dimensions of their intersections with the slope and normalizes these values using the maximum detected

Ad Hoc Kinematic Analysis Process
After defining the input data and the preprocessing calculation, the algorithm performs the ad hoc Kinematic Analysis (KA) for each point of the point cloud representing the rock slope ( Figure 1).
In particular, for each point of the cloud, the algorithm performs the following: 1. Selects all the points of the model contained in the moving scan-sphere, which has a scan-radius defined by the user (see Figure 3a-c); 2.
Calculates the mean orientation of every portion of the slope using the normal vectors of the points ('normal2attitude.m' Matlab function) and creates a moving scan-plane (see Figure 3c,d); 3.
Verifies if the discontinuity planes and intersections cross the scan-plane ('intSPdisc.m' and 'intSPinters.m' Matlab functions), considering the defined cutoff value of the discontinuities; 4.
Performs the KA of the possible Mode Of Failure (planar failure, flexural toppling, wedge failure, or direct toppling) considering the scan plane's attitude and of the intersected discontinuity plane and intersection ('KinematicAnalysis.m' Matlab function). If the scan plane has an overhanging attitude, the algorithm considers its attitude as vertical (dip equal to 90 • ) and facing toward the outside the rock slope.

ROKA Results and Their Visualization
The outputs of ROKA algorithm consist of (a) the identification of the discontinuities and intersections that could induce failures in all the different positions of the slope, (b) their classification based on their criticality, and (c) a stability classification of every point of the slope.
(a) To indicate the discontinuity planes and intersections that could activate an MOF (the critical discontinuities), the ROKA algorithm exports them as DXF files into specific directories representing all the different MOFs that could be activated, allowing a further simple 3D visualization of these results. (b) To classify the critical discontinuities, ROKA considers the dimensions of their intersections with the slope and normalizes these values using the maximum detected value. Then, every discontinuity is ranked by its critical value on a scale between 0 and 1 (0-100%), where one is assigned to the critical discontinuity that forms the longest intersection with the slope. (4) the number of the discontinuities and (5) discontinuities intersections that intersect the local surface of the slope; (6) the critical value of the discontinuity planes that could cause planar failure or (7) flexural toppling considering the local orientation of the slope; (8) the critical value of the discontinuity intersections that could cause wedge failure; and (9) direct toppling.
The visualization of the results of the ROKA algorithm allows understanding which part of the rock slope can be more critical, considering the possible MOFs and the dimension of the critical discontinuities and indicating the areas of the rock slope where an eventual further field-based analysis must be focused (e.g., concerning a detailed characterization of discontinuity properties that cannot be defined with remote sensing techniques).

Case A: Northern Apennines-Calcareous Turbidites
Case study A is represented by a 50 m length road cut, which is located along the upper Staffora Valley's right side in the Northern Apennines (44.733 • latitude, 9.247 • longitude). The slope is made up by Monte Antola flysch formation (Early Campanian-Early Maastrichtian) composed by dominantly turbiditic calcareous graded beds, calcareous sandstones, sandstones and marlstones, with rare intervening thin, carbonate-free graded layers or uniformly fine-grained shale ( Figure 5).  (5) discontinuities intersections that intersect the local surface of the slope; (6) the critical value of the discontinuity planes that could cause planar failure or (7) flexural toppling considering the local orientation of the slope; (8) the critical value of the discontinuity intersections that could cause wedge failure; and (9) direct toppling. The visualization of the results of the ROKA algorithm allows understanding which part of the rock slope can be more critical, considering the possible MOFs and the dimension of the critical discontinuities and indicating the areas of the rock slope where an eventual further field-based analysis must be focused (e.g., concerning a detailed characterization of discontinuity properties that cannot be defined with remote sensing techniques).

Case A: Northern Apennines-Calcareous Turbidites
Case study A is represented by a 50m length road cut, which is located along the upper Staffora Valley's right side in the Northern Apennines (44.733°latitude, 9.247° longitude). The slope is made up by Monte Antola flysch formation (Early Campanian-Early Maastrichtian) composed by dominantly turbiditic calcareous graded beds, calcareous sandstones, sandstones and marlstones, with rare intervening thin, carbonate-free graded layers or uniformly fine-grained shale ( Figure 5). The thickness of the more competent beds varies from a few centimeters to tens of meters, with an average thickness of 1.5 to 2 m. With an E-W direction and south-facing, the road cut is affected by frequent small rockfalls (in general ≤1 m 3 ) that imperil the road SP-131.
This outcrop was studied using the Terrestrial Digital Photogrammetry approach with a Canon EOS 5D 12.8 Megapixel camera (Canon Inc., Tokyo, Japan) from very close distances (≈5 m). The thickness of the more competent beds varies from a few centimeters to tens of meters, with an average thickness of 1.5 to 2 m. With an E-W direction and south-facing, the road cut is affected by frequent small rockfalls (in general ≤1 m 3 ) that imperil the road SP-131.
This outcrop was studied using the Terrestrial Digital Photogrammetry approach with a Canon EOS 5D 12.8 Megapixel camera (Canon Inc., Tokyo, Japan) from very close distances (≈5 m).
The 166 full-resolution pictures (with a mean resolution of ≈1.4 mm/pixel) were elaborated by the Metashape software (Agisoft LLC., St. Petersburg, Russia) [41], and a 3D point cloud representing the outcrop with 97 million points and a mean points density of 199,000 pts/m 2 was developed ( Figure 6). The 166 full-resolution pictures (with a mean resolution of ≈1.4 mm/pixel) were elaborated by the Metashape software (Agisoft LLC., St. Petersburg, Russia) [41], and a 3D point cloud representing the outcrop with 97 million points and a mean points density of 199,000 pts/m 2 was developed ( Figure 6). The 3D model was georeferenced in local coordinates using Ground Control Points positioned onto an L-shape marker (Figure 7). The relative accuracy of the model was ±0.5° in attitude and 0.1% in length. The 3D model was georeferenced in local coordinates using Ground Control Points positioned onto an L-shape marker (Figure 7). The relative accuracy of the model was ±0.5 • in attitude and 0.1% in length.
If we consider that the error of the traditional manual techniques of measure is 3-5 • for the attitude [42], the obtained values can be regarded as fully acceptable [23,43].
The discontinuities considered by the ROKA calculation were mapped manually onto the 3D point cloud through the CloudCompare software (with the methodology described in [23,44,45]. This methodology consists of selecting the points onto the 3D model representing the discontinuity traces and planes and then estimating the orientation of their best-fit plane. A total of 1919 discontinuities were detected (Figure 8). The 3D model was georeferenced in local coordinates using Ground Control Points positioned onto an L-shape marker (Figure 7). The relative accuracy of the model was ±0.5° in attitude and 0.1% in length. If we consider that the error of the traditional manual techniques of measure is 3-5° for the attitude [42], the obtained values can be regarded as fully acceptable [23,43].
The discontinuities considered by the ROKA calculation were mapped manually onto the 3D point cloud through the CloudCompare software (with the methodology described in [23,44,45]. This methodology consists of selecting the points onto the 3D model representing the discontinuity traces and planes and then estimating the orientation of their best-fit plane. A total of 1919 discontinuities were detected ( Figure 8). showing the circular-shape considered for the discontinuities following the Baecher's disc model [33].
The features of the discontinuities detected onto the DOM (position, dimension and orientation) were exported into an XLSX file. Then, they were plotted onto a stereographic diagram to identify the main sets and perform the traditional KA.
The analysis reveals the presence of five main sets ( Figure 9): one related to bedding and the others related to fractures. showing the circular-shape considered for the discontinuities following the Baecher's disc model [33].
The features of the discontinuities detected onto the DOM (position, dimension and orientation) were exported into an XLSX file. Then, they were plotted onto a stereographic diagram to identify the main sets and perform the traditional KA.
The analysis reveals the presence of five main sets ( Figure 9): one related to bedding and the others related to fractures.
After the analysis of the discontinuities orientation, the stereographic projection was used for the KA of the main modes of failure, planar failure, flexural toppling, wedge failure, and direct toppling ( Figure 10).
The features of the discontinuities detected onto the DOM (position, dimension and orientation) were exported into an XLSX file. Then, they were plotted onto a stereographic diagram to identify the main sets and perform the traditional KA.
The analysis reveals the presence of five main sets ( Figure 9): one related to bedding and the others related to fractures.  After the analysis of the discontinuities orientation, the stereographic projection was used for the KA of the main modes of failure, planar failure, flexural toppling, wedge failure, and direct toppling ( Figure 10).
The traditional KAs was performed considering a mean slope orientation of 181°/56° (dip direction/dip), a friction angle of 30°, and lateral limits of 20°. The principal results show the following: • Part of the K4 discontinuity set and some random discontinuities could act as a sliding surface ( Figure 10a   The orientation of the intersection between K4 fractures and bedding could activate the direct toppling MOF (Figure 10d).
The number of the critical discontinuities and intersections calculated considering all the discontinuities as infinite planes are reported in Tables 2-4.  Before the application of the ROKA algorithm, the point cloud was subsampled to 599,973 points, obtaining a mean point density of 718 points/m 2 , and the discontinuities were modeled as circular, according to the Baecher's disc model ( [33] Figure 2), considering their orientation, position, and dimension. Finally, the "real" intersections were calculated (as shown in Figure 4).
Successively, using a moving scan-sphere with a defined scan-radius, the local orientation of the outcrop was calculated (Figure 11a,b), and its overhanging attitudes were defined (Figure 11c). 599,973 points, obtaining a mean point density of 718 points/m 2 , and the discontinuities were modeled as circular, according to the Baecher's disc model ( [33] Figure 2), considering their orientation, position, and dimension. Finally, the "real" intersections were calculated (as shown in Figure 4).
Successively, using a moving scan-sphere with a defined scan-radius, the local orientation of the outcrop was calculated (Figure 11a and b), and its overhanging attitudes were defined (Figure 11c). In this study, due to the small dimensions of the analyzed portion of the road cut, a 10-cm scan radius was used to obtain a representative "local" orientation of the slope, highlighting all the small changes.
Subsequently, the algorithm tested whether the discontinuities and their intersections intersect the moving circular scan-plane (Figure 12), representing the slope's local geometry. In this study, due to the small dimensions of the analyzed portion of the road cut, a 10-cm scan radius was used to obtain a representative "local" orientation of the slope, highlighting all the small changes.
Subsequently, the algorithm tested whether the discontinuities and their intersections intersect the moving circular scan-plane (Figure 12), representing the slope's local geometry. The middle of the outcrop shows no intersection because it is highly vegetated (see Figure 6) and, therefore, no visible discontinuities were mapped (Figure 8).
If the intersection test is positive, ROKA algorithm performs the Kinematic Analysis (KA) of the possible MOFs (planar failure, flexural toppling, wedge failure, and direct toppling- Figure 13) considering the local orientation of the slope and the attitude of the detected discontinuity planes and intersections. The middle of the outcrop shows no intersection because it is highly vegetated (see Figure 6) and, therefore, no visible discontinuities were mapped (Figure 8).
If the intersection test is positive, ROKA algorithm performs the Kinematic Analysis (KA) of the possible MOFs (planar failure, flexural toppling, wedge failure, and direct toppling- Figure 13) considering the local orientation of the slope and the attitude of the detected discontinuity planes and intersections.
tween the scan-plane and the discontinuity intersections. The middle of the outcrop shows no intersection because it is highly vegetated (see Figure 6) and, therefore, no visible discontinuities were mapped (Figure 8).
If the intersection test is positive, ROKA algorithm performs the Kinematic Analysis (KA) of the possible MOFs (planar failure, flexural toppling, wedge failure, and direct toppling- Figure 13) considering the local orientation of the slope and the attitude of the detected discontinuity planes and intersections. In Figure 13 (see percentage values), a classification of the different critical discontinuities, based on the length of their intersection with the slope (see Section 2.5), is re- In Figure 13 (see percentage values), a classification of the different critical discontinuities, based on the length of their intersection with the slope (see Section 2.5), is reported. In this way, the most unstable portions of the slope are highlighted, allowing eventually to focus the further detailed stability analyses (e.g., field-based roughness, aperture, filling, seepage, wall strength characterization and definition).
The results of algorithm KA are reported in Tables 4 and 5. The comparison of the results of the kinematic analyses performed by the traditional method and ROKA algorithm will be discussed further in the Discussion section.

Case B: Subvertical Rock Cliff of Quartzites-Ormea (CN, Italy), Ligurian Alps
The second study rock slope is a sub-vertical cliff that is composed mainly by quartzites and is located along the right side of the Tanaro Valley, in the Ligurian Alps (44.147 • latitude, 7.919 • longitude), close to the village of Ormea (CN, Italy). The rock cliff is ≈80 m high and ≈100 m wide, and it is characterized by frequent instability events (rock falls) that imperil some houses and infrastructures below it (road and bridge). For this reason, it was studied using the unmanned aerial vehicle-based digital photogrammetry.
The dataset of the slope was already presented and validated in [23]. The 3D point cloud has 98 million points, with a point density of 1000 pts/m 2 and a mean planimetric and vertical accuracy of 3.3 and 0.9 cm. One thousand thirty-six discontinuities were manually mapped (Figure 14). Three main sets of discontinuity ( Figure 15) were identified: the bedding and the K2 and K3 fracture sets.  Three main sets of discontinuity ( Figure 15) were identified: the bedding and the K2 and K3 fracture sets.  The traditional KA was performed using the orientations of mapped discontinuities. A mean attitude of the rock slope of 300 • /75 • (dip direction/dip) indicates that all the considered MOFs can be activated (Figure 16). The traditional KA was performed using the orientations of mapped discontinuities. A mean attitude of the rock slope of 300°/75° (dip direction/dip) indicates that all the considered MOFs can be activated ( Figure 16). In particular, the results reported in Tables 6 and 7 show that the most probable MOF is the wedge failure, with 17% of all the possible discontinuity intersections that could activate it.  In particular, the results reported in Tables 6 and 7 show that the most probable MOF is the wedge failure, with 17% of all the possible discontinuity intersections that could activate it. After the traditional KA, the information of the discontinuity mapped onto the 3D model and the point cloud representing the outcrop were formatted as required by the ROKA algorithm. Due to its big dimension and resolution (98 million points and 1000 pts/m 2 ), the point cloud was subsampled to 574,708 points (≈20 cm point spacing; ≈25 pts/m 2 ). Considering the visible dimension of the discontinuities represented by a Baecher's disc, 4667 "real" intersections (as described in Section 2.3) were identified.
Subsequently, the ROKA algorithm was performed using a scan-radius of 1 m (therefore with a scan-sphere and a scan-plane of around 4.19 m 3 and 3.14 m 2 , respectively). The different orientations of the slope calculated using the specified scan-radius are shown in Figure 17.  In Figure 18, the frequency of the discontinuities and discontinuities intersections are highlighted. In Figure 18, the frequency of the discontinuities and discontinuities intersections are highlighted.
The results of the ROKA algorithm KA for the intersected discontinuities and discontinuities intersections are indicated in Figures 19 and 20.   The statistics of the results obtained by the ROKA algorithm are reported in Tables 8  and 9.   The statistics of the results obtained by the ROKA algorithm are reported in Tables 8 and 9.

Discussion
The comparison between the results of the traditional kinematic analysis (Markland test) and the results of the ROKA analysis performed on two different rock slopes highlights the main advantages of the last methods.
The first is the strong reduction of the discontinuity intersections that must be verified. While traditional methods calculate all the possible intersections between the supposed infinite planes of discontinuity, the ROKA algorithm calculates only the "real" discontinuity intersections, considering not only their orientation but also their position and dimension as mapped onto the DOM. In the two examined study cases, this reduction is particularly evident; from 1,838,093 and 564,373 to 4349 and 4667 intersections for cases A and B, respectively, with a difference of more than 99%. Consequently, the ROKA algorithm detects a lower number of critical discontinuity intersections with respect to the traditional KA (Table 10). The second main advantage of the ROKA method is the more accurate determination of the different local slope orientations using a movable scan-sphere, the radius of which can be defined by the operator. On the contrary, the traditional KA generally considers only one orientation, i.e., the slope's mean orientation ( Figure 21).
The comparison between the results of the traditional kinematic analysis (Markland test) and the results of the ROKA analysis performed on two different rock slopes highlights the main advantages of the last methods.
The first is the strong reduction of the discontinuity intersections that must be verified. While traditional methods calculate all the possible intersections between the supposed infinite planes of discontinuity, the ROKA algorithm calculates only the "real" discontinuity intersections, considering not only their orientation but also their position and dimension as mapped onto the DOM. In the two examined study cases, this reduction is particularly evident; from 1,838,093 and 564,373 to 4349 and 4667 intersections for cases A and B, respectively, with a difference of more than 99%. Consequently, the ROKA algorithm detects a lower number of critical discontinuity intersections with respect to the traditional KA (Table 10). The second main advantage of the ROKA method is the more accurate determination of the different local slope orientations using a movable scan-sphere, the radius of which can be defined by the operator. On the contrary, the traditional KA generally considers only one orientation, i.e., the slope's mean orientation ( Figure 21). This different geometrical modeling of the slope influences the results of the stability analysis. The differences in the number of discontinuities that can cause planar failure and flexural toppling for the study cases A and B are synthesized in Table 11. Table 11. Comparison of the results of the different KA of planar and flexural toppling modes of failure. This different geometrical modeling of the slope influences the results of the stability analysis. The differences in the number of discontinuities that can cause planar failure and flexural toppling for the study cases A and B are synthesized in Table 11. The traditional KA tends to underestimate the number of critical discontinuity planes and the number of discontinuity sets because the different slope orientations are not considered (see Tables 2, 4, 6 and 8). For example, in the case A, whereas the traditional KA considers only the discontinuity set K2 as critical for planar failure (see Table 2), the ROKA algorithm considers critical also the sets K2, K3, K4, and K5 (see Table 4).
Another substantial advantage of the ROKA method consists of evaluating the intersection dimension between the slope and a critical discontinuity. In this way, it is possible to filter the results by discarding the smallest fractures and visualizing and classifying each critical discontinuity for every portion of the slope (Figures 22 and 23 The traditional KA tends to underestimate the number of critical discontinuity planes and the number of discontinuity sets because the different slope orientations are not considered (see Tables 2, 4, 6, and 8). For example, in the case A, whereas the traditional KA considers only the discontinuity set K2 as critical for planar failure (see Table 2), the ROKA algorithm considers critical also the sets K2, K3, K4, and K5 (see Table 4).
Another substantial advantage of the ROKA method consists of evaluating the intersection dimension between the slope and a critical discontinuity. In this way, it is possible to filter the results by discarding the smallest fractures and visualizing and classifying each critical discontinuity for every portion of the slope (Figures 22 and 23).  This can be an important advantage for the stability analysis of complex or very large natural outcrops, where the localization of the critical discontinuities and critical portion of the slope can be crucial. ROKA can highlight where to focus possible detailed field investigations, as the definition of the discontinuity parameters that the remote sensing techniques could not totally or partially assess (e.g., roughness, aperture, seepage, infilling, wall strength).
Big discontinuity datasets, such as those acquired by the automatic discontinuity planes detection algorithms (e.g., [46,47]), can be positively processed by ROKA. As shown in [23], automatic discontinuity planes detection algorithms require a difficult and timeexpensive validation process that could be done manually or automatically. The manual validation is a time-consuming process that increases the time needed for the automatic operation. ROKA permits simplifying the validation process with the possibility to detect and visualize, even in a very large dataset, only the discontinuities that could activate an MOF and filter them by their critical value. This can be an important advantage for the stability analysis of complex or very large natural outcrops, where the localization of the critical discontinuities and critical portion of the slope can be crucial. ROKA can highlight where to focus possible detailed field investigations, as the definition of the discontinuity parameters that the remote sensing techniques could not totally or partially assess (e.g., roughness, aperture, seepage, infilling, wall strength).
Big discontinuity datasets, such as those acquired by the automatic discontinuity planes detection algorithms (e.g., [46,47]), can be positively processed by ROKA. As shown in [23], automatic discontinuity planes detection algorithms require a difficult and time-expensive validation process that could be done manually or automatically. The manual validation is a time-consuming process that increases the time needed for the automatic operation. ROKA permits simplifying the validation process with the possibility to detect and visualize, even in a very large dataset, only the discontinuities that could activate an MOF and filter them by their critical value.

Conclusions
In this paper, a new open-source algorithm for the ROck slope Kinematic Analysis (ROKA) based on remotely sensed data is proposed and applied to two different case studies in two different locations in Italy, (a) Northern Apennines and (b) the Ligurian Alps, and two different geological contexts, (a) a carbonate flysch (Monte Antola formation, Early Campanian-Early Maastrichtian) and (b) and a quartzose sandstones and

Conclusions
In this paper, a new open-source algorithm for the ROck slope Kinematic Analysis (ROKA) based on remotely sensed data is proposed and applied to two different case studies in two different locations in Italy, (a) Northern Apennines and (b) the Ligurian Alps, and two different geological contexts, (a) a carbonate flysch (Monte Antola formation, Early Campanian-Early Maastrichtian) and (b) and a quartzose sandstones and conglomerate (Verrucano Lombardo formation, Late Permian and Ponte di Nava Quartzites, Early Triassic).
This method requires the modeling of the studied slope through a 3D point cloud [48], which can be acquired from terrestrial or aerial photogrammetry or laser scanner surveys, and the mapping on the Digital Outcrop Model (DOM) of all the visible discontinuities, which can be performed manually or automatically with the appropriate software (e.g., [46,47]). The surface of the discontinuities is assumed planar and circular, following the Baecher's disc model [33], of diameter equal to their maximum visible extension.
The main advantages and limitations of the ROKA algorithm can be synthesized as shown in Table 12.

Advantages
Limitations Only the "real" intersections between the 3D mapped discontinuities are considered, evaluating their effective position, orientation, and dimension ( Figure 4); All the "real" orientations of the slope are calculated, using a scan-sphere (Figure 3), and these local attitudes are used to perform the KA only when a discontinuity, or two discontinuities that form a wedge, intersect the local portion of the slope; The KA can be validated by visualizing the results onto the 3D model of the rock slope ( Figure 22); The critical discontinuities can be located in the 3D space and classified according to the extent of their intersection with the slope (Figure 23), giving great help to the user that must perform the stability analysis The Baecher's disc model considered by ROKA can lead to miscalculation of the possible dimensions and positions of the discontinuities intersections (wedge) due to the presence of discontinuities of different shapes; Pitfalls of the remote sensing technique used to build the slope 3D model (e.g., unsatisfying resolutions or occlusion effects) can cause inaccuracies in the detection and mapping of the discontinuities; The possible modes of failure (planar failure, wedge failure, flexural toppling, and direct toppling) are determined considering only the discontinuities and slope orientations, without considering other important features as aperture, roughness, filling, and seepages of the discontinuities.
Unlike traditional methods, the ROKA algorithm considers only the real intersections of the discontinuities with each other and with the slope; this procedure can evaluate the discontinuities' position and all the different orientations of the slope, obtaining more valid results. As already shown by some author [11,12,37], one of the main pitfalls of the traditional KA is the impossibility of considering all the spatial variation of the discontinuity sets that can strongly affect the results of the KA. With this procedure, the number of discontinuities to be evaluated is considerably reduced, while the relative number of discontinuities and sets potentially critical in the different portions of the slope increases.
Moreover, the proposed algorithm gives the possibility of visualizing the critical discontinuity planes and intersections onto the point cloud, allowing the validation of the results and highlighting the slope's main critical portions.
The ROKA algorithm can also help the evaluation and validation of the stability analysis results of automatic discontinuities detection methods that can produce many artefacts, as discontinuity planes that do not exist [23].
In conclusion, the use of remotely sensed discontinuity data and the presented opensource ROKA algorithm can improve the stability analysis of rock slope affected by a fracture network. In particular, mapping the discontinuities' position and orientation and calculating the local orientation of slope surface, the ROKA algorithm allows performing a specific kinematic analysis that indicates the "real" discontinuities and intersections critical for the stability, also identifying the most dangerous portions of the slope.
Future developments will concern: • The creation of a MATLAB Graphical User Interface (GUI) to make the algorithm more user-friendly; • The generation of a MATLAB standalone executable/application file to allow the use of the algorithm also on a computer without MATLAB installed; • The transcription of the algorithm in other coding languages (e.g., R, Python).