QDC-2D: A Semi-Automatic Tool for 2D Analysis of Discontinuities for Rock Mass Characterization

: Quantitative characterization of discontinuities is fundamental to deﬁne the mechanical behavior of discontinuous rock masses. Several techniques for the semi-automatic and automatic extraction of discontinuities and their properties from raw or processed point clouds have been introduced in the literature to overcome the limits of conventional ﬁeld surveys and improve data accuracy. However, most of these techniques do not allow characterizing ﬂat or subvertical outcrops because planar surfaces are difﬁcult to detect within point clouds in these circumstances, with the drawback of undersampling the data and providing inappropriate results. In this case, 2D analysis on the fracture traces are more appropriate. Nevertheless, to our knowledge, few methods to perform quantitative analyses on discontinuities from orthorectiﬁed photos are publicly available and do not provide a complete characterization. We implemented scanline and window sampling methods in a digital environment to characterize rock masses affected by discontinuities perpendicular to the bedding from trace maps, thus exploiting the potentiality of remote sensing techniques for subvertical and low-relief outcrops. The routine, named QDC-2D (Quantitative Discontinuity Characterization, 2D) was compiled in MATLAB by testing a synthetic dataset and a real case study, from which a high-resolution orthophoto was obtained by means of Structure from Motion technique. Starting from a trace map, the routine semi-automatically classiﬁes the discontinuity sets and calculates their mean spacing, frequency, trace length, and persistence. The fracture network is characterized by means of trace length, intensity, and density estimators. The block volume and shape are also estimated by adding information on the third dimension. The results of the 2D analysis agree with the input used to produce the synthetic dataset and with the data collected in the ﬁeld by means of conventional geostructural and geomechanical techniques, ensuring the procedure’s reliability. The outcomes of the analysis were implemented in a Discrete Fracture Network model to evaluate their applicability for geomechanical modeling.


Introduction
In the last two decades, remote sensing techniques have become widespread and have attracted many researchers for their application on rock mass investigations. The main goal is to overcome the limits of the conventional geostructural and geomechanical surveys to characterize rock masses, related to logistics, limited access, vegetation, adverse weather conditions, and human and instrument errors [1][2][3][4]. In this perspective, terrestrial/aerial laser scanning and photogrammetry are widely accepted because of their capability to acquire high-resolution point clouds in a reasonable time and in safe conditions, with relatively low costs [5][6][7]. Great attention has been paid to discontinuities in 3D models, covered by the Plio-Pleistocene deposits, is available at the site. In addition, 0.5 steps located in the eastern side of the rock mass were found to be essential for acterization of the bedding surfaces, which contribute to the formation of potent stable discrete blocks. Further, before implementing the study, field surveys wer out to ensure that Pietra Piatta site was representative of the whole study area.

Generation of the Digital Outcrop Model (DOM)
A DOM of the study site was produced by means of Unmanned Aerial Vehic acquisition and Structure from Motion (SfM) processing techniques. The UAV s the study site was carried out in five steps: (1) flight mission planning; (2) positioning and coordinates' acquisition of Ground Control Points (GCPs); (3) flight and image collection; (4) Structure from Motion (SfM) processing and generation of the dense point cl (5) building of the orthomosaic.
To achieve the optimal coverage of the investigated area, an automatic flig with front and side overlap of 75% and flight altitude of 18 m, was set in the p phase. The nadir photogrammetric survey was performed on 12 December 2019 morning, in order to avoid sunlight reflections from the sea, using a quadcopter DJI Inspire 2 equipped with a 20.8-Megapixel (MP) resolution camera, an integr ellite positioning system, and a remote flight controller. The system's specificat the details of the survey are summarized in Table 1. In addition, three GCPs wer ally positioned on the terrain before the flight such that they could be easily det the photos, and their coordinates were acquired by means of a Stonex SIII Di Global Positioning System. Geomorphology of the area is characterized by a series of marine terraces subparallel to the coastline [50,51], gently dipping to NE and linked by small scarps. They are, at places, carved by water ways locally called lame [52]: slightly incised in the limestone bedrock and flat-bottomed, typically dry, valleys that constitute the main drainage network during exceptional rainfall events. The mentioned morpho-structures are the result of the superimposition of the tectonic uplift of the Apulian platform and the absolute sea-level changes, starting from the middle Pleistocene [53]. Platforms and cliffs form the coast up to 20 m high, linked by embayments constituted by coastal erosion deposits (pocket beaches).
From a geologic standpoint, the site belongs to the eastern part of the Murge plateau, an emerged block of the Apulian foreland characterized by a 3-km-thick Cretaceous succession related to a wide carbonate platform, overlain by upper Pliocene to Lower Pleistocene transgressive deposits of shallow marine waters [54,55]. The lithofacies outcropping in the study area are composed of whitish to greyish limestones and dolostones belonging to the Calcare di Bari Fm., which are discontinuously overlain by yellowish calcarenites belonging to the Calcarenite di Gravina Fm. While the latter has a massive structure, the former lithofacies is constituted by thin to medium bedded layers, crossed by a network of subvertical discontinuities and locally folded.
The fracture pattern, together with the marine and karst processes [56,57] strongly modelling this sector of Apulia, contributes to the geomorphologic evolution of the coastal area: Detachment niches along the subvertical walls of the cliffs and boulders at their base, occasionally visible below the sea level, indicate local failures of the rock mass, and further Remote Sens. 2021, 13, 5086 4 of 28 potential instabilities cannot be excluded. These are among the most frequent geological hazards in coastal karst settings and are partly favored by the diffuse presence of karst conduits and caves, further weakening the carbonate rock mass [58][59][60][61].
Thus, a geostructural-geomechanical characterization of the site, with particular emphasis on the identification and characterization of the discontinuity sets and the estimation of the rock block shape, size, and volume, as well as the failure modes, is crucial for the management of prevention and mitigation measures at Polignano a Mare, one of the most important touristic sites of Apulia.
The rock mass considered for the 2D analysis is a 6400-m 2, low-relief platform developing from 10 m above the sea level to the current coastline, locally named Pietra Piatta. The choice of this area to perform the study was dictated by the lack of the Mediterranean vegetation and of anthropogenic elements as well, considering the disturbance produced by the residential center and by the vegetative area in nearby sectors. Due to erosional processes, the calcarenite facies crops out only in the external part of Pietra Piatta; therefore, a detailed view of the fracture traces in the Calcare di Bari Fm., which is generally covered by the Plio-Pleistocene deposits, is available at the site. In addition, 0.5-to 2-m steps located in the eastern side of the rock mass were found to be essential for the characterization of the bedding surfaces, which contribute to the formation of potentially unstable discrete blocks. Further, before implementing the study, field surveys were carried out to ensure that Pietra Piatta site was representative of the whole study area.

