Due to the strong distortion of photos taken by UAVs, it is difficult to effectively apply the traditional extraction method based on geometric and texture features to extract feature points. The scale-invariant feature transform (SIFT) algorithm is characterized by the three properties of scaling, rotation, and affine invariance, via which it is capable of resisting certain illumination changes and viewpoint transformation. Hence, the SIFT algorithm is adopted for feature point extraction in the proposed method. The main concept of the SIFT algorithm is that the scale-space representation of UAV aerial images is established, after which the extreme points of images are searched in the scale space and extracted as feature points.

Image matching is conducted to reconstruct 3D information from multiple 2D images. However, the process of image matching using only SIFT feature points is slow. The location data of GPS coordinates in the images collected by the UAV and the attitude angle data provided by the inertial measurement unit (IMU) can assist in the construction of the topological structure between images. Next, the nearest-neighbor method is utilized to find the corresponding relationships between the feature points of images and establish a set of matching feature points that meet the geometric constraints. A large number of coordinate points constitutes a 3D point cloud of the target object in space.

The image points in the photo are projected into spatial coordinates according to the principle of camera imaging. The error function is defined as the sum of squares of the reprojection errors. The objective function is defined as follows: g ( C p , X ) = ∑ i = 1 n ∑ j = 1 m v i j f ( P ( C i , X j ) , q i j ) 2 (1) where C p = { C 1 , C 2 , C 3 , ... C n } are camera parameters, X = { X 1 , X 2 , X 3 , ... X m } are the coordinates of space points, the variable v i j represents the visibility of space point X i in camera C i , n is the number of images, m is the number of feature points obtained by precise matching, and the function f ( P ( C i , X j ) , q i j ) 2 represents the projection error of point X j in camera C i .

Finally, sparse beam adjustment is used for step-by-step iteration to minimize the reprojection error between the projected points and the points on the observed images, thereby calculating coordinates of the 3D point cloud in an optimal camera pose and camera scene.

A point cloud model commonly includes massive target points in a 3D region, and lacks topological information. Therefore, the principal problem of processing point cloud data is to establish a topological relationship among discrete points and realize the fast search of the nearest adjacent points.

A k-d tree (referred to as a k-dimensional tree) is a data structure that represents spatial partitions, and is mainly applied to search key data in multidimensional space (such as range searches and nearest-neighbor searches). In this study, the k-nearest neighbors are searched by the k-dimensional tree.

The problem of determining the normal of a point on a surface is similar to the problem of estimating a section of a normal of a surface. Therefore, the problem can be transformed into one concerned with least-squares plane-fitting estimation. In this study, the surface normal is evaluated by analyzing the eigenvectors and eigenvalues (or principal component analysis, PCA) (Riquelme et al., 2014; Robson et al., 2016; Guo et al., 2017) of the covariance matrix created from the nearest points. The covariance matrix C corresponding to each point P i can be defined as C = 1 k ∑ i = 1 k ( P i − P ¯ ) . ( P i − P ¯ ) T , C . V j → = λ j . V j → , j ∈ { 0,1,2 } , (2) where k is the number of the adjacent points of P i , P ¯ represents the 3D centroid of point sets, and λ j and V j → are the eigenvalue and eigenvector of the covariance matrix, respectively. The normal vector can be determined by the eigenvector that corresponds to its minimum eigenvalue. After all the normal vectors are obtained, the subsequent step is to extract rock mass discontinuities.

Given a point cloud set P = { P 1 , ... , P N } and the normal vector N = { n 1 , ... n N } of all points, the output result is a series of parameters ψ = { ψ 1 , ... ψ N } of the plane model. In this paper, local sampling is proposed to acquire new candidate planes in each iteration. The RANSAC algorithm is then applied to determine the parameters of the plane model with the highest score (i.e., the largest number of inliers). All the candidate planes are placed in the set C , and a new evaluation method is used to calculate the score m of the best plane. Moreover, | m | is the number of points in a candidate plane, | c | is the number of candidate planes, and p ( | m | , | c | ) is the probability of ignoring the planes with a higher score. When p ( | m | , | c | ) is large enough, the extracted plane is the best, and the remaining points will be used for the subsequent iteration. When p ( τ , | c | ) is sufficiently large, the iteration process is terminated. Finally, τ (default) is the number of minimum points on a plane.

Consider a point cloud with N points and a plane with n points, and k is the number of points in the minimum point cloud set that determines a plane. Provided that any subset with k points will generate a plane model, then the probability of detecting the plane model in one iteration is as follows. P ( n ) = ( n k ) / ( N k ) ≈ ( n N ) k (3)

When s candidate planes are detected, the probability of detecting the plane Ψ is as follows. P ( n , s ) = 1 − ( 1 − P ( n ) ) s (4)

The threshold value p t is artificially set. The number of planes T that meet the requirement P ( n , T ) ≥ p t can then be obtained by solving s, as follows. T ≥ ln ( 1 − p t ) ln ( 1 − P ( n ) ) (5)

