1. Introduction
Rock slopes are geological entities with unique structural characteristics, displaying significant discontinuities, heterogeneity, and anisotropy. Due to complex tectonic movements at different historical periods, surface weathering, and stress release processes, rock masses often develop various discontinuities of different scales and properties, such as fractures, joints, bedding planes, and faults. The distribution and evolution of discontinuities in 3D space govern the deformation and failure characteristics of rock masses. It is one of the important factors that affect the strength, permeability, and stability of rock masses [
1,
2]. The characterization of rock mass discontinuities involves identifying the location, size, and some typical parameters of the discontinuities. Several geometric parameters can quantitatively describe rock discontinuities, such as orientation, roughness, spacing, persistence, and more [
3]. Among these, the orientation (dip angle and dip direction) is particularly critical for analyzing the stability of jointed rock masses [
4]. For instance, orientation is a critical input parameter in block theory [
5], rock mass classification [
6], and the use of the limited equilibrium method [
7]. Hence, it is extremely important to accurately map these discontinuities and determine the orientation of rock discontinuities in rock engineering [
1,
8].
Traditionally, geologists use tools such as geological compasses and tape measures to manually measure the discontinuity plane of rock masses using methods like scanline surveys or window sampling to collect their geometric information [
9,
10]. Although this method is straightforward and user-friendly, it has several significant drawbacks: (1) The measurement process relies entirely on manual execution, which results in low efficiency [
11]; (2) The working conditions are easily constrained by terrain and environmental factors [
12]; (3) In complex terrains, especially in fragmented rock areas and steep slopes, surveyors may be unable to access specific areas for measurement [
13]; (4) Investigators need to have relevant professional backgrounds, and the back of these makes it easy for human errors to occur [
14,
15]. These limitations constrain the ability to conduct an in-depth analysis of the structural characteristics of rock slopes.
In modern geotechnical engineering and geological exploration, Light Detection and Ranging (LiDAR) and digital photogrammetry are widely used remote sensing techniques for acquiring high-quality 3D point cloud data, capable of capturing high-resolution 3D information [
16,
17]. The origins of 3D laser scanning technology can be traced back to the mid-20th century. This technology uses 3D laser scanners to rapidly scan and measure target objects, directly capturing the distance, direction, and intensity values of the reflected laser from the object. Due to its ability to collect real-time, high-precision 3D point clouds of rock surfaces without physical contact, this technology eliminates the subjective factors associated with manual contact measurements. Therefore, compared to traditional manual measurement methods, 3D laser scanning offers a more objective and accurate means of data acquisition. These techniques have seen rapid advancement in geological and geotechnical research fields [
18]. Rock mass discontinuities can be categorized into planar and linear types. Linear discontinuities primarily manifest as fractures, with their spatial distribution occurring within the rock mass. On the other hand, planar discontinuities are exposed fault surfaces that appear on the rock mass surface after sliding has occurred. This study focuses on the characterization of planar discontinuities. In recent years, many researchers have utilized 3D point cloud data to analyze discontinuities, including discontinuity recognition [
19,
20,
21,
22], slope stability analysis [
23,
24], discontinuity and rapid characterization [
25,
26]. The widespread use of remote sensing technology in the analysis of rock discontinuities has demonstrated its feasibility and effectiveness.
Currently, the main algorithms used in this field fall into three categories: model fitting, region-growing, and unsupervised clustering. The Random Sample Consensus (RANSAC) algorithm and Randomized Hough Transform (RHT) are two typical model-fitting methods that can detect planes in 3D point clouds. Ge et al. [
27] employ the region-growing method to determine the geometric characteristics, including orientation, spacing, and roughness, of rock discontinuities. They perform growth through the differences in normal vectors, thereby obtaining sets of discontinuities with different orientations. The performance of this method is highly dependent on the selection of initial seed points. Unsupervised clustering classifies point clouds based on the similarity of their underlying structure, geometric features, or spatial distribution. Unsupervised clustering methods typically surpass the model-fitting approach in extracting discontinuities from point cloud data, especially when dealing with complex geometries and a wide range of shapes and sizes [
22]. K-means [
28,
29] and its variants [
30,
31] are common distance-based clustering methods. These methods directly assign groups to each point, which can lead to inaccurate division of boundary points. Slob et al. [
32] proposed the automatic identification of discontinuities using fuzzy k-means clustering. This method introduced fuzzy theory and addressed the boundary grouping problem more effectively by calculating membership degrees. However, a significant drawback of such algorithms is the need to preset the number of clusters, making the clustering results highly influenced by subjective factors. To address this issue, scholars have proposed several algorithms that can automatically determine the number of clusters, including the surface classification algorithm [
33], improved K-means [
34], firefly algorithm [
35], unsupervised K-means [
36], and optimized K-means [
37]. However, partition-based clustering methods, like region-growing techniques, are affected by the choice of initial cluster centers, which affects the clustering results. Density-based clustering algorithms can avoid this issue. The Density-Based Spatial Clustering of Applications with Noise (DBSCAN) determines cluster boundaries based on the density around data points and is commonly used to further extract single discontinuity clusters from main discontinuity sets [
38]. Mean Shift (MS) is a window-based method that enables samples to converge at density peaks. Compared to DBSCAN, this algorithm handles high-dimensional feature data more effectively and requires only one parameter, the bandwidth. It is robust against outliers and is frequently used for identifying main discontinuity sets.
Daghigh et al. [
39] conducted a qualitative evaluation of RHT, region-growing, and RANSAC on a rock slope point cloud dataset. The results indicate that these methods exhibit poor performance when segmenting the discontinuity plane, often detecting spurious planes and over-segmenting. Due to the curvature and roughness of rock mass surfaces, point clouds often exhibit undulations [
34]. Some discontinuities exhibit significant local undulations or even fragmentation, which can lead to over-segmentation. Therefore, recognizing discontinuities is a complex process that typically requires the use of multiple algorithms in combination. Leng et al. [
40] combined Hough-Transform and Region-Growing (HT-RG) methods to segment the rock mass planes. This method can detect discontinuities of different sizes and has good segmentation accuracy. However, this approach may miss detecting small and low-density plane surfaces and is prone to over-segmentation. Hu et al. [
41] combined RANSAC with Voxel-Based Region-Growing (RAN-RG) for plane extraction from a point cloud. This method effectively improved the issue of detecting spurious planes, but still experienced over-segmentation in areas with significant local surface undulations. Yu et al. [
42] introduced the Supervoxel concept to improve RAN-RG (Supervoxel-RG), which to some extent improved over-segmentation. However, it sometimes results in holes within the segmented plane. To address the issue of over-segmentation, Daghigh et al. [
22] refined the surface variation of the rock point cloud by modifying the search neighborhood size and performing downsampling. However, some discontinuities with undulations or small offsets still remain. Chen et al. [
26] presented a novel approach that considers points with high curvature values as noise and proposed a modified MS algorithm that incorporates five denoising operations. While this method can improve over-segmentation to some extent, it inevitably overlooks some small, highly undulating discontinuities.
In summary, existing methods for characterizing rock discontinuities are not effective in improving the issue of over-segmentation, which impacts the accuracy of the results. To tackle this problem, we propose a method to improve over-segmentation and extract discontinuities along with their orientations from 3D point cloud data of rock surfaces.
2. Methods
In this section, we introduce a workflow that integrates unsupervised clustering with RANSAC model fitting, utilizing a watershed-like algorithm to effectively address the issue of over-segmentation in rock discontinuities. The methodology is structured into five main steps:
Step 1: Adaptive resampling of the point cloud data based on density ratio, transforming the dataset through adaptive scaling according to density ratio. This task is detailed in
Section 2.2.
Step 2: Estimating normal vector using Principal Component Analysis (PCA) and orienting them in the positive direction. This task is outlined in
Section 2.3.
Step 3: Determination of the main discontinuity sets, which includes: (1) Applying a watershed-like algorithm to identify the non-planar areas of the discontinuities. This part is described in
Section 2.4.2. (2) Segmenting the dominant discontinuity sets using the Multi-threshold Mean Shift (MTMS) algorithm. This part is covered in
Section 2.4.3.
Step 4: Recognizing individual discontinuity clusters within each set using the DBSCAN algorithm. This task is explained in
Section 2.5.
Step 5: Fitting a plane to each cluster from Step 4 using the RANSAC algorithm and calculating the orientation of each discontinuity. This task is described in
Section 2.6.
A summary of this procedure is illustrated in
Figure 1.
2.1. Description of the Experimental Datasets
The experimental dataset used in this study originates from a quarzitic rock mass at a road cut slope in Ouray, CO, USA (
Figure 2a). It was acquired from the Rock bench repository [
43] using an Optech Ilris 3D Laser Scanner. It consists of 1,515,722 points with an average point spacing of approximately 2 cm.
Figure 2b shows its corresponding trimmed point cloud, which contains 1,024,521 points. This dataset is a popular resource in the fields of geotechnical and geological engineering for estimating the orientations of structural surfaces. In previous studies, two sets of ground truth orientations have been documented. Numerous researchers have created both semi-automatic and fully automatic methodologies to determine the orientations of structural surfaces within this point cloud.
2.2. Density-Ratio-Based Resampling Method
Point cloud data obtained from laser scanning often suffer from issues like uneven data density distribution due to factors such as vegetation cover, occlusion, and noise. Density-based clustering algorithms, such as DBSCAN, are highly sensitive to variations in point cloud density and struggle to detect all clusters in datasets with significantly varying densities. Therefore, in our method, we draw inspiration from an image processing technique known as histogram equalization. This technique applies nonlinear stretching to an image, redistributing the grayscale values so that the grayscale levels within a certain range become approximately equal. The most common method for equalization is the Cumulative Distribution Function (CDF).
Based on this principle, we substitute the grayscale values with the local point density of the point cloud and use a density estimation algorithm as the density estimator to estimate the local density of each sampled point. We then resample the data according to the CDF estimated over a neighborhood (
). The effect of this process is to achieve a roughly uniform spatial distribution of the data [
44]. This step enhances the accuracy of the DBSCAN algorithm in recognizing single discontinuity clusters, as described in
Section 2.5. The key symbols and notations utilized in this study are provided in
Table 1.
The specific process is as follows.
- (1)
The original dataset is . When it is divided into equal intervals, the boundary positions are ;
- (2)
At each boundary position , calculate ;
- (3)
Calculate the resampled dataset: .
The theoretical basis for steps (2) and (3) is as follows.
Let the CDF at the sampling point
be
. We define the cumulative distribution function. We resample the dataset
into
using the follow mapping function:
This is a probability integral transform, where , representing a uniform distribution function.
Specifically, let
be the average density of
within a neighborhood of size
, then its density in the complete space is expressed as:
When
is large and varies slowly around
, this expression can be approximated as:
Therefore,
, i.e., when resampling with a smaller
, the resampling result is scaled by different degrees according to the density ratio between the small neighborhood
and the large neighborhood around the sampling point. In other words, when the sampling point
is in a low-density region, the points in that area are less thinned out, while in a high-density region, the points are thinned out more significantly. This process is summarized in
Table 2.
2.3. Normal Vector Calculation
Accurate normal estimation in a point cloud is essential for effectively segmenting the desired plane features. Principal Component Analysis (PCA) is a statistical method that focuses on variation to uncover significant patterns in a dataset through an orthogonal transformation. PCA is widely used to identify local planar surfaces and linear points within a 3D point cloud. Using KD-tree nearest neighbor search, extract the neighborhood point sets (
) for each sample point within the search radius (
), respectively. Calculate the centroids of each neighborhood point set (
). Construct the covariance matrices (
) relative to the centroids, which reflect the correlation of each coordinate component in the point sets. The 3×3 covariance matrix for a sample point
can be written as:
Utilize Singular Value Decomposition (SVD) on the covariance matrix of each point. The covariance matrix
can be represented as:
Among them, are the eigenvalues of , and , , are the corresponding eigenvectors. After normalizing the eigenvalues and sorting them by size, we have , where , and represent the fitting deviations in the directions of the corresponding eigenvectors of the local plane. The eigenvector corresponding to the smallest eigenvalue, , indicates the direction of the least variation in the point set distribution. This direction can be considered the normal direction of the local surface.
Since PCA can only provide a direction without distinguishing between “inward” and “outward”, an additional post-processing step is required to correct any normal vectors that are flipped by 180°. This study uses a method based on the continuity of normal directions to adjust and ensure the consistency of the normal vector orientations.
2.4. Determination of the Main Discontinuity Sets
Identifying group discontinuities is fundamental for pinpointing individual discontinuities, which is essential for accurately identifying rock discontinuities. This phase focuses on categorizing each point into a discontinuity group based on distribution patterns.
MS is a feature space analysis method that estimates the probability density distribution of feature data through Kernel Density Estimation (KDE). It clusters data points by moving them in the direction of increasing probability density until they converge at the local density maximum. The method can process point clouds of any shape without needing a predefined number of clusters. It requires only one parameter, the bandwidth size.
The bandwidth determines the step size for each adjustment of the probability density function.
For MS, tuning the bandwidth is the most critical hyperparameter.
2.4.1. The Impact of Bandwidth on MS
The bandwidth, as the sole input parameter for MS, has a crucial impact on the segmentation effect.
Figure 3a illustrates the kernel density of a one-dimensional dataset. The dashed lines indicate the kernel function for each sample, while the solid line depicts the probability density of the dataset, which is the aggregate of the sample kernel functions. It can be observed that the probability density of this dataset has only one local maximum, so the current dataset is segmented into one cluster using MS with a bandwidth of 6. If we set the bandwidth to approach zero infinitely, the sample kernel functions will be as shown in
Figure 3b, approximating impulse functions, and the probability density will be as shown in
Figure 3c. The dataset will have five local maxima and be segmented into five clusters. Therefore, the bandwidth determines the probability density of the dataset. If the bandwidth is too large, MS may overly smooth the data, failing to capture local details and features, resulting in under-segmentation. Conversely, if the bandwidth is too small, it will include fewer points, leading to the amplification of some local details and features, causing over-segmentation.
Figure 4 illustrates the effect of bandwidth on different types of discontinuities.
Figure 4a shows a representative area of the experimental dataset with flat discontinuities. The spatial orientation of plane 1, plane 2, and plane 3 is illustrated in
Figure 4b, where each plane belongs to a different set of discontinuities. By comparing segmentation effects with different bandwidths in this area, we find that the segmentation effect in
Figure 4c with a bandwidth of 0.25 is ideal because it effectively segments the area into three discontinuities while preserving edge details. For the non-planar discontinuities targeted in this study (
Figure 4d),
Figure 4e simulates their planes in 3D space. It is clear that the non-planar discontinuities fit well with the plane model, suggesting that most points on the non-planar discontinuities should theoretically be classified into the same discontinuities set.
Figure 4f shows the segmentation effect of using different bandwidths for MS in this area. When the bandwidth is 0.25, the points are over-segmented. To avoid over-segmentation, the bandwidth is increased, and the result obtained with a bandwidth of 0.35 is considered the most suitable representation because most points in this area are assigned to the same discontinuities set. However, for flat discontinuities, when the bandwidth is 0.35 (see
Figure 4c), plane1 and plane3 are assigned to the same discontinuities set, resulting in under-segmentation.
Based on the comprehensive analysis above, it is evident that the optimal bandwidth for non-planar discontinuities and flat discontinuities differs. However, most existing MS methods use a global fixed bandwidth [
22,
26]. To ensure accuracy, this often leads to over-segmentation in areas with non-planar discontinuities. On the other hand, expanding the threshold range can reduce accuracy when segmenting flat discontinuities. This makes it difficult to balance both types of discontinuities within the same dataset, leading to over-segmentation or under-segmentation. Therefore, it is necessary to identify non-planar areas and set different bandwidths for points located in these areas and those in flat discontinuities. To tackle this issue, we propose a new method for detecting individual discontinuities, which primarily includes: (1) applying a watershed-like algorithm to identify non-planar areas of discontinuity. (2) using the Multi-threshold Mean Shift (MTMS) algorithm to segment the advantageous discontinuity sets. Detailed descriptions of these improvements are provided in the following subsections.
2.4.2. Watershed-like Algorithm for Identifying Non-Planar Areas
The Watershed algorithm is an image segmentation technique that separates and labels different regions within an image. It is based on the grayscale values or gradient information of the image, treating the depth image as a topographic map [
45]. The algorithm assigns a gradient value to each pixel by calculating gradient operators, and simulating peaks and valleys (basins). As illustrated in
Figure 5, assuming there are small holes on the surface of the minimal regions of the terrain, water seeps in through these holes and gradually floods the surrounding areas at a uniform rate, forming corresponding catchment basins. The boundary line of each catchment basin is known as the watershed line [
46,
47,
48].
Figure 6 demonstrates the Watershed algorithm’s segmentation process in an image [
49].
Figure 6a shows a limestone-blasted grayscale image, treated as a gradient image, where the grayscale values at corresponding coordinates represent elevation.
Figure 6b depicts the 3D simulation of the water catchment basin in this image, assuming a gray threshold of 50 as the watershed line.
Figure 6c shows the 3D simulation post-flooding, and
Figure 6d presents the segmented 2D cross-section.
Current researchers are exploring the application of watershed-like concepts or modified algorithms to point cloud data, as the density or other attributes of point clouds may, to some extent, resemble “height” or “intensity”. By simulating this analogy, “terrain” in 3D space can be defined. Based on this idea, this paper proposes a watershed-like algorithm that uses non-planar indicators as gradient operators to restrict the flooding process of the watershed algorithm, thereby better identifying non-planar regions of discontinuities.
Calculate the Non-Planar Indicator
In
Section 2.3, we used PCA to compute the local covariance matrix eigenvalues
,
,
as shown in Equation (2). According to the theory of Demantké et al., as shown in
Figure 7, the eigenvalues can represent the dimensional properties of the local neighborhood within the point cloud [
50]. When
, the local neighborhood of the point cloud tends to be rod-shaped; when
, the local neighborhood of the point cloud tends to be planar; when
, the local neighborhood of the point cloud tends to be a three-dimensional surface.
Therefore, we construct a local neighborhood dimension feature model based on these dimensional characteristics [
51]:
In the study by Demantké et al. [
50], the dimensional features
,
and
can, respectively, represent the probability that the local neighborhood of the point is a one-dimensional line, a two-dimensional plane, or a three-dimensional surface [
50]. The planar dimension feature
indicates how well the local neighborhood of the point cloud aligns with a plane. The closer the point is to a planar structure, the larger its corresponding dimension feature
. Therefore, this paper defines an indicator to describe the curvature of the local neighborhood:
From Equation (7), it can be observed that the smaller the value of , the closer the local neighborhood shape is to a plane. When the local neighborhood is a curved surface, the value corresponding to is a fixed value related to the curvature of the surface.
To distinguish between non-planar and flat areas, as shown in
Figure 8, surfaces with different degrees of curvature have varying normal direction deviations within the same neighborhood range. The more curved the surface, the greater the deviation in the normal direction. Therefore, we use the angle between the normal vectors of the sampling point and its neighboring points as a reference factor.
To simplify the complex curvature calculation process, we establish a local neighborhood
, centered on the sampling point. The size of
is customized according to the scale of the rock mass. Given that the number of neighborhood points within
is
, the curvature at the sampling point is defined as follows:
where
represents the deviation between the sampling point and the normal vectors of the neighboring points within the local neighborhood. It can be observed that as
increases, indicating a greater degree of bending in the local structure, the curvature value correspondingly increases.
Combining the two-dimensional feature
, this paper defines a non-planar indicator:
where
is an indicator that measures the degree of bending within the local neighborhood, and
describes the curvature of the local neighborhood. The shape characteristics of the point cloud’s local neighborhood are determined based on the non-planar indicator
. When the local neighborhood is close to being planar, the non-planar indicator approaches 0. As
increases, it indicates a greater degree of bending in the local neighborhood of the point cloud.
Fitting Spatial Gradient Fields
Due to the similarity in curvature among closely spaced points, we sample the point cloud data to reduce computational load and improve efficiency. A voxelized grid method is used to create a voxel grid, which consists of 3D boxes in space, to encompass the input point cloud. The voxel size is determined by the user based on specific requirements. To sample feature points, select point
located at the centroid of each voxel, and calculate the non-planar indicator
for each feature point in the experimental dataset.
Figure 9 is an HSV simulation diagram of the experimental dataset colored according to the non-planar indicator
. It can be observed that the values of
in this dataset primarily range
, The smaller the
value, the flatter the region where the point is located; the larger the
value, the more curved the point’s location is. We use the ascending
values to simulate the elevation from valleys (basins) to peaks in the high-dimensional features, with the peaks representing the non-planar areas we aim to identify.
The watershed algorithm typically processes gradient data [
52]. In this study, the Sobel operator is used to calculate the gradient data of the voxels
. Assuming
,
,
,
,
,
,
,
,
,
,
,
,
,
represent the points in the 14-neighborhood of voxel
in the directions: top-left, top-right, bottom-left, bottom-right, front-top, back-top, front-bottom, back-bottom, top, bottom, left, right, front, and back, respectively. The gradient values of voxel V in the x, y and z directions, denoted as
,
and
, are given by:
The gradient magnitude
is:
For each point
, in the point cloud data, assuming the index of its voxel is
, we approximate the gradient value
of the point
as being equal to the gradient value of the voxel it resides in, that is
Since point cloud data exist in 3D space, it is not feasible to simulate ridges and valleys on the same horizontal plane as one would with images. To facilitate the simulation of the watershed inundation process, we use the normal direction of each point as the vertical direction of the peaks. Through this approach, we define the “terrain” in the multidimensional space, allowing for more convenient simulation and analysis.
To determine the frequency distribution of each gradient level, we need to utilize a sorting algorithm to arrange the experimental dataset by gradient magnitude. Each point is then assigned to an ordered queue corresponding to its gradient level. Assume the number of gradient levels is
:
where
is the maximum gradient value,
is the minimum gradient value,
is the interval of the gradient levels, and INT denotes the rounding operation.
The gradient level
at point
is shown as
From Equation (13), it can be seen that the gradient level interval determines the number of gradient levels N, which in turn determines the number of points submerged during each step of the flooding process. According to Equation (14), the gradient level for each point can be calculated. The flooding process submerges points sequentially according to their gradient levels.
Flooding
First, let us clarify some basic concepts in this method:
Minimum Region: The minimum region is an independent connected area where the gradient values within the region are lower than the gradient values of the surrounding neighboring points.
Water Catchment Basin: The water catchment basin is an area composed of a collection of points within the experimental dataset that form the minimum region . Within the water catchment basin, water flow from any location will slide down along a certain path to reach the minimum region .
Geodesic Distance: Assuming is a connected region, the geodesic distance is along the geodesic path from to within .
Geodesic Influence Region: Let set contain set , where consists of multiple connected regions . Within set , if there exists a point set whose geodesic distance to connected region is shorter than its distance to any other connected region within set , then is considered the geodesic influence region of connected region .
The flooding process is essentially a recursive process. Suppose the experimental dataset is
, and
is the threshold set for point
at level
,
The flooding process starts from the threshold set
, which is the set of points with the smallest gradients. Let the collection of first-level basins be:
Let n be the gradient level, and the recursive process of the
-th level basin
is
The set of catchment basins with
as the watershed line is
For the experimental dataset
, the set of non-planar areas
(slope peaks) is
According to Equation (18), the watershed line
determines the range of slope peaks and valleys. For the experimental dataset, we set the number of gradient levels to
,
Figure 10 shows a comparison of the flooding results for different watershed lines, where the red area represents the non-planar area (slope peaks). It can be observed in
Figure 10a–f that as the value of
decreases, the range of the identified non-planar area gradually increases. However, in
Figure 10g,h, even though
continues to decrease, the range of the identified non-planar area remains almost unchanged. The study found that there exists a critical value for the watershed line, which distinguishes between planar areas and non-planar areas—essentially, the cliff edge of the peak. Due to the significant gradient difference between planar areas and non-planar areas, the identification results remain nearly constant when
continues to decrease. In the experimental dataset,
represents the critical value of the watershed line, where the division at the peak of the slope accurately reflects the location of non-planar area. Therefore, we use
as the watershed line for the experimental dataset.
2.4.3. The Multi-Threshold Mean Shift (MTMS) Algorithm
Currently, most discontinuity recognition methods based on MS employ a fixed bandwidth. However, using a single global bandwidth for rock slopes with numerous non-planar areas is not optimal, as a single global value may result in over-segmentation or under-segmentation.
Section 2.4.1 provides a detailed explanation of the impact of bandwidth on segmentation performance. MTMS is an improvement over the traditional MS method. This approach assigns different bandwidths to points in planar and non-planar areas during the shifting process. By doing so, it controls the step size of the probability density function adjustment for points in different regions. This method achieves a better balance between local and global features. The specific process is shown in
Figure 11.
Map the estimated normal vectors from
Section 2.3 to the feature space, extending from a 2D Gaussian surface to a 3D Gaussian sphere (
Figure 11a), as the input for MTMS. Define a flag value to mark whether a point is in a non-planar area. Search for the points located in the non-planar area defined in
Section 2.4.2 in spherical coordinates, setting their flag to 1, while setting the flag to 0 for the remaining points. The resulting binary map is shown in
Figure 11b.
Based on the flag, construct a trigger-based bandwidth allocation mechanism as detailed in
Figure 11c. Randomly select a point from all unprocessed points as the center point, check its flag label, and if flag = 1, set its bandwidth to
; if flag = 0, set its bandwidth to
. The kernel function of MS assigns certain weights to points around each point. In the feature space, it weights the points within a specified bandwidth region around the current center point to estimate the local density. The density calculation formula is
where
is the kernel function weight formula,
is the number of features,
is the bandwidth of the current search point,
is a normalization constant that ensures
= 1.
Calculate the gradient direction of the dataset density, which is the direction of the fastest increase in local density, to obtain the drift vector. The gradient of the density function is:
where
is a real number, so the drift vector is
Translate the sample point x iteratively according to the drift vector. The center point for the next iteration is
Repeat the iteration several times until , so that the samples converge to the local density peak, resulting in multiple sets of dominant points with similar principal directions.
2.5. Recognition of Single Discontinuity Cluster
To distinguish each individual planar surface among many planes with similar orientations in the primary cluster, secondary clustering is required. Within each discontinuity set, each single discontinuity creates a dense region in space, resulting in a density difference between the discontinuities. This characteristic makes them ideal for separation using density-based methods [
26]. DBSCAN is a typical density-based clustering algorithm. It has performed well in previous studies for identifying individual discontinuities [
21,
22,
26,
53]. In this study, the DBSCAN algorithm is utilized, with parameters configured for the neighborhood search radius (
) and the minimum number of points needed to form a cluster (
). The algorithm starts with a randomly chosen starting point
and searches within a radius of
. If there are at least
points within this radius, the point (
) is classified as a core point. If the number of points is fewer than
, the point is considered a boundary point. Core points, along with all points within the radius
, are grouped into a cluster. This process is repeated for each unprocessed point, ultimately yielding the second-stage clusters that belong to each structural plane. the results of identifying the first set of discontinuities are illustrated in
Figure 12.
2.6. Discontinuity Plane Fitting
Given that rock discontinuity surfaces are typically rough and uneven, the orientation of a discontinuity is determined by the orientation of its best-fitting plane. RANSAC is an iterative method designed to estimate surface model parameters from a dataset. It identifies the parameters of a geometric shape by utilizing the smallest number of points needed to define that shape within a point cloud [
54,
55,
56].
It accounts for the variability of the points, offering a more accurate estimation. The RANSAC algorithm is widely applied in shape detection [
57,
58] and has been effectively used in previous research for plane fitting [
34,
59].
The RANSAC algorithm repeatedly applies a plane equation to fit the discontinuity, aiming to find the best parameters. The equation for the discontinuity plane is expressed as:
In this equation, A, B, C, and D are the coefficients, with A, B, and C representing the three components of the plane’s normal vector. D denotes the perpendicular distance from the origin to the plane. The dip angle and dip direction of the planes are determined using Equations (25) and (26).
2.7. Qualitative Segmentation Results
Daghigh et al. [
39] reviewed the experimental findings of reference [
60] and reference [
25], highlighting various examples of deficiencies associated with rectangles and identifying the type of deficiency (OS: over-segmentation; US: under-segmentation) (
Figure 13a,b). To optimize segmentation results, Daghigh et al. [
22] modified the KNN value and utilized VG downsampling to reduce the surface variation in the point cloud. Nevertheless, certain discontinuities with undulations or minor offsets remain (as shown in
Figure 13c, highlighted by the circled and rectangular areas). The visual segmentation results from these studies, along with the current study, are displayed in
Figure 13. It can be seen that the results of the present study (
Figure 13d) show an ideal segmentation effect in the areas marked with deficiencies, demonstrating that our method for determining the main discontinuity sets can effectively avoid over-segmentation.
For the extraction results of discontinuity orientations, the experimental dataset used in this study has several sets of true orientation values that have been published by other researchers. Riquelme et al. [
53] and Chen et al. [
60] manually measured the orientations of several typical discontinuities in this dataset. Daghigh et al. [
39] conducted a critical evaluation of these measurements and identified significant errors in the data obtained by Riquelme et al. [
53] and Chen et al. [
60]. As a result, Daghigh et al. [
39] re-evaluated the discontinuity orientations by employing the Segment tool in CloudCompare. This method entails drawing a closed polygon around the points related to a particular discontinuity plane and then using the best-fit plane tool to ascertain the plane’s orientation based on the chosen points. This paper uses these results as the true values for discontinuity orientations.
Figure 14 shows the numbering of the discontinuities used for comparison, and
Table 3 compares the measurement results of our method with those of reference [
22] and the true values. It can be observed that our method has smaller errors compared to reference [
22], validating that our method provides satisfactory characterization results.
4. Conclusions
This paper primarily investigates the problem of discontinuity characterization in unstructured point clouds of rock masses. To address the issue of over-segmentation that often arises when characterizing complex rock mass discontinuous surfaces, a new method based on laser point cloud data is proposed. The method is divided into five main stages: (1) adaptive resampling; (2) normal vector calculation; (3) determination of the main discontinuity sets; (4) recognition of single discontinuity cluster; and (5) discontinuity plane fitting.
The primary novel contributions of the proposed work are: (1) adaptive resampling of point cloud data based on density ratios, which has the capacity to process datasets with varying densities; (2) the watershed-like algorithm can identify the non-planar areas on the rock mass surface; (3) proposed MTMS, which includes a trigger-based bandwidth allocation mechanism, effectively avoiding over-segmentation in non-planar areas (4) accurate calculation of the orientation of discontinuities; (5) fewer critical hyperparameters require adjustment compared to other methods.
To verify the reliability of the results, we applied this method to three real rock mass datasets and compared it with previous research results and manual measurements. For the experimental dataset, it can be seen that the results of the present study achieve satisfactory segmentation in the set of marked deficiency areas, demonstrating that our method for determining the main discontinuity sets effectively avoids over-segmentation. When comparing the extracted discontinuity orientations from this dataset with the ground truth published by other researchers, our method shows smaller errors. For Case A, the identification results from our method align closely with the color-coded ground truth point cloud. Moreover, in areas previously recognized by scholars as prone to over-segmentation, our method provides better segmentation results, proving its high accuracy and effectiveness in mitigating over-segmentation. Furthermore, we validated our method on Case B by comparing the computed discontinuity orientations with measurements from CloudCompare. The average errors for the dip direction and dip angle across eight sample surfaces were all within 2°, indicating good generalizability of our method. We also analyzed the computational efficiency of our method by comparing its computation time with several representative methods on the same dataset. The results demonstrate that our method has a satisfactory computational efficiency, meeting the practical needs of engineering applications, especially for complex rock masses with many non-planar areas.