Generation of the Digital Outcrop Model (DOM)
A DOM of the study site was produced by means of Unmanned Aerial Vehicle (UAV) acquisition and Structure from Motion (SfM) processing techniques. The UAV survey at the study site was carried out in five steps: (1) flight mission planning; (2)  To achieve the optimal coverage of the investigated area, an automatic flight mode, with front and side overlap of 75% and flight altitude of 18 m, was set in the planning phase. The nadir photogrammetric survey was performed on 12 December 2019, at early morning, in order to avoid sunlight reflections from the sea, using a quadcopter platform DJI Inspire 2 equipped with a 20.8-Megapixel (MP) resolution camera, an integrated satellite positioning system, and a remote flight controller. The system's specifications and the details of the survey are summarized in Table 1. In addition, three GCPs were manually positioned on the terrain before the flight such that they could be easily detected on the photos, and their coordinates were acquired by means of a Stonex SIII Differential Global Positioning System.
The SfM technique was carried out using Agisoft Metashape Professional [62] to process the images and obtain a 3D rock mass model. During the photo importation phase, the software automatically detected the camera calibration and location parameters (camera focal length, coordinates of the image principal point, and lens' distortion coefficients). The images were georeferenced in a WGS84/UTM 33 N metric coordinate system. The three GCPs were semi-automatically identified on the photos and their position was validated by the operator. Taking into consideration potential imprecisions in the acquisition of the GPS coordinates from the drone, these GCPs were used as a constraint to optimize the georeferencing of the model. The Root Mean Square Errors of the GCPs are reported in Table 2. Successively, the SfM algorithm recognized multiple key points in each picture and matched them in the overlapping photos. Then, 248/248 photos were aligned with the "high-accuracy" alignment option and optimized by means of sparse bundle adjustment algorithm [63], while the key point matches (tie points) were positioned in a 3D environment, thus obtaining a sparse point cloud. Later, a Multi-View Stereo (MVS) algorithm was applied to generate a "high-quality" 3D dense point cloud (98,375,478 points). The dense point cloud was cleaned of unwanted elements such as points belonging to the sea water moving on the borders of the model, and objects used to perform the photogrammetric surveys, by means of segmentation process. We chose to remove the disturbance in this phase because the creation of masks on the unwanted elements in the preliminary phase of the SfM processing would have required a large amount of time, considering the high number of photos.
The final step consisted of the generation of a mesh of the model, using the "highquality" option available in Agisoft Metashape (9,785,754 faces) and in the extraction of a 2D orthomosaic from the mesh, with a 4.71-mm pixel size.

Conventional Characterization of Discontinuity Sets
Geostructural and geomechanical surveys were performed at Pietra Piatta site on 4 November 2020 to carry out a quantitative analysis of the discontinuities in the Calcare di Bari Fm. We chose to adopt window sampling techniques rather than scanline methods to avoid orientation biases, since discontinuities subparallel to the scanline are difficult to detect [8,64]. A preliminary visual inspection in the field helped to estimate the main discontinuity sets. Successively, two 100-m 2 -wide squares were created on the rock mass with a tape, to collect information on the discontinuities intersecting or contained in the survey window, following [65]. These areas, corresponding to two of the five sectors later analyzed with the MATLAB routine (sectors A and C in Figure 2), were selected for the presence of easily recognizable planar surfaces. Moreover, for each joint belonging to the analyzed discontinuity set, the strike, spacing, persistence, opening, and filling were measured by means of a Wilkie-type compass and a measuring tape. Roughness and wall strength were estimated using, respectively, a profilometer (Barton Comb) and a sclerometer (L-type Schmidt hammer). The strike data of the detected discontinuities were processed using the Dips software [66] by means of equal-angle stereographic projections (lower hemisphere) to identify the main DSs and calculate their statistical parameters. Successively, the geometrical data (i.e., spacing and trace length) collected in the field were processed on spreadsheets.
2D orthomosaic from the mesh, with a 4.71-mm pixel size.

Conventional Characterization of Discontinuity Sets
Geostructural and geomechanical surveys were performed at Pietra Piatta site November 2020 to carry out a quantitative analysis of the discontinuities in the Calca Bari Fm. We chose to adopt window sampling techniques rather than scanline met to avoid orientation biases, since discontinuities subparallel to the scanline are diffic detect [8,64]. A preliminary visual inspection in the field helped to estimate the main continuity sets. Successively, two 100-m 2 -wide squares were created on the rock mass a tape, to collect information on the discontinuities intersecting or contained in the su window, following [65]. These areas, corresponding to two of the five sectors later lyzed with the MATLAB routine (sectors A and C in Figure 2), were selected for the ence of easily recognizable planar surfaces. Moreover, for each joint belonging to th alyzed discontinuity set, the strike, spacing, persistence, opening, and filling were m ured by means of a Wilkie-type compass and a measuring tape. Roughness and strength were estimated using, respectively, a profilometer (Barton Comb) and a scle eter (L-type Schmidt hammer). The strike data of the detected discontinuities were cessed using the Dips software [66] by means of equal-angle stereographic projec (lower hemisphere) to identify the main DSs and calculate their statistical param Successively, the geometrical data (i.e., spacing and trace length) collected in the were processed on spreadsheets.

Manual Mapping of the Fracture Traces
The high-resolution orthophoto (1.25-GB-sized TIF image) produced by mea SfM technique was imported into a GIS environment by means of QGis open source ware, in a WGS 84/UTM 33N metric coordinate system. After a meticulous visual in tion of the model, five sectors were considered to be representative of the rock mass; t fore, a vector layer made up of five square windows was created on the orthophoto. cessively, for each sector, the traces of the discontinuities were manually drawn as ylines keeping the same scale, which was established according to the resolution o orthophoto, in order to reduce sampling biases in the different windows. The traces d on each polygon are reported in Table 3. The nodes of the polylines were extracted u

Manual Mapping of the Fracture Traces
The high-resolution orthophoto (1.25-GB-sized TIF image) produced by means of SfM technique was imported into a GIS environment by means of QGis open source software, in a WGS 84/UTM 33N metric coordinate system. After a meticulous visual inspection of the model, five sectors were considered to be representative of the rock mass; therefore, a vector layer made up of five square windows was created on the orthophoto. Successively, for each sector, the traces of the discontinuities were manually drawn as polylines keeping the same scale, which was established according to the resolution of the orthophoto, in order to reduce sampling biases in the different windows. The traces drawn on each polygon are reported in Table 3. The nodes of the polylines were extracted using the "extract vertices" geometry tool available in QGis, and their coordinates were obtained through the field calculator in the attribute table. The "clip" vector algorithm (input: layer of the traces, overlay layer = square window) allowed us to extract the traces of each sector, as well as their nodes coordinates. Finally, the five layers were exported into CSV format and converted to text files with a spreadsheet. The code for the 2D quantitative analysis of discontinuities was written in MATLAB environment after a preliminary literature review of the discontinuity properties that must be defined for geostructural and geomechanical characterization of rock masses and of the methods to define them as well. According to [65], the quantitative description of the discontinuities requires the determination of their orientation, spacing, persistence, roughness, wall strength, aperture, filling, and seepage, as well as the number of discontinuity sets and the block size. In this study, we focused on the determination of the geometrical properties and on the number of discontinuity sets from orthophotos, derived from remote sensing techniques by implementing the methods presented in the literature in a digital environment. The parameters taken into account to perform the 2D analysis, the methods, and the source publications of the formulas are summarized in Table 4. A synthetic dataset was created to gradually build the code. Three discontinuity sets (consisting of about 300 traces, in order to achieve statistical stability) were generated with strike and length following a normal distribution and a negative exponential law, respectively. The details of the synthetic dataset are illustrated in Figure 3. Successively, the 2D quantitative analysis was performed on the trace map (more than 1000 traces in sector C) of the case study and validated by comparing the results of the conventional geostructural and geomechanical surveys. To this aim, 50 discontinuities measured in the field in the same area were reported in Dips software [66] and classified into two main joint sets by means of stereographic projections (lower hemisphere). For each discontinuity set, the measurements of the spacings and trace lengths were elaborated in a spreadsheet to derive their mean values.
To test the potentiality of the 2D analysis for geological modelling, a Discrete Fracture Network model was generated on a 20-m-high cliff located at Lama Monachile site, about 90 m east from the study site. The top of the rock cliff is not visible because of the presence of the Calcarenite di Gravina Fm. and the overlying buildings. For this reason, the fracture pattern of the Calcare di Bari Fm. analyzed at Pietra Piatta site was used as an analogue model to generate the fracture sets by means of the FracMan software [72], after validation in the field of its representativeness for the cliff. The volume (grid) was generated from the top and bottom surfaces of the mesh of the rock mass using a stochastic distribution of the bedding surfaces, which was determined from measurements on a vertical scanline on the point cloud. Successively, for each discontinuity set, the mean strike, standard deviation, minimum, maximum, mean spacing, and trace length derived from the MATLAB routine were used to generate the fractures in the DFN model. 90 m east from the study site. The top of the rock cliff is not visible because of the pres of the Calcarenite di Gravina Fm. and the overlying buildings. For this reason, the fra pattern of the Calcare di Bari Fm. analyzed at Pietra Piatta site was used as an anal model to generate the fracture sets by means of the FracMan software [72], after valid in the field of its representativeness for the cliff. The volume (grid) was generated the top and bottom surfaces of the mesh of the rock mass using a stochastic distrib of the bedding surfaces, which was determined from measurements on a vertical sca on the point cloud. Successively, for each discontinuity set, the mean strike, standar viation, minimum, maximum, mean spacing, and trace length derived from the MAT routine were used to generate the fractures in the DFN model. igure 3. To the left, table illustrating the experimental dataset used to create the routine; to the right, graphical represenation of the three joint sets.

Workflow of the MATLAB Routine
The MATLAB routine performs the 2D quantitative analysis on orthophotos by bining digital scanline and areal sampling methods. In detail, rectangular/square dows [8,67] are applied to identify the main discontinuity sets according to the orient of the discontinuities, and to estimate their trace length and persistence as well. Cir windows [69,73] are performed to characterize the fracture network by identifyin mean trace length, intensity, and density estimators of the whole dataset. Areal sam is preferred to scanline methods to avoid sampling inaccuracies such as orientation length biases. In particular, orientation bias consists of the underestimation of the i sity of discontinuities, which are not perpendicular to a scanline [8,64,74], and length refers to the undersampling of short discontinuities with respect to the longer ones 77]. Since areal measures allow analyzing larger areas than scanlines [78], curtailmen caused by the loss of information of discontinuities extending beyond the sampling dows [8] is also reduced. However, scanline techniques are indispensable to defin properties, such as normal spacing and fracture frequency of the discontinuity sets [ Aimed at creating a tool for a 2D quantitative analysis of discontinuities, availab the scientific community and easily reproducible on different case studies, the MAT routines are accessible through an editable template.
The 2D quantitative analysis is achieved in four steps ( Figure 4): 1. graphical representation of the discontinuities; 2. semi-automatic classification of the discontinuity sets; 3. characterization of one discontinuity set; and 4. characterization of the fracture network.
The routine reads the first line of the template, in which the "STEP" comma defined. This allows performing the analyses described in Sections 2.6.1 and 2.6. simply adding 1, 2, 3, or 4 next to "STEP". Precise instructions for the users are in the

Workflow of the MATLAB Routine
The MATLAB routine performs the 2D quantitative analysis on orthophotos by combining digital scanline and areal sampling methods. In detail, rectangular/square windows [8,67] are applied to identify the main discontinuity sets according to the orientation of the discontinuities, and to estimate their trace length and persistence as well. Circular windows [69,73] are performed to characterize the fracture network by identifying the mean trace length, intensity, and density estimators of the whole dataset. Areal sampling is preferred to scanline methods to avoid sampling inaccuracies such as orientation and length biases. In particular, orientation bias consists of the underestimation of the intensity of discontinuities, which are not perpendicular to a scanline [8,64,74], and length bias refers to the undersampling of short discontinuities with respect to the longer ones [75][76][77]. Since areal measures allow analyzing larger areas than scanlines [78], curtailment bias caused by the loss of information of discontinuities extending beyond the sampling windows [8] is also reduced. However, scanline techniques are indispensable to define 1D properties, such as normal spacing and fracture frequency of the discontinuity sets [8,79].
Aimed at creating a tool for a 2D quantitative analysis of discontinuities, available to the scientific community and easily reproducible on different case studies, the MATLAB routines are accessible through an editable template.
The 2D quantitative analysis is achieved in four steps ( Figure 4): 1. graphical representation of the discontinuities; 2.
semi-automatic classification of the discontinuity sets; 3.
characterization of one discontinuity set; and 4. characterization of the fracture network.
The routine reads the first line of the template, in which the "STEP" command is defined. This allows performing the analyses described in Sections 2.6.1 and 2.6.2, by simply adding 1, 2, 3, or 4 next to "STEP". Precise instructions for the users are in the help and instruction files (Supplementary Material).  Step 1: workflow of the MATLAB routine developed for the 2D analysis of the discontinuities.

Characterization of the Discontinuity Sets
Orientation This property is commonly expressed in terms of dip (maximum inclination of a discontinuity to the horizontal) and dip direction (direction of the horizontal trace of the line of dip, measured clockwise from north). However, since this study dealt with discontinuity traces mapped on an orthophoto, the orientation was given by their strike (trace of the intersection of an inclined plane with a horizontal reference plane).
Two methods are used to identify and classify the discontinuity sets according to the orientation of the traces.

•
Method 1: fitting of the composite Gaussian curve The orientations of the polylines are calculated from the strike of their segments and weighted according to their lengths by means of a square/rectangular window applied on the selected dataset. The number of discontinuity sets is identified by means of an innovative method, which semi-automatically plots the strike distribution of the traces on an appositely developed Graphic User Interface (GUI) ( Figure 5). The MATLAB routine automatically extracts the distribution of each discontinuity set, with peaks corresponding to their mean strike. Three sliders allow adapting the strike distributions of each discontinuity set by interactively changing their mean, standard deviation, and amplitude. A line parallel to the x axis is used to filter out the "noise", derived from random discontinuities, which might be present in traces from real rock masses. After the manual adjustment of the curves, an optimization algorithm is run so that the sum of the curves and noise fits the raw orientation histogram curve. In detail, the best solution to fit the synthetic curve on the original curve is obtained by minimizing the differences between the curve of the real data and the experimental composite curve by means of the fminunc (Find minimum of unconstrained multivariable) algorithm, which is based on the BFGS Quasi-Newton method [80][81][82][83]: where the two terms are the orientation histogram of the dataset and the sum of the experimental curves, respectively. Successively, the MATLAB routine calculates the intersections between each normal curve, corresponding to the limits of the orientation of each joint set, to classify them. Therefore, for each DS, a file with the coordinates of the polylines' nodes is created. Step 1: workflow of the MATLAB routine developed for the 2D analysis of the discontinuities.

Characterization of the Discontinuity Sets Orientation
This property is commonly expressed in terms of dip (maximum inclination of a discontinuity to the horizontal) and dip direction (direction of the horizontal trace of the line of dip, measured clockwise from north). However, since this study dealt with discontinuity traces mapped on an orthophoto, the orientation was given by their strike (trace of the intersection of an inclined plane with a horizontal reference plane).
Two methods are used to identify and classify the discontinuity sets according to the orientation of the traces.

•
Method 1: fitting of the composite Gaussian curve The orientations of the polylines are calculated from the strike of their segments and weighted according to their lengths by means of a square/rectangular window applied on the selected dataset. The number of discontinuity sets is identified by means of an innovative method, which semi-automatically plots the strike distribution of the traces on an appositely developed Graphic User Interface (GUI) ( Figure 5). The MATLAB routine automatically extracts the distribution of each discontinuity set, with peaks corresponding to their mean strike. Three sliders allow adapting the strike distributions of each discontinuity set by interactively changing their mean, standard deviation, and amplitude. A line parallel to the x axis is used to filter out the "noise", derived from random discontinuities, which might be present in traces from real rock masses. After the manual adjustment of the curves, an optimization algorithm is run so that the sum of the curves and noise fits the raw orientation histogram curve. In detail, the best solution to fit the synthetic curve on the original curve is obtained by minimizing the differences between the curve of the real data and the experimental composite curve by means of the fminunc (Find minimum of unconstrained multivariable) algorithm, which is based on the BFGS Quasi-Newton method [80][81][82][83]: where the two terms are the orientation histogram of the dataset and the sum of the experimental curves, respectively.
Successively, the MATLAB routine calculates the intersections between each normal curve, corresponding to the limits of the orientation of each joint set, to classify them. Therefore, for each DS, a file with the coordinates of the polylines' nodes is created.

Figure 5.
Step 1: Optimization process for the semi-automatic identification of the discontinuity sets forming the synthetic dataset by means of Gaussian composite curve fitting. The classification limits for each discontinuity set are represented by the black lines. The red and orange lines represent the smoothed curve and the estimation of the three joint sets, respectively; the purple, green, and blue curve represent J1, J2, and J3, respectively.

Method 2: Hough transform
A second method to classify the discontinuity sets is based on the Hough transform, a technique introduced in [84] for machine analysis of bubble chamber photographs and later extended and adapted to image analysis and computer vision to detect shapes such as lines, circles, and ellipses in images [85][86][87][88]. Based on the Hough transform method, the MATLAB routine converts each joint in a new frame with the x and y axes corresponding to the orientation (θ) and distance from the origin (r) of each polyline, respectively ( Figure  6). In this new frame, each joint is defined by a point, and joints with similar orientation are aligned on the vertical axis, allowing us to identify the main discontinuity sets ( Figure  7).  Step 1: Optimization process for the semi-automatic identification of the discontinuity sets forming the synthetic dataset by means of Gaussian composite curve fitting. The classification limits for each discontinuity set are represented by the black lines. The red and orange lines represent the smoothed curve and the estimation of the three joint sets, respectively; the purple, green, and blue curve represent J1, J2, and J3, respectively.

Method 2: Hough transform
A second method to classify the discontinuity sets is based on the Hough transform, a technique introduced in [84] for machine analysis of bubble chamber photographs and later extended and adapted to image analysis and computer vision to detect shapes such as lines, circles, and ellipses in images [85][86][87][88]. Based on the Hough transform method, the MATLAB routine converts each joint in a new frame with the x and y axes corresponding to the orientation (θ) and distance from the origin (r) of each polyline, respectively ( Figure 6). In this new frame, each joint is defined by a point, and joints with similar orientation are aligned on the vertical axis, allowing us to identify the main discontinuity sets (Figure 7).

Spacing
The normal spacing of the discontinuity set (i.e., distance between two adjacent discontinuities belonging to the same set measured along a sampling line orthogonal to the mean direction of the set) can be calculated with two different methods according to their persistence. In a 2D analysis, the discontinuity persistence, which is the areal extent or size of a discontinuity within a plane [51], can be expressed as the limit length ratio along a given line on a joint plane [89]: where L s is the length of a straight line segment S and l s i is the length of the i-th joint segment in S. In other words, persistent joints are represented by continuous traces, while non-persistent joints are formed by more segments separated by rock bridges.
later extended and adapted to image analysis and computer vision to detect shapes such as lines, circles, and ellipses in images [85][86][87][88]. Based on the Hough transform method, the MATLAB routine converts each joint in a new frame with the x and y axes corresponding to the orientation (θ) and distance from the origin (r) of each polyline, respectively ( Figure  6). In this new frame, each joint is defined by a point, and joints with similar orientation are aligned on the vertical axis, allowing us to identify the main discontinuity sets ( Figure  7).

Spacing
The normal spacing of the discontinuity set (i.e., distance between two adjacent discontinuities belonging to the same set measured along a sampling line orthogonal to the mean direction of the set) can be calculated with two different methods according to their persistence. In a 2D analysis, the discontinuity persistence, which is the areal extent or size of a discontinuity within a plane [51], can be expressed as the limit length ratio along a given line on a joint plane [89]: where Ls is the length of a straight line segment S and l is the length of the i-th joint segment in S. In other words, persistent joints are represented by continuous traces, while non-persistent joints are formed by more segments separated by rock bridges.
• Spacing for non-persistent joints After importing the file of one DS, the scanline method is applied to derive the normal spacing and frequency, with the assumption of non-persistent discontinuities. One or more scanlines can be plotted according to the user's setting. Specifically, one linear scanline can be automatically plotted in the middle of the frame orthogonal to the mean strike

•
Spacing for non-persistent joints After importing the file of one DS, the scanline method is applied to derive the normal spacing and frequency, with the assumption of non-persistent discontinuities. One or more scanlines can be plotted according to the user's setting. Specifically, one linear scanline can be automatically plotted in the middle of the frame orthogonal to the mean strike of the discontinuity set. To find the best scanline, the routine plots several randomly oriented scanlines and selects the one intersecting the highest number of joints. Alternatively, the user can pick the two endpoints or pick one endpoint and set the mean strike of the scanline directly on the trace map. In addition, a series of scanlines parallel to the reference one can be traced, with the possibility to choose the number of samples and the ∆x-∆y interval. Considering the spatial variation of the strike of the traces, even if belonging to the same discontinuity set, Terzaghi's correction [64] is applied to each line intersecting the scanline to avoid orientation bias: where S, S app , and θ are the corrected spacing, apparent spacing, and the minimum angle between the scanline and the mean strike of the discontinuity set. After the calculation of the apparent (not corrected) and normal spacings, the histogram and cumulative distributions of the spacings are plotted on the interface. The normal spacing is then calculated as the mean distance of the intersections between the scanline and the traces (Figures 8 and 9). The mean fracture frequency is obtained as the inverse of the mean normal spacing.
Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 28 The normal spacing is then calculated as the mean distance of the intersections between the scanline and the traces (Figures 8 and 9). The mean fracture frequency is obtained as the inverse of the mean normal spacing. Figure 8. Normal set spacing calculated along a scanline (red) perpendicular to the mean direction of the joint set (black) (modified after [4]).

•
Spacing for persistent joints • The Hough transform method is applied to calculate the spacing of persistent joints, assuming infinite lengths ( Figure 6). For each two adjacent joints, the spacing is calculated as the mean difference along the r axis. In the Cartesian frame, the mean spacing is given by: where S2 and S1 are the orthogonal distances between adjacent traces and θ2 and θ1 are the angles between the joints and the scanline. The normal spacing is then calculated as the mean distance of the intersections between the scanline and the traces (Figures 8 and 9). The mean fracture frequency is obtained as the inverse of the mean normal spacing. Figure 8. Normal set spacing calculated along a scanline (red) perpendicular to the mean direction of the joint set (black) (modified after [4]).

•
Spacing for persistent joints • The Hough transform method is applied to calculate the spacing of persistent joints, assuming infinite lengths ( Figure 6). For each two adjacent joints, the spacing is calculated as the mean difference along the r axis. In the Cartesian frame, the mean spacing is given by: where S2 and S1 are the orthogonal distances between adjacent traces and θ2 and θ1 are the angles between the joints and the scanline. • Spacing for persistent joints The Hough transform method is applied to calculate the spacing of persistent joints, assuming infinite lengths ( Figure 6). For each two adjacent joints, the spacing is calculated as the mean difference along the r axis. In the Cartesian frame, the mean spacing is given by: where S 2 and S 1 are the orthogonal distances between adjacent traces and θ 2 and θ 1 are the angles between the joints and the scanline.

Trace Length
The mean trace length of the discontinuity set is obtained by calculating the mean length of the polylines in the dataset. In addition, the histogram and cumulative frequency of the trace lengths are plotted on the interface.

Persistence
The method proposed in [67] is used to estimate the persistence of the discontinuity set by plotting a rectangular window on the trace map. The percentage of coverage area can be chosen by the user or drawn on the trace map. The discontinuities transecting the window n 2 (two intersections between the polyline and the window), contained in the window (n 0 ) (both endpoints are located in the window), and the total discontinuities (n tot ) are automatically counted, so that the mean persistence is calculated as: where w and h are, respectively, the length and height of the window and Φ is the acute angle between the mean strike of the discontinuity set and the height of the window ( Figure 10). The persistence is calculated in two steps: (1) over the entire dataset ( Figure 11a); (2) by dividing the dataset with a grid and calculating the persistence for each element of the grid (Figure 11b,c).

Figure 9.
Step 3: linear scanline on JS3. (a) automatic generation of a linear scanline (black, dotted line) pe the mean strike of the discontinuity set; the discontinuities and their intersections with the scanline are, res resented by the blue lines and the red crosses; (b) rose diagram of the JS; (c) table of the calculated parame grams of the calculated parameters; (e) cumulative distributions of the calculated parameters.

Trace Length
The mean trace length of the discontinuity set is obtained by calc length of the polylines in the dataset. In addition, the histogram and cum of the trace lengths are plotted on the interface.

Persistence
The method proposed in [67] is used to estimate the persistence of set by plotting a rectangular window on the trace map. The percentage can be chosen by the user or drawn on the trace map. The discontinuiti window n2 (two intersections between the polyline and the window), window (n0) (both endpoints are located in the window), and the total di are automatically counted, so that the mean persistence is calculated as μ = where w and h are, respectively, the length and height of the window a angle between the mean strike of the discontinuity set and the height of ure 10). The persistence is calculated in two steps: (1) over the entire dataset ( Figure 11a); (2) by dividing the dataset with a grid and calculating the persistence of the grid (Figure 11b,c).
The results of both methods are provided by the routine.