Because the value of P ( n ) is commonly small, its logarithm ln ( P ( n ) ) can be expanded by its Taylor series: ln ( 1 − P ( n ) ) = − P ( n ) + O ( P ( n ) 2 ) . The substitution of ln ( P ( n ) ) into Eq. 5 yields the following. T ≈ − ln ( 1 − p t ) P ( n ) (6)

The complexity of an algorithm is closely related to the sampling method. The sampling method used in this study is detailed as follows. Shape is a local phenomenon, and the closer two points are, the more likely they belong to the same plane. The sampling efficiency can be greatly improved by utilizing this characteristic. Research has revealed that it may be effective to increase the number of inliers within the model by utilizing the locality of the shape to the sample (Myatt et al., 2002). Generally, in random sampling, a circle with a given radius is used to randomly select sample points, but the radius needs to be determined in advance according to the density and distribution of the points. However, the density and distribution of outliers vary greatly for different models; even at different locations on the same model, the density of outliers can vary dramatically. Therefore, a method is presented in this paper to adapt to the density of outliers.

The octree structure is an effective method by which to establish spatial proximity between sampling points. First, the point p 1 is chosen without restriction to create the candidate plane, and then a set C that includes p 1 is randomly selected from the structural layers of the octree structure. Finally, the remaining K − 1 sample points are selected from the set C . The probability of finding the plane ψ containing n points in this manner can be calculated as follows: P l o c a l ( n ) = P ( p 1 ∈ Ψ ) P ( p 2 ... p k ∈ C ) (7) where n / N is the first probability value, and the second probability value relies on the choice of the elements in the set C . The set C is considered to be an optimization if abundant points on the plane ψ are included in the   set. Most points on a plane, excluding boundary points and edge points, have neighbors that belong to the plane. Generally speaking, although the adjacent points on a plane cannot be determined by the element set of the octree structure, the candidate plane model includes a large number of points to ensure that it is more representative of the actual data. Therefore, the number of these neighbors must be as large as the number of elements in the octree structure, excluding a few points. To facilitate the analysis, it is assumed that the set C embodies all points p i on the plane ψ ( p i ⊂ ψ ) , the number of all points on the plane ψ is half of the number of points in the set C , and the other half of the points in the set C contains outliers or noise points. The probability of choosing a large set C is conservatively evaluated to be 1 / d , where d refers to the depth of the octree structure. Therefore, the conditional probability of choosing points p 2 and p 3   ( p 2 , p 3 ⊂ ψ ) from the set C set can be calculated by Eq. 8. Then, by introducing Eq. 8 into Eq. 7 and Eq. 9 can be obtained. ( | C | / 2 k − 1 ) ( | C | k − 1 ) ≈ ( 1 2 ) k − 1 (8) P l o c a l ( n ) = n N d 2 k − 1 (9)

The evaluation function σ p is utilized to evaluate the extracted candidate plane model. This evaluation function mainly includes the following three aspects. 1) After the selected plane is created, the points whose distance to the plane is less than the distance threshold value ε are regarded as the points on the candidate plane. 2) The points that meet the distance requirement will be further filtered. When the angles between the normal vectors of the detected points and the normal vectors of the candidate planes are less than the angle threshold value α , these points will be selected as the points in the plane. 3) A new threshold value β , which represents continuity, is added to the proposed method. For the points that have met the first two requirements, only those satisfying the continuity requirement can be selected as inliers on the plane.

In short, for a plane ψ , its evaluation function σ p can be expressed as follows. σ p ( Ψ ) = | P ψ | (10)

For example, points P ψ on a plane model ψ can be defined in two steps: P Ψ = { p | p ∈ P ∩ d ( Ψ , p ) < ε ∩ arccos ( | n ( p ) . n ( Ψ , p ) | ) < α } (11) P ψ = max c o m p o n e n t ( Ψ , P ψ ) (12) where d ( ψ , p ) is the Euclidean distance from point p to plane ψ , n ( p ) is the normal vector of point p , n ( ψ , p ) is the relationship between the normal vector of the plane model ψ and the projection of the normal vector of point p on the plane ψ ,  and  max c o m p o n e n t ( ψ , p ψ ) refers to the point set where the projected points on the plane ψ can form the largest connected part.

The discontinuity of a plane can be represented as follows. A x + B y + C z + D = 0 (13)

The normal vector of the plane is as follows. N = ( N x , N y , N z ) = ( A , B , C ) / A 2 + B 2 + C 2 (14)

The orientation of the discontinuity can be calculated by the following equations.When N z > 0 , β = cos − 1 ( N z ) (15) If N x ≥ 0 , α = cos − 1 N y N z 2 + N y 2 (16) If N x < 0 , α = 2 π − cos − 1 N y N z 2 + N y 2 (17) When N z < 0 , β = cos − 1 ( − N z ) (18) If − N x ≥ 0 , α = cos − 1 − N y N z 2 + N y 2 (19) If − N x < 0 , α = 2 π − cos − 1 N y N z 2 + N y 2 (20)