Characterization of the Fracture Network
Intensity, Density, and Trace Length Estimators Information on the fracture network is collected through the method Figure 10. Method for persistence calculation (modified after [90]).
spacing, trace length, and persistence. The mean intensity, density, and trace length estimators of the synthetic dataset, calculated by means of circular windows, are 177.15 m -1 , 25.38 m -2 , and 1.14 m, respectively. A graphical representation of the results is given by means of intensity and density maps ( Figure 12).
The block volume, volumetric joint count, and shape factor were calculated as 0.01 m 3 , 33.81 m -1 , and 40.58, according to the mean strike and normal set spacings. Moderately flat blocks were identified on the diagram in Figure 13.  The results of both methods are provided by the routine.

Characterization of the Fracture Network Intensity, Density, and Trace Length Estimators
Information on the fracture network is collected through the method proposed in [91] by means of circular window sampling. A grid of circles is plotted on the polylines to compute the mean trace length P 11 , the intensity P 21, and the density P 20 estimators: Σm Σr 2 (6) where n is the mean number of intersections between the polylines and the circles and m and m are the total and mean number of endpoints of the polylines contained in the circles, respectively. Additionally, a persistence map can be plotted with the possibility of choosing the number of circles, with each pixel representing the persistence calculated in one circle.

Block Volume and Shape
Formulas illustrated in [70,71] are used to derive the block volume V B , volumetric joint count J v , and block shape factor β of the blocks delimited by discontinuities. As stated in [70,71], rock blocks are formed from the intersection of at least three discontinuity sets with different directions. Since rock volumes cannot be determined by means of 2D analyses from plan views, a third dimension is needed to apply this procedure. If only two joint sets are detected on the trace map, the third dimension can be represented by the bedding. If the rock mass is characterized by more than three discontinuity sets, a rough estimation of the volumes can be achieved by considering the prevailing ones.
The block volume can be estimated from the mean spacings of the discontinuity sets: where S1, S2, S3 are the mean spacings of the three discontinuity sets delimiting the rock volume, and γ 1 , γ 2 , γ 3 are the angles between the discontinuity sets. Since the method is applied to joint sets perpendicular to the strata (γ 1 = γ 2 = 90 • ), only the angle between the two joint sets (γ 3 ) is taken into account.
For rock masses characterized by only one or two discontinuity sets, rock blocks can be formed when additional random joints cross the volume. In this case, the equivalent block volume can be calculated as: V B ≈ 50 S1 3 for only one joint set with mean spacing S1 (8) V B ≈ 5 S1 2 S2 for two joint sets with mean spacing S1 and S2 (9) The volumetric joint count is calculated as the number of joints per unit volume: where S i is the mean spacing of the i-th discontinuity set. The third dimension is needed to estimate J v . The parameters α 2 (medium spacing/smallest spacing) and α 3 (largest spacing/smallest spacing) are automatically calculated to determine the block shape factor β and to plot it on the chart shown in [70,71] in order to classify the block shape:

Results of the Synthetic Dataset
The synthetic dataset was classified into three discontinuity sets by means of the Hough Transform and composite Gaussian curve fitting (Tables 5 and 6).  Figures 9 and 11 illustrate, respectively, the results of the scanline and rectangular window methods on the discontinuity set JS3 to calculate the mean orientation, normal spacing, trace length, and persistence.
The mean intensity, density, and trace length estimators of the synthetic dataset, calculated by means of circular windows, are 177.15 m -1 , 25.38 m -2 , and 1.14 m, respectively. A graphical representation of the results is given by means of intensity and density maps ( Figure 12).
The block volume, volumetric joint count, and shape factor were calculated as 0.01 m 3 , 33.81 m -1 , and 40.58, according to the mean strike and normal set spacings. Moderately flat blocks were identified on the diagram in Figure 13.  Step 4: semi-automatic estimation of block shape and volume of the experimental dataset, according to formulas proposed in [70,71].

Application to the Case Study
The results of the mean strike, spacing, trace lengths, and persistence, both for the MATLAB routine and field surveys, are shown in Figures 14 and 15. Three scanlines at different positions were traced on both the trace maps of JS1 and JS2 to investigate the effect of scanline positions on the estimation of the normal set spacing (Figure 16), as well as the fracture density and intensity maps (Figure 17). A sensitivity analysis was performed to estimate the effects of the sample size and number of circles to calculate the

Application to the Case Study
The results of the mean strike, spacing, trace lengths, and persistence, both for the MATLAB routine and field surveys, are shown in Figures 14 and 15. Three scanlines at different positions were traced on both the trace maps of JS1 and JS2 to investigate the effect of scanline positions on the estimation of the normal set spacing (Figure 16), as well as the fracture density and intensity maps (Figure 17). A sensitivity analysis was performed to estimate the effects of the sample size and number of circles to calculate the persistence, trace length intensity, trace intensity, and trace density ( Figure 18). An equivalent block volume, equal to 0.03 m 3 , was obtained from the spacings of the two main discontinuity sets and the mean layer thickness (S0, calculated on exposed subvertical surfaces in the point cloud of the rock mass), using Equation (6). The reported volumetric joint count is 9.94 m -1 , while the block shape factor (β = 28.98) corresponds to a compact, slightly flat block shape, in agreement with field observations. The results of the Matlab routine are reported in Table 7. discontinuity sets and the mean layer thickness (S0, calculated on exposed subvertical surfaces in the point cloud of the rock mass), using Equation (6). The reported volumetric joint count is 9.94 m -1 , while the block shape factor (β = 28.98) corresponds to a compact, slightly flat block shape, in agreement with field observations. The results of the Matlab routine are reported in Table 7.
Moreover, the routine was applied on the trace maps of sectors A, B, D, and E to identify potential differences in the number of discontinuity sets, as well as their mean strikes. Figure 19 shows the differences between the mean joint sets detected in sector C and sector B.
The statistics of strike, spacing, and trace length of the joint sets extracted from the MATLAB routine ( Figure 20) were used to generate a realistic Discrete Fracture Network model of the Calcare di Bari Fm along the cliff located at Lama Monachile site. The data obtained from sector C were considered to be representative for the cliff after validation in the field (Figure 21).  discontinuity sets and the mean layer thickness (S0, calculated on exposed subvertical surfaces in the point cloud of the rock mass), using Equation (6). The reported volumetric joint count is 9.94 m -1 , while the block shape factor (β = 28.98) corresponds to a compact, slightly flat block shape, in agreement with field observations. The results of the Matlab routine are reported in Table 7.
Moreover, the routine was applied on the trace maps of sectors A, B, D, and E to identify potential differences in the number of discontinuity sets, as well as their mean strikes. Figure 19 shows the differences between the mean joint sets detected in sector C and sector B.
The statistics of strike, spacing, and trace length of the joint sets extracted from the MATLAB routine ( Figure 20) were used to generate a realistic Discrete Fracture Network model of the Calcare di Bari Fm along the cliff located at Lama Monachile site. The data obtained from sector C were considered to be representative for the cliff after validation in the field (Figure 21).          Moreover, the routine was applied on the trace maps of sectors A, B, D, and E to identify potential differences in the number of discontinuity sets, as well as their mean strikes. Figure 19 shows the differences between the mean joint sets detected in sector C and sector B.
The statistics of strike, spacing, and trace length of the joint sets extracted from the MATLAB routine ( Figure 20) were used to generate a realistic Discrete Fracture Network model of the Calcare di Bari Fm along the cliff located at Lama Monachile site. The data obtained from sector C were considered to be representative for the cliff after validation in the field (Figure 21).

Discussion
With regards to the synthetic dataset, the Gaussian fitting method recognized 3/3 discontinuity sets, with a difference of 2° for the strike of each discontinuity set and up to 1° for the standard deviation ( Figure 5). On the contrary, some difficulties were found for the classification by means of the Hough transform method, because discontinuities with too similar orientations create unprecise clusters in the theta-r diagram (Figure 7). We found out that the histogram fitting method for classification is reliable for a high number of traces (at least 50), since it is based on statistic procedures, regardless of the orientations of the discontinuity sets. The Hough transform method is recommended for lower numbers of discontinuities, provided that the mean strikes of the sets are not too close ( Figure  22).

Discussion
With regards to the synthetic dataset, the Gaussian fitting method recognized 3/3 discontinuity sets, with a difference of 2 • for the strike of each discontinuity set and up to 1 • for the standard deviation ( Figure 5). On the contrary, some difficulties were found for the classification by means of the Hough transform method, because discontinuities with too similar orientations create unprecise clusters in the theta-r diagram (Figure 7). We found out that the histogram fitting method for classification is reliable for a high number of traces (at least 50), since it is based on statistic procedures, regardless of the orientations of the discontinuity sets. The Hough transform method is recommended for lower numbers of discontinuities, provided that the mean strikes of the sets are not too close ( Figure 22).

Discussion
With regards to the synthetic dataset, the Gaussian fitting method recognized 3/3 discontinuity sets, with a difference of 2° for the strike of each discontinuity set and up to 1° for the standard deviation ( Figure 5). On the contrary, some difficulties were found for the classification by means of the Hough transform method, because discontinuities with too similar orientations create unprecise clusters in the theta-r diagram (Figure 7). We found out that the histogram fitting method for classification is reliable for a high number of traces (at least 50), since it is based on statistic procedures, regardless of the orientations of the discontinuity sets. The Hough transform method is recommended for lower numbers of discontinuities, provided that the mean strikes of the sets are not too close ( Figure  22).  Based on these observations, we processed the traces of the study site with the Gaussian fitting method and achieved good results with respect to the geostructuralgeomechanical analysis carried out on the study site. The mean strikes of the discontinuity sets J1 and J2 obtained from the MATLAB routine were 34 • and 124 • . These results were in agreement with the data collected by means of field surveys (32 • and 123 • for J1 and J2, respectively). The same standard deviation value of the strike of J2 was found both for the MATLAB routine and for the on-site surveys, while a difference of about 4 • was detected for J1. The standard deviations differed 4.2 • for J1 and 0 • for J2. This difference can be related to a major dispersion of the poles of the discontinuities, which was not detected in the field because of the lower number of sampled discontinuities (50).
With regards to the spacing estimations, the mean values were 6 cm (J1) and 2 cm (J2), higher than the data measured in the field. Several tests were carried out in order to obtain results in agreement with the field observations because the calculated mean spacing can vary significantly with the position and strike of the scanline (Figure 16). In this perspective, we remark that it is fundamental to use a proper dataset as input: unmapped traces in the orthophoto or objects/vegetation covering the joints could determine fewer intersections between the traces and the scanline, leading to an overestimation of the spacing. In addition, the possibility to apply linear samples in different locations of the dataset may help to understand the spatial variability of the calculated parameter, to identify potential changes of the stress field, or to detect sectors with different geomechanical behavior.
It is remarkable that both the spacing and trace lengths of J1 and J2 followed an exponential negative distribution, indicated by very good R-Squared values (in the range of 0.91-0.98, Figure 20), as found by [8,75,92,93] for a variety of rock masses.
A difference of 0.5 m was found in the estimation of the mean persistence of J2. It is believed that this difference is attributable to the different sampling methods for its estimation. Indeed, the persistence of discontinuities is one of the most difficult parameters to estimate in the field [65] and, during the conventional geostructural surveys, this property was roughly approximated from the trace lengths, according to [8]. However, this approach can lead to size and censoring bias [79,94], thus undersampling smaller discontinuities and censoring surfaces, which extend beyond the sampling area. Indeed, since the traces belonging to J1 were not truncated or censored because both the majority of their terminations fell within the area sampled in the field, no difference was found with respect to the results of the MATLAB routine. Based on these observations, it is believed that the results of the MATLAB routine are more reliable because the estimation of the persistence by means of [67]'s method takes into account the number of discontinuities transecting and contained in the sampling window, thus avoiding censoring bias. However, care should be taken when choosing the size of the coverage area. As pointed out by [8], this method is not suitable when all the discontinuities transect the window (the persistence would be infinite) and when no discontinuities are contained in the window (the persistence would be zero). For this reason, the sampling window should be chosen so that at least one trace is located inside it. In addition, we observed how the persistence varied with different dimensions of the sampling area (Figure 18a): Although two datasets are not enough to identify potential mathematical relations, it is evident that significantly different results were obtained by changing the area covered by the sampling window, especially for long traces (e.g., J2 curve in Figure 18a). As a matter of fact, [79] observed that persistence can vary up to more than 50% of the real value by changing the rectangular sampling window location and size.
The additional parameters calculated from the routine such as mean intensity, density, and trace length estimators give useful information for the identification of more fractured zones, represented by the yellow pixels in Figure 17c. Moreover, the fracture abundance of one discontinuity set constituted by non-parallel, subparallel, or non-persistent traces can be expressed by means of fracture intensity [69]. In fact, the estimation of fracture abundance for non-parallel traces by means of spacing is rather unclear because the distance is not unique. With regards to the length, intensity, and density estimators, sensitivity analyses allowed us to observe that more accurate results can be achieved by choosing about 6-8 circles: A minor number would give unreliable results, while a larger number could cause unnecessary long computational times (Figure 10b,d).
Concerning the block volume, this method relies on the assumption that the blocks are determined by the totally persistent discontinuities; therefore, the block volume could be underestimated. More precise results can be obtained directly from measurements on 3D point clouds, with respect to 2D analyses.
The script for the semi-automatic identification of discontinuity sets can be used in different zones of an orthophoto to identify potential deformation zones or changes of the stress field. Figure 19 depicts a shear deformation zone in sector B detected from two additional joint sets (with mean strikes N-S and E-W) and higher standard deviations of the strike of J1 and J2 with respect to sector C, which are the result of non-linear shear structures.
The presented approach helps to calculate the geometrical parameters of the discontinuity sets affecting a rock mass in a less time-consuming, more precise, and safer manner compared to the conventional geostructural and geomechanical surveys. However, a complete characterization of rock masses also requires information on roughness, wall strength, aperture, filling, and seepage of the discontinuities that cannot be estimated from 2D analyses. The determination of wall strength, nature of filling, and seepage requires direct measurements at the site, while only centimetric apertures could be measured on high-resolution remote sensing products [29]. In addition, recent advances in the literature proposed methods to determine the roughness/undulation of discontinuities from point clouds [29,35,95,96].
In this perspective, the MATLAB routine can contribute to identify the main discontinuity sets from discontinuity traces in low-relief rock masses characterized by discontinuities perpendicular to the strata and combined with point clouds, to localize accessible and representative surfaces in order to measure the non-geometrical properties by means of conventional techniques. Finally, the combination of geometrical and non-geometrical properties of discontinuities and point clouds or meshes can be used to create discrete models for quantitative stability analyses in rock masses by means of discontinuum approaches. It is specified that interactions among individual fractures or discontinuity sets (fracture topology) need to be defined to characterize rock mass permeability, especially when dealing with fault damage zones and reservoir modelling [41,[97][98][99][100][101][102][103]. In this perspective, the MATLAB routine developed could be improved through the characterization of fracture nodes and terminations [104,105] to define the network connectivity and its influence on the rock mass physics [106,107].

Conclusions
This study illustrated a new MATLAB tool for 2D semi-automatic analyses of discontinuities from high-resolution orthophotos obtained by means of remote sensing techniques, aimed at characterizing rock masses. Although similar tools were presented in the literature [29,40], they are not publicly available or do not provide a complete description of discontinuities for rock mass characterization. Therefore, the routine was compiled by adapting and updating the standard methods (i.e., scanline, rectangular, and circular window sampling) for a complete and detailed analysis of discontinuity traces, which is preferable to 3D point cloud analyses for low-relief rock masses or man-made excavations, in a user-friendly digital environment.
The code was initially built on a synthetic dataset and successively tested and validated on a case study by comparing the results of the conventional geostructural and geomechanical surveys carried out in the same area. The routine was developed in the form of consecutive steps, which can be singularly run depending on the objective of the analysis. In addition, the calculations do not require high-performance computers but can be run on standard laptops in a few seconds.
A new feature allows us to semi-automatically identify and classify the mean discontinuity sets in a Graphic User Interface, by fitting Gaussians curves on the strike histograms of the traces or by using the Hough transform, according to the number and approximate orientation of the discontinuity sets. In addition, the normal spacing of the discontinuity sets can be calculated both for persistent and non-persistent joints. The procedure for the classification and characterization of the discontinuity sets, as well as the estimation of the jointing degree of the analyzed area, can be easily repeated in different parts of an orthophoto to identify potential changes of the mean discontinuity sets, implicating modifications in the stress field, which can be further investigated on site for structural analyses.
Future developments may concern the improvement of the unit block volume calculation through direct measurements on the areas delimited by discontinuities rather than from spacing and orientation of the main discontinuity sets. Moreover, an automaticsemiautomatic method for drawing the discontinuities on orthophotos could exceptionally decrease the time required for data pre-processing.
Eventually, future research topics could deal with the conversion of the proposed Hough Transform method in the 3D space to detect discontinuity surfaces from point clouds or triangulated surfaces. The results could be compared with the methods available in the literature to validate the technique's reliability or identify potential discrepancies [108][109][110]. In addition, the Gaussian fitting method could be implemented to have 3D histograms of the orientations (x = dip, y = dip direction, z = frequency) and a composite Gaussian surface instead of a curve, to detect discontinuity planes from 3D models. Finally, discontinuity sets' spacings could be calculated both from the 3D Hough Transform Method and from the 3D Gaussian fitting method and compared with each other. Data Availability Statement: The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.