An MFF-SLIC Hybrid Superpixel Segmentation Method with Multi-Source RS Data for Rock Surface Extraction

Abstract

Multi-Source RS data integration is a crucial technology for rock surface extraction in geology. Both Terrestrial laser scanning (TLS) and Photogrammetry are primary non-contact active measurement techniques. In order to extract comprehensive and accurate rock surface information by the integration of TLS point cloud and digital images, the segmentation based on the integrated results generated by registration is the crux. This paper presents a Multi-Features Fusion for Simple Linear Iterative Clustering (MFF-SLIC) hybrid superpixel segmentation algorithm to extract the rock surface accurately. The MFF-SLIC algorithm mainly includes three contents: (1) Mapping relationship construction for TLS point cloud and digital images; (2) Distance measure model establishment with multi-features for initial superpixel segmentation; (3) Hierarchical and optimized clustering for superpixels. The proposed method was verified with the columnar basalt data, which is acquired in Guabushan Geopark in China. The results demonstrate that the segmentation method could be used for rock surface extraction with high precision and efficiency, the result of which would be prepared for further geological statistics and analysis.

1. Introduction

Rock masses behavior is governed by properties of rock surface including fractures, joints, faults and other geological structures. For a rock engineering research, it is necessary to make a detailed investigation for geological conditions. With the distribution, combination condition and geometrical information of rock mass, both engineering geological evaluation and rock mass classification would be proceeded well. Therefore, it is worthy for transportation engineering, hydropower engineering and mining engineering to extract rock surface accurately, efficiently and comprehensively, which is necessary for engineering design, exploration, construction and evaluation.

However, for those dangerous and inaccessible regions, such as alpine and gorge, the traditional contact measurements could not be proceeded efficiently and safely with a limited time. Terrestrial Laser Scanning (TLS) and close-range photogrammetry, which are the mainly two non-contact active measurements, could collect point cloud or digital images data and achieve rock surface extraction more comprehensively and conveniently in a virtual digital environment [1,2,3,4].

TLS could obtain an accurate 3D model of an object in high efficiency, and get the spatial coordinates of an object with high accuracy and speed, which could be used for several purposes [5]. With the advantages of high spatial resolution, high temporal resolution and uniform accuracy, the TLS data is more appropriate for many applications, especially for inaccessible and dangerous regions [3,6]. However, as the TLS data is characterized as unstructured 3D point clouds, it is poor at expressing the object features. As another non-contact active measurement technique, photogrammetry is used to process stereoscopic images delivering RGB and spectral channel information, which could be considered largely complementary to TLS [7,8]. So, the integration of the multi-source RS data has complementary advantages, which is valuable for geological investigation.

Registration is one crucial problem that the research of rock surface extraction based on the integration of the multi-source RS data involves, as the correlation of the TLS point cloud and digital images could be obtained by registration. About the registration of the two data for geological objects, the authors have provided a detailed method in another article [9]. For the integration of multi-source RS data, it is also necessary to comprehensively analyze and evaluate multi-feature of the data and make a correct and efficient segmentation, which is the key of accurate, efficient and full extraction for rock surface. So, this paper mainly discusses the segmentation problem which the TLS point cloud and digital images integration is involved in.

The first step to be finished is the construction of the direct mapping relationship between the two data sources. Texture mapping is a simple way for 3D visualization of the textured point cloud. However, the true mapping relationship for the point cloud and digital images could not be obtained. Therefore, exploring a new method for direct mapping relationship construction is the precondition of rock surface extraction. Furthermore, in order to get an efficient and accurate segmentation result, the difference as well as the correlation between the characteristics of geological objects in the Multi-source RS data should be taken full use of.

At present, there is only so much research on point cloud or image segmentation, which would be valuable for the segmentation of multi-source RS data integration. The point cloud segmentation methods include edge-based segmentation [10,11,12], model fitting-based segmentation [13,14,15], region-based segmentation [16,17,18,19,20] and feature clustering-based segmentation [21,22]. For the point cloud, the only spatial information is insufficient for segmentation. Segmentation based on the depth image generated from point cloud would result in an expensive computational cost and information loss, which may affect the precision of segmentation result. For the image segmentation, there are mainly three groups: threshold-based segmentation [23,24,25], edge-based segmentation and region-based segmentation [26,27]. Instead of the above segmentation methods in pixels, some scholars put forward the idea of superpixel segmentation, which is a preprocessing stage to group pixels into “superpixels” [28]. The superpixel is composed of a set of pixels with some similar features, including texture, brightness and color. The image feature description with superpixels could decrease the edge extraction difficulty of complex objects, and improve the efficiency of image segmentation. Generally, superpixel segmentation algorithms could be classified into: (1) graph-based, including Normalized cut, Superpixel Lattice and Entropy Rate; (2) gradient-ascent-based, including Watersheds, Mean-shift, Quick-shift, Turbopixels and SLIC (Simple Linear Iterative Clustering). All the above algorithms show different performance in efficiency, precision, boundary adherence and expansibility [29]. In comparison of segmentation algorithms in pixels, superpixel segmentation has a better consideration of the association between local regional features, merging and clustering based on which could improve the efficiency and precision of segmentation. However, as most of existing algorithms are proposed for image data, which only applies the 2D image features, the precision could not be guaranteed for those objects with the similar image features but different spatial or other features.

Based on the above analysis, this paper develops a new hybrid superpixel segmentation algorithm, MFF-SLIC, for accurate, efficient and full extraction of rock surface. This algorithm firstly obtained the mapping relationship between the TLS point cloud and digital images according to the registration result; and then explored a new distance measure model for initial superpixel segmentation. The new distance measure model integrated multi-features of both two types’ data, including color and pixel location of digital images as well as geological properties and spatial position of TLS point cloud. With the above distance measure model, the geological object could be divided into a number of superpixels which have certain semantic features; next, superpixels would be clustered hierarchically and optimally (i.e., the initial segmentation result would be further processed) to provide credible data for fine extraction of rock surface. At last, with the columnar basalt data (acquired in Guabushan Geopark in Nanjing, China), the segmentation method proposed in this paper was verified. The results show that the MFF-SLIC hybrid superpixel segmentation algorithm could achieve rock surface extraction accurately based on the integration of the Multi-source RS data, even without any manual interactions, which provides a valuable data basis for further geological statistics and analysis.

The rest of this paper is organized as follows: Section 2 reviews the registration overview about multi-source RS data and the existing 2D SLIC algorithm for image segmentation. Section 3 describes the detailed methodology of the MFF-SLIC hybrid superpixel segmentation algorithm for rock surface extraction. In Section 4, an experiment with the columnar basalt data is conducted in order to illustrate the feasibility and applicability in geology of the algorithm. Section 5 draws the conclusion.

2. Related Work

2.1. Registration for the Multi-Source RS Data

Segmentation needs to take use of the direct correlation between the TLS point cloud and digital images, which could be obtained according to the registration result. Existing registration algorithms for point cloud and images mainly includes four groups: (1) Registration by use of artificial targets (e.g., retro-reflective targets); (2) Registration through corresponding features respectively extracted from images and point cloud [30,31]; (3) Registration with the point cloud respectively generated by digital images and acquired by TLS [32,33,34]; (4) Registration with intensity image and digital image generated by point cloud [35,36,37]. For geological objects, as the indeterminate distribution and irregular appearance of rock mass, it is difficult to directly extract corresponding primitives including points, regular line segments and planar features from the point cloud and images.

A G-Super4PCS algorithm for registration of TLS point cloud and digital images is provided [9]. With the Structure from Motion (SfM) technique, a point cloud could be generated from digital images, which has different scale, position and attitude with the TLS point cloud. The G-Super4PCS algorithm firstly performs key-scale rough estimation combining spin images and cumulative contribute rate, then defines a generalized super 4-points congruent base as the registration primitive, and introduces the rock surface features contraints to improve the efficiency and accuracy of its extraction. On the above basis, two key point clouds respectively extracted from the TLS point cloud and generated from digital images are used for rough estimation of rigid transformation parameters, and then the original TLS point cloud and the dense SfM point cloud from digital images are used for fine registration.

2.2. D SLIC Algorithm for Image Segmentation

The SLIC is a gradient-ascend-based segmentation algorithm, which generates superpixels by clustering pixels with the distance measure according to their proximity and color similarity in the image plane [38]. An example is shown in Figure 1. Figure 1a shows the origin image, and Figure 1b shows the corresponding superpixel segmentation result.

The distance measure is constructed by a 5D vector including the color vector (L, a, b) and the pixel position coordinates (x, y). Instead of the Euclidean distance, the 2D SLIC algorithm introduced a normalized distance measure considering superpixel size. Suppose the size of an image I is W by H, then the number of pixels N is described as Equation (1).

$$N=W\times H,$$

(1)

This algorithm specifies K superpixels with approximately equal size as input. For image I with N pixels, the size of each superpixel is $N/K$. Therefore, the grid interval S of superpixels is computed as Equation (2).

$$S\approx \sqrt{N/K},$$

(2)

where S denotes the size of the superpixel. For the distance measure normalization, the value of S is used to determine weight coefficients of different feature components.

Choose K superpixel cluster centers $C_k={\left[L_k,a_k,b_k,x_k,y_k\right]}^T$ (k = [1, K]) at regular grid intervals S. The computational formula of distance measure is shown in Equation (3).

$$\begin{array}{ccc}\hfill D_S& =& d_{Lab}+m{d}_{xy}/S\hfill \\ \hfill d_{Lab}& =& \sqrt{{(L_k-L_i)}^2+{(a_k-a_i)}^2+{(b_k-b_i)}^2}\hfill \\ \hfill d_{xy}& =& \sqrt{{(x_k-x_i)}^2+{(y_k-y_i)}^2}\hfill \end{array},$$

(3)

where $D_S$ denotes the normalized distance measure; $d_{Lab}$ denotes the Euclidean distance of the color vector $\left(L_k, a_k,b_k\right)$ and $\left(L_i, a_i,b_i\right)$ in CIELAB color space; $d_{xy}$ denotes the Euclidean distance of pixel $\left(x_k, y_k\right)$ and pixel $\left(x_i, y_i\right)$ in the image plane; m is a variable to control the compactness of a superpixel. The greater the value of m, the more important spatial proximity and the more compact the cluster. Experience indicates that the optimal range of values for m is [1,20], which can get a good tradeoff between color similarity and spatial proximity [38].

2.3. Segmentation for the Multi-Source RS Data

Existing segmentation methods for image and point cloud data mainly concentrate on planar (or other geometrically primitive) feature extraction, most of which are computationally expensive and often utilize 2D or 3D information alone. Based on the above consideration, several segmentation schemes combining 2D features from digital image with 3D features from point cloud data are proposed.

Reference [39] builds 3D models of urban environments based on 3D range scans and 2D images. This method utilizes a range segmentation algorithm to extract planar regions and linear features so that a 3D range scan could be converted to a set of bounded planes and finite lines. However, the applicability of the algorithm is limited due to planar regions and linear features are only common for urban environments, but geological objects are often characterized as irregular features. Reference [40], for segmenting a colored laser point cloud, an extension method for camera-only graph-based segmentation is proposed. In this method, segment unions based on spatial proximity is proposed, and a dynamic segment union criterion based on color and surface normal produces a quality segmentation. However, this method requires a fixed positional relationship between camera and laser scanner in order for the co-registration. Reference [41] takes both colorimetric and geometric data as input, clusters pixels and divides the colorimetric data by use of SLIC algorithm, then produces normal vector for each superpixel with a plane-fitting technique, and last performs classification of superpixels by a multi-class support vector machine technique. This scheme in reference [41] needs to place several colored markers in the scene prior to data acquisition, which is infeasible for some geological objects, such as alpine and gorge region which is often dangerous and inaccessible. Besides, although the SLIC algorithm with high efficiency is utilized in this reference, it still performs segmentation only by color vector and planar position of pixels, and ignores the local rock surface attitude, which is an important feature for geological objects. It may affect the rock surface extraction precision, especially for those regions with similar color and texture.

The MFF-SLIC algorithm developed in this paper performs superpixel segmentation not only by color and planar position features but also by spatial and some geological features, which guarantees the accuracy of complicated rock surface extraction.

3. Methodology

Comparing with other superpixel segmentation algorithms, the advantages of 2D SLIC algorithm are as follows: (1) with equally sized compact superpixels; (2) with low computational cost and few input parameters; (3) with an optimal tradeoff between under-segmentation error and boundary recall. However, from the definition of distance measure, it’s apparent that only the color and planar position of pixels are involved in, and the spatial geometric and some other attribute features which would have more of an effect on the segmentation results, are ignored. For geological objects, attitude is an important attribute to reflect the difference of the rock surface. Besides color and planar position, rock surface attitude is an essential component for distance measure to get a reliable segmentation result. Moreover, in some special circumstances, the spatial location may have more impact on geological objects segmentation. For example, pixels with similar color and proximal planar position but very far apart in the direction perpendicular to the image plane may be assigned into the same superpixel with 2D SLIC algorithm. Therefore, a new hybrid superpixel segmentation algorithm, MFF-SLIC, has been developed in this paper. The algorithm firstly describes a 1-N mapping relationship model for the TLS point cloud and digital images based on registration result. Instead of traditional distance measuring in 2D SLIC algorithm, a new multi-feature fusion distance measure model is defined, which combines color and position of pixels in digital images, local attitude of rock surface and spatial location of the TLS point cloud; Then an initial segmentation result consisting of a set of superpixels is obtained; the superpixels need to be further processed by hierarchical and optimized clustering in order to get the final result. The flowchart of MFF-SLIC algorithm is shown as Figure 2.

In the process of initial segmentation, clustering centers are optimized to the position with the minimum gradient in a local neighborhood so that they are not edge pixels or noise pixels. The gradient calculation formula for $3\times 3$ neighboring window is shown as Equation (4):

$$G\left(x,y\right)={\Vert I\left(x+1,y\right)-I\left(x-1,y\right)\Vert }^2+{\Vert I\left(x,y+1\right)-I\left(x,y-1\right)\Vert }^2,$$

(4)

where $I\left(x,y\right)$ means CIELAB color vector $\left(L,a,b\right)$ of pixel $\left(x,y\right)$; $\Vert ·\Vert$ denotes the norm of the color vector.

As the connectivity of pixels is not explicitly mandatory, some isolated pixels may exist in the end of the iteration for superpixel segmentation. A further connectivity processing is necessary to merge the isolated pixels into the nearest superpixels.

3.1. 1-N Mapping Relationship between TLS Point Cloud and Digital Images

As the point in the TLS point cloud may be projected to more than one image, this paper defines a 1-N mapping relationship between the TLS point cloud and digital images. The coordinate transformation mathematical model is the key to get the 1-N mapping relationship.

Suppose one point $P_i={\left(X_{\mathrm{TLS}},Y_{\mathrm{TLS}},Z_{\mathrm{TLS}}\right)}^T$ in the TLS point cloud, and the corresponding point $P_i^{\prime }={\left(X_{\mathrm{TLS}}^{\prime },Y_{\mathrm{TLS}}^{\prime },Z_{\mathrm{TLS}}^{\prime }\right)}^T$ in the SfM point cloud generated by digital images could be obtained according to the rigid transformation model of registration, which is shown as Equation (5).

$$P_i^{\prime }=\left[\begin{array}{c}{X}_{\mathrm{TLS}}^{\prime }\\ Y_{\mathrm{TLS}}^{\prime }\\ Z_{\mathrm{TLS}}^{\prime }\end{array}\right]=\mu \mathit{R}\left[\begin{array}{c}{X}_{\mathrm{TLS}}\\ Y_{\mathrm{TLS}}\\ Z_{\mathrm{TLS}}\end{array}\right]+\mathit{T}=\mu \mathit{R}{P}_i+\mathit{T},$$

(5)

where R and T respectively denotes the rotation and translation matrix; $\mu$ means the scale factor.

Let the rotation matrix ${\mathit{R}}_j$ and translation matrix ${\mathit{T}}_j$ of image $I_j$, then in the camera coordinate system of image $I_j$, coordinates of the point ${P_i^{\left(j\right)}}^{\prime }$ which corresponds to the point $P_i$ could be calculated by Equation (6):

$$P_i^{\left(j\right)\prime }=\left[\begin{array}{c}{X_{\mathrm{TLS}}^{\left(j\right)}}^{\prime }\\ {Y_{\mathrm{TLS}}^{\left(j\right)}}^{\prime }\\ {Z_{\mathrm{TLS}}^{\left(j\right)}}^{\prime }\end{array}\right]={\mathit{R}}_j\left[\begin{array}{c}{X}_{\mathrm{TLS}}^{\prime }\\ Y_{\mathrm{TLS}}^{\prime }\\ Z_{\mathrm{TLS}}^{\prime }\end{array}\right]+{\mathit{T}}_j={\mathit{R}}_j{P}_i^{\prime }+{\mathit{T}}_j,$$

(6)

All elements of ${P_i^{\left(j\right)}}^{\prime }$ are divided by $\left(-{Z_{\mathrm{TLS}}^{\left(j\right)}}^{\prime }\right)$, and denotes as Equation (7):

$$p=-\left[\begin{array}{c}{X_{\mathrm{TLS}}^{\left(j\right)}}^{\prime }/{Z_{\mathrm{TLS}}^{\left(j\right)}}^{\prime }\\ {Y_{\mathrm{TLS}}^{\left(j\right)}}^{\prime }/{Z_{\mathrm{TLS}}^{\left(j\right)}}^{\prime }\end{array}\right]=\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right],$$

(7)

Let f denote the camera focal length, and $k_1$ as well as $k_2$ denote radial distortion parameters, then coordinates of the point $p^{\prime }$ in image coordinate system corresponding to the point $P_i$ could be obtained by Equation (8).

$$p^{\prime }=\left[\begin{array}{c}x\\ y\end{array}\right]=f·r\left(p\right)·\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right]=f·r\left(p\right)·p,$$

(8)

where,

$$r\left(p\right)=1+k_1·{\Vert p\Vert }^2+k_2·{\Vert p\Vert }^4,$$

(9)

The image coordinate system o-xy takes the image principle point as the origin, the x axis is to the right, and the y axis is vertical upward. However, instead of the above image coordinate system, the pixel coordinate system o′-uv is commonly used in the process of calculation, the origin of which is located in the upper-left corner of the image. The relation between image coordinate system and pixel coordinate system is shown as Figure 3.

In Figure 3, it can be concluded that for the point $P_i$ in the TLS point cloud, its corresponding point $p^{″}$ in pixel coordinate system is expressed as Equation (10).

$$p^{″}=\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{c}\left(W/2\right)+x\\ \left(H/2\right)-y\end{array}\right],$$

(10)

where W and H respectively means the width and height of an image.

Both the intrinsic and extrinsic parameters for each image as well as transformation parameters between the TLS point cloud and SfM dense point cloud could be obtained with the G-Super4PCS registration algorithm introduced in reference [9]. According to all of the above process, the 1-N mapping relationship between the Point cloud and images could be concluded. The sketch map is shown as Figure 4, where P denotes one 3D point in the TLS point cloud, and $\left[p_1\dots p_N\right]$ respectively denotes the corresponding image point in digital images $\left[I_1\dots I_N\right]$.

3.2. Multi-Feature Fusion Distance Measure Model of MFF-SLIC

In the initial segmentation stage of MFF-SLIC algorithm, K clustering centers are initialized first, and the distance between each pixel in search area and clustering centers is respectively calculated by multi-feature fusion distance measure model based on a local k-means method, then each pixel would be assigned to the clustering center with minimum distance. Finally, the image I is segmented into K superpixels with special sematic characteristics by iteratively.

During the above process, instead of the traditional k-means method searching the whole image, the local k-means defines a local square search area with $2S\times 2S$. The difference of the two methods are shown as Figure 5.

The main function of the distance measure model is to distinguish dissimilar pixels and assign similar pixels into the same superpixel. The quality of the distance measure model determines the accuracy of superpixel segmentation.

With the 1-N mapping relationship of the TLS point cloud and digital images, not only the 2D image features such as color and planar position of pixels but also the 3D TLS point cloud features including spatial geometric as well as other characteristics should be considered for distance measure model. For the geological objects, the local attitude of rock surface is an important component for distance measure model construction. Besides, the 3D spatial location of pixels contributes to distinguishing pixels with similar texture and local attitude but are very far apart from each other in order that they would not be assigned into the same superpixel. Based on the above analysis, this paper defines a multi-feature fusion distance measure model according to 1-N mapping relationship between the 3D TLS point cloud and digital images, which is used to calculate the distance between unclassified pixels and clustering centers. The flowchart of distance measure model construction is shown as Figure 6.

3.2.1. Color Feature Calculation of Digital Image

The transformation from RGB color space to Lab color space needs to use intermediate variables XYZ. Therefore, the color feature calculation of digital images could be divided into two sections.

1. RGB to XYZ [42]

Suppose the three original color variables r, g and b, and the range of values is [0, 255]. A function is defined as Equation (11),

$$\gamma \left(x\right)=\{\begin{array}{ll}{\left(\frac{x+0.055}{1.055}\right)}^{2.4},& x>0.04045\\ \frac{x}{12.92},& \mathrm{other}\end{array},$$

(11)

where $\gamma \left(x\right)$ is used to improve contrast by non-linear hue adjustment. The three original variables r, g and b are transformed by Equation (11),

$$\left[\begin{array}{c}R\\ G\\ B\end{array}\right]=\left[\begin{array}{c}\gamma \left(r/255\right)\\ \gamma \left(g/255\right)\\ \gamma \left(b/255\right)\end{array}\right],$$

(12)

Let ${\left[R,G,B\right]}^{\mathrm{T}}$ is transformed to ${\left[X,Y,Z\right]}^{\mathrm{T}}$ by matrix M, which is shown as Equation (13),

$$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\mathit{M}\ast \left[\begin{array}{c}R\\ G\\ B\end{array}\right],$$

(13)

where,

$$\mathit{M}=\left[\begin{array}{ccc}0.4124& 0.3576& 0.1805\\ 0.2126& 0.7152& 0.0722\\ 0.0193& 0.1192& 0.9505\end{array}\right],$$

(14)

2. XYZ to Lab [42,43]

Suppose a function $f\left(y\right)$ as Equation (15),

$$f\left(y\right)=\{\begin{array}{cc}{y}^{1/3},& y>{\left(\frac{6}{29}\right)}^3\hfill \\ \frac{1}{3} {\left(\frac{29}{6}\right)}^{2}y+\frac{4}{29},& \mathrm{other}\hfill \end{array},$$

(15)

where $f\left(y\right)$ is also a correction function, then the calculation formula of color feature in CIELAB color space is described as Equation (16),

$$\left[\begin{array}{c}L\\ a\\ b\end{array}\right]=\left[\begin{array}{c}116f\left(Y/Y_0\right)-16\\ 500\left[f\left(X/X_0\right)-f\left(Y/Y_0\right)\right]\\ 200\left[f\left(Y/Y_0\right)-f\left(Z/Z_0\right)\right]\end{array}\right],$$

(16)

where $X_0,Y_0$ and $Z_0$ is respectively used for the normalization of $X,Y$ and $Z$. Generally, the value of ${\left[X_0,Y_0,Z_0\right]}^{\mathrm{T}}$ is set to ${\left[95.047,100,108.883\right]}^{\mathrm{T}}$.

3.2.2. Local Attitude Feature Calculation of TLS Point Cloud

The attitude of rock surface is usually described as dip angle $\theta$ and dip α, the value range of which is respectively $\left[0^{\circ },{90}^{\circ }\right]$ and $\left[0^{\circ },{360}^{\circ }\right]$. The values of $\theta$ and α could be calculated with the normal vector according to space geometry relation of rock surface. The normal vector of any point in a point cloud is calculated by the local plane fitted by points in neighborhood. The normal vector of the point is approximately equal to the one of the fitted local plane. According to the above, the local attitude of the point cloud could be obtained by the normal vector of the corresponding point.

The geometric relation of normal vector and local attitude is shown as Figure 7. In Figure 7a, the normal vector of the point cloud is denoted as $\mathit{n}=\mathit{O}\mathit{A}={\left(n_x,n_y,n_z,\right)}^{\mathrm{T}}$, and θ and ${\mathit{O}\mathit{A}}^{\prime }$ respectively means the dip angle and the direction of the dip in the local rock surface. In Figure 7b, the projection of the normal vector in horizontal plane is denoted as ${\mathit{n}}^{\prime }={\mathit{O}\mathit{A}}^{\prime }$, and α and ${\mathit{O}\mathit{A}}^{\prime }$ respectively means the magnitude and direction of the dip. The local attitude could be calculated as the following two ways because of the different values of the normal vector.

1. All three components of the normal vector are not equal to 0.

Under this condition, the relation of the normal vector $\mathit{n}$, dip α and dip angle θ is described as Equation (17).

$$\{\begin{array}{l}{n}_x=n_z·\tan \theta ·\cos \alpha \\ n_y=n_z·\tan \theta ·\sin \alpha \end{array},$$

(17)

A calculation formula for dip α and dip angle θ derived from Equation (17) by using the elimination method is shown as Equation (18).

$$\{\begin{array}{l}\theta ={\tan }^{-1}\left(\sqrt{n_x^2+n_y^2}/\left|n_z\right|\right)\\ \alpha ={\tan }^{-1}\left(n_y/n_x\right)\end{array},$$

(18)

where the sign of $n_x,n_y$ and $n_z$ determines the quadrant of dip α in the horizontal plane. The four geographic directions of east, south, west and north are respectively represented by E, S, W and N, which is shown as Figure 6b, then the relation of three components of the normal vector and the dip is listed in

Table 1. The relation of the normal vector and the dip.

${\mathit{n}}_{\mathit{x}}$	${\mathit{n}}_{\mathit{y}}$	${\mathit{n}}_{\mathit{z}}$	Dip (Direction)	${\mathit{n}}_{\mathit{x}}$	${\mathit{n}}_{\mathit{y}}$	${\mathit{n}}_{\mathit{z}}$	Dip (Direction)
>0	>0	>0	NE	<0	>0	>0	NW
>0	>0	<0	SW	<0	>0	<0	SE
>0	<0	>0	SE	<0	<0	>0	SW
>0	<0	<0	NW	<0	<0	<0	NE

.

2. No less than one component of the normal vector is equal to 0.

For some certain special rock surfaces, such as those developed horizontally and vertically, it is not suitable for Table 1. The calculation of local attitude of rock surface could be achieved by

Table 2. The relation of the normal vector and the dip for horizontally and vertically developed rock surface.

$\mathbf{The}\mathbf{Normal}\mathbf{Vector}\mathit{n}=\left({\mathit{n}}_{\mathit{x}},{\mathit{n}}_{\mathit{y}},{\mathit{n}}_{\mathit{z}}\right)$		Dip (Magnitude)	Dip (Direction)	Dip Angle
$\left(n_x=0,n_y\ne 0,n_z\ne 0\right)$	$\left(n_y\times n_z>0\right)$	--	N	Equation (18)
	$\left(n_y\times n_z<0\right)$	--	S
$\left(n_x\ne 0,n_y=0,n_z\ne 0\right)$	$\left(n_x\times n_z>0\right)$	--	E
	$\left(n_x\times n_z<0\right)$	--	W
$\left(n_x\ne 0,n_y\ne 0,n_z=0\right)$	--	Equation (18)	--	90°
$\left(n_x=n_y=0,n_z=1\right)$		--	--	0°
$\left(n_x=n_z=0,n_y=1\right)$		--	--	90°
$\left(n_x=1,n_y=n_z=0\right)$		--	--	90°

.

3.2.3. Distance Measure Model Construction of the MFF-SLIC

Suppose the fusion eigenvector for pixels $\left(I_i, I_j\right)$ in image $I$ with $M\times N$ is denoted as ${\left(L,a,b,x,y,\theta ,\alpha ,X,Y,Z\right)}^{\mathrm{T}}$, then for digital images, the calculation formula for distance of color components ${\left(L_i,a_i,b_i\right)}^{\mathrm{T}}$ and ${\left(L_j,a_j,b_j\right)}^{\mathrm{T}}$ is represented as Equation (19).

$$d_{Lab}=\sqrt{{\left(L_j-L_i\right)}^2+{\left(a_j-a_i\right)}^2+{\left(b_j-b_i\right)}^2},$$

(19)

The calculation formula for distance of pixel planar position components is described as Equation (20).

$$d_{xy}=\sqrt{{\left(x_j-x_i\right)}^2+{\left(y_j-y_i\right)}^2},$$

(20)

For the TLS point cloud, the local attitude (the dip angle ${\theta }_i$ and ${\theta }_j$, and the dip ${\alpha }_i$ and ${\alpha }_j$) of each point could be calculated as Equation (21).

$$\{\begin{array}{c}{d}_{\theta }=\left|{\theta }_j-{\theta }_i\right|\\ d_{\alpha }=\left|{\alpha }_j-{\alpha }_i\right|\end{array},$$

(21)

The distance for spatial location components of points in the TLS point cloud could be calculated as Equation (22).

$$d_{XYZ}=\sqrt{{\left(X_j-X_i\right)}^2+{\left(Y_j-Y_i\right)}^2+{\left(Z_j-Z_i\right)}^2},$$

(22)

All above distances defined in various standards have different effects on the construction of multi-feature fusion distance measure model, which should be normalized according to some certain criterions. In this paper, the normalized distance measure model constructed by the eigenvector integrated with multiple features is described as Equation (23).

$$D_S=d_{Lab}/m+d_{xy}/S+d_{\theta }/\Delta \theta +d_{\alpha }/\Delta \alpha +d_{XYZ}/\Delta XYZ,$$

(23)

where S means the grid interval of superpixels; m means the compactness control variable of superpixels; $\Delta \theta$ and $\Delta \alpha$ respectively means the maximum distance of local dip angle and dip.

The TLS point cloud and digital images acquired independently are not always associated in a one-to-one relationship. The TLS point cloud is composed of a set of points distributing as a certain density, while the digital image consists of a series of pixels distributing densely and uniformly. For those pixels without corresponding points in the TLS point cloud, the spatial features could be obtained by use of the inverse distance weighting method, the expression of which is shown as Equation (24).

$$\tilde{Z}\left(x_0,y_0\right)={{\displaystyle \sum }}_{i=1}^n{\omega }_{i}Z\left(x_i,y_i\right),\mathrm{where}{\omega }_i=\frac{1}{d_i}/\left({{\displaystyle \sum }}_{i=1}^n\frac{1}{d_i}\right),$$

(24)

3.3. Hierarchical and Optimized Clustering for Superpixels

The result of superpixel segmentation is a set of superpixels composed by pixels with similar features such as the local attitude. Although these superpixels could express the boundary of rock surface well, it is difficult to accurately describe the whole structure of rock surface. In order to get the optimal results, it is necessary for initial superpixels to make a further process. This paper provides a hierarchical and optimized clustering strategy for superpixels, which takes superpixels as nodes, constructs the Region Adjacency Graph (RAG) and Nearest Neighbor Graph (NNG) by defining a regional dissimilarity function for superpixels based on the attitude of rock surface, and achieves superpixels merging. The flowchart of hierarchical and optimized clustering for superpixels is shown as Figure 8. In this paper, the attitude of a superpixel corresponds to the histogram peak of the dip distribution of all pixels contained in the current superpixel.

3.3.1. Definition of Dissimilarity between Superpixels

The superpixel segmentation is a process to assign a superpixel label for each pixel, while the clustering of superpixels refers to reassign a new label for each superpixel, i.e., superpixels with similar characteristics would be divided into the same group. In this paper, regional dissimilarity is defined as the criterion of superpixels clustering.

Suppose the image I is segmented into n superpixels, which is represented as the following set S:

$$S=\left\{S_1,S_2,\cdots ,S_n\right\},$$

(25)

For the ith superpixel $S_i\left(i=1,\cdots ,n\right)$, the corresponding set of pixels is denoted as Equation (26).

$$S_i=\left\{p_1^{\left(i\right)},p_2^{\left(i\right)},\cdots ,p_m^{\left(i\right)}\right\},$$

(26)

where m means the number of pixels in the ith superpixel; $p_k^{\left(i\right)}$ denotes the kth pixel in the ith superpixel, $1\le k\le m$.

The regional dissimilarity relation between superpixels $S_i$ and $S_j$ could be described as Equation (27).

$$\left(S_i,S_j\right)=1-\left(S_i\cap S_j/S_i\cap S_j\right),$$

(27)

where the range of value is $\left[0,1\right]$, and the more similar the two superpixels, the closer the value and 0, else, the closer the value and 1.

In order to describe the regional dissimilarity in a global sense, Equation (27) should be combined with Equation (28).

$$\mathrm{sgn}\left(S\right)=\{\begin{array}{cc}1,\hfill & S_i \mathrm{is} \mathrm{next} \mathrm{to} S_j\\ +\infty ,\hfill & else\end{array},$$

(28)

Suppose the dip angle and dip sets of pixels in superpixel $S_i$ are denoted as Equation (29), and the corresponding histograms could be obtained.

$$\begin{array}{c}\left\{{\theta }_1^{\left(i\right)},{\theta }_2^{\left(i\right)},\cdots ,{\theta }_m^{\left(i\right)}\right\}\\ \left\{{\alpha }_1^{\left(i\right)},{\alpha }_2^{\left(i\right)},\cdots ,{\alpha }_m^{\left(i\right)}\right\}\end{array},$$

(29)

Denote the attitude of superpixel $S_i$ as $\left({\theta }_i,{\alpha }_i\right)$, then values of which are respectively equal to peaks of above histograms. The regional dissimilarity function between superpixels $S_i$ and $S_j$ is defined as Equation (30).

$$\mathrm{diff}\left(S_i,S_j\right)=\mathrm{diff}\left({\theta }_i,{\theta }_j\right)+\mathrm{diff}\left({\alpha }_i,{\alpha }_j\right),$$

(30)

3.3.2. Topological Description for Superpixels with RAG

Theoretically, the clustering of superpixels based on the regional dissimilarity needs an exhaustive comparison between superpixels, the complexity of which is $O\left(n^2\right)$, here n means the number of superpixels. However, only the adjacent superpixels are possible to be merged. Therefore, this paper introduces RAG to describe the topological relation between superpixels, on the basis of which only adjacent superpixels would be compared, and the complexity is effectively reduced to linear. Figure 9a,b respectively shows an example of superpixel segmentation result for an image as well as the corresponding RAG.

The RAG of the superpixel set $\left\{S_1,S_2,\cdots ,S_n\right\}$ could be described by nodes and edges. Each node represents a superpixel, and each edge connects two adjacent superpixels. An example of RAG composed of nine superpixels is shown as Figure 10. In Figure 10b, $\left\{V_1,V_2,\cdots ,V_9\right\}$ means the set of nodes, i.e., the set of superpixels, and $\left\{E_1,E_2,\cdots ,E_9\right\}$ means the set of edges. In Figure 10a, the superpixel $S_1$ corresponds to the node $V_1$ in Figure 10b; adjacent superpixels of $S_1$ is $S_2$, $S_5$ and $S_6$, and $V_2$, $V_5$ and $V_6$ are corresponding nodes, and edge $E_1$, $E_2$ and $E_3$ respectively connect $S_1$ to $S_2$, $S_5$ and $S_6$.

The weight of each edge is calculated according to the regional dissimilarity function. Suppose attitudes of two adjacent superpixels $S_i$ and $S_j$ are denoted as $\left({\theta }_i,{\alpha }_i\right)$ and $\left({\theta }_j,{\alpha }_j\right)$, the number of pixels which contain are respectively denoted as $m_i$ and $m_j$. Let

$${\eta }_{ij}=\left(m_i·m_j\right)/\left(m_i+m_j\right),$$

(31)

The calculation formula of weight for edge $E_k$ connecting two adjacent superpixels is written as Equation (32).

$$W_{E_k}=\frac{1}{{\eta }_{ij}·\left(\left({\theta }_i-{\theta }_j\right)/{\Delta }_{{\theta }_{\max }}+\left({\alpha }_i-{\alpha }_j\right)/{\Delta }_{{\alpha }_{\max }}\right)},$$

(32)

where ${\Delta }_{{\theta }_{\max }}$ and ${\Delta }_{{\alpha }_{\max }}$ respectively means the maximum of the dip angle and the dip between all adjacent superpixels.

For the RAG in Figure 10b, if the length of the edge is used to describe its weight, it can be concluded that the weight of the bold black edge $E_{11}$ is minimum, which indicates that the node $V_5$ and $V_6$, i.e., the superpixel $S_5$ and $S_6$ should be merged. The updated RAG is shown as Figure 11.

3.3.3. Fast Clustering for Superpixels Based on Nearest Neighbor Graph

The clustering of superpixels based on RAG takes a bottom-up way and a hierarchical merge strategy. During the above process, the merging as the minimum of edge weight involves in repeated update and storage of RAG, the computation of which is too expensive. Therefore, this paper combines NNG with RAG in order to get the required results with higher efficiency and less computation.

From Figure 10, it is obvious that one superpixel node may correspond to several adjacent nodes and edges, and these edges have no direction so that the RAG is a scalar graph. The NNG keeps nodes the same as the RAG, while scalar edges are converted to vector ones in the NNG. Figure 12 shows the conversion process between RAG and NNG.

In the NNG, all adjacent superpixels are connected by scalar edges, and the length of each edge reflects the regional dissimilarity of the adjacent superpixels. Compare weights of all scalar edges, preserve the one with minimum weight, and specify the direction from the current node to the adjacent one with minimum weight. As one node corresponds to no more than one edge, the total number of edges is no more than the number of nodes. For example, in Figure 12, node $V_1$, $V_2$ and $V_7$ respectively correspond to edge $E_2$, $E_6$ and $E_{14}$; node $V_3\leftrightarrow V_4$, $V_5\leftrightarrow V_6$ and $V_8\leftrightarrow V_9$ respectively share edge $E_7$, $E_{11}$ and $E_{15}$. In comparison of RAG and NNG, it can be concluded that the number of edges in NNG is less than the one in RAG although they have the same number of nodes, and edges of NNG are vector. An important characteristic of NNG is that for any NNG, there exists at least one bidirectional edge which connects the two adjacent superpixels with the minimal regional dissimilarity. Therefore, the clustering of superpixels could be achieved by searching all bidirectional edges in NNG.

Based on all the above analysis, superpixels could be clustered hierarchically and efficiently by combination of RAG and NNG.

4. Experiment and Discussion

In order to verify the segmentation algorithm based on the registration result obtained by the G-Super4PCS algorithm proposed by Zhang et.al [9], this paper also takes the columnar basalt data from Guabushan Geopark in Nanjing, China as the experimental data. The digital images are acquired by Fujifilm X-T10 digital camera, and the TLS point cloud is acquired by FARO Focus^3D X330 terrestrial laser scanner. The test area is respectively shown as Figure 13a (the digital image) and Figure 13b (the point cloud).

Firstly, the 1-N mapping relationship between the point cloud and images is obtained based on the registration result. According to mapping relationships between the TLS point cloud and different images, colors and textures of images could be added to the TLS point cloud. The correspondences of the TLS point cloud and images (a)~(g) in 3D visualization is shown as Figure 14.

Secondly, distance measure model calculation is an important part of MFF-SLIC algorithm. For color feature calculation of digital image, RGB color space should be transformed into Lab color space first according to Section 3.2.1. There are 251,861 3D points in the experimental TLS point cloud. Calculate the local attitude of each 3D point in the point cloud according to Section 3.3.2. Figure 15a shows an intuitional distribution of local attitudes for 251,861 points (red points in the figure), which is a pole plot projected by equal angle. The range values of the dip angle and dip are respectively $\left[0^{\circ },{90}^{\circ }\right]$ and $\left[0^{\circ },{360}^{\circ }\right]$. The radius of the pole plot is ${90}^{\circ }$, so the distance between the blue center of the circle and each red point denotes the dip angle of the red point, and correspondingly, its dip is equal to the angle between due North and the direction in which the red point is connected to the blue center of the circle. For example, dips angle and dips of two yellow cross marks in the pole plot are respectively $\left[{65}^{\circ },{23}^{\circ }\right]$ and $\left[{40}^{\circ },{273}^{\circ }\right]$

Figure 15b shows the contour plot generated from local attitudes for 251,861 points. The contour plot describes poles density distribution characteristics and laws by use of isopleth based on the pole plot. The color gradations in Figure 15b represent fissure developmental state.

In this paper, the overlapping area of Figure 14a is taken as an example for the following segmentation and clustering experiments of superpixels.

In digital images, attitudes of those pixels corresponding to 3D points in the TLS point cloud could be obtained according to the 1-N mapping relationship between the two data sources. For those pixels without corresponding 3D points in the TLS point cloud, attitudes could be calculated via interpolation with Equation (24) in Section 3.2.3. The 2D visualization of attitudes for a digital image is shown as Figure 16.

The next is to construct the multi-feature fusion distance measure model with all feature components of both the two data sources, and perform MFF-SLIC hybrid super pixel segmentation. The number of superpixels is a key parameter which affects the result of segmentation. This paper makes a comparison between 2D SLIC algorithm and MFF-SLIC algorithm by setting various superpixel number. Figure 17 shows the comparison result when the number of superpixels is set to 500. In this figure, the red and green lines respectively represent superpixel segmentation boundaries of the two algorithms. It is observed from the two local rectangle regions A and B that the minute fissure of rock stratum is processed more accurately with the MFF-SLIC algorithm, while under-segmentation could not be avoided well for the 2D SLIC algorithm, such as the fissures marked with cyan, the color and texture of which are similar to the surrounding area. Therefore, it is difficult to distinguish without the local geological attitude information. It is obvious that the MFF-SLIC algorithm considering local geological features has a better performance on quality segmentation. Moreover, superpixels segmented by the use of MFF-SLIC algorithm are all with geological attributes, which would make the next clustering process achieve higher efficiency and accuracy.

Figure 18 shows the comparison of segmentation results with 1000 superpixels. It can be seen that MFF-SLIC still performs better boundary consistency.

According to the above experimental results, it can be indicated that compared with 2D SLIC, the MFF-SLIC algorithm has precise and reliable advantages. It is worth mentioning that the difference between the two algorithms will not be obvious when the number of superpixels is too large. Figure 19 shows the comparison of segmentation results for the two algorithms when the number of superpixels is set to 2500. The size of the superpixel is so small that it is difficult to evaluate the performance of the two algorithms. Therefore, the number of superpixels may affect both efficiency and precision of the algorithm, which should be reasonably set according to image size and superpixel size. Too few superpixels would cause under-segmentation, which may decrease the rock surface extraction precision. Conversely, too many superpixels would cause over-segmentation, which may reduce segmentation efficiency. Besides, the size of superpixels is so small that the attitude difference between superpixels is not obvious, so it is difficult to distinguish valid rock surface.

Histograms for attitude distributions of some part superpixels are shown as Figure 20, the left column shows dip angle histogram distributions, and the right column shows dip histogram distributions. The horizontal axis represents the value of the dip angle or the dip, and the vertical axis represents pixel numbers corresponding to different dip angles or dips in a superpixel. All histograms indicate that attitudes of pixels in a superpixel are approximately consistent, i.e., most pixels in a superpixel have similar attitude characteristics, which demonstrates that the result of superpixel segmentation is reliable and accurate.

The test image size is 784 by 845 pixels, which is approximately determined by overlapping area, and the size and the number of superpixels are experimentally set to 442 (pixels) and 1500 (superpixels). The corresponding RAG and NNG for superpixel segmentation result are shown as Figure 21. The blue nodes represent superpixels, each red line segment connects two adjacent superpixels, and the weight of the red line segment is reflected by dissimilarity of the adjacent superpixels.

With the RAG and NNG, the hierarchical and optimized clustering of superpixels could be achieved according to Section 3.3. Figure 22 shows the superpixels clustering result for the test image. After 20 iterations, the number of superpixels reduced to 118.

The final segmentation results could be obtained through a further congruency analysis to all overlapping areas between the two data sources. Therefore, 61 structural planes are successfully extracted according to the segmentation results, which are marked in different colors, and the 3D visualization is shown as Figure 23.

In order to verify the reliability of the MFF-SLIC hybrid superpixel segmentation algorithm, this paper makes a statistic and comparison to rock surface extraction results respectively obtained by the proposed segmentation method and the manual measuring method by geologists. Comparing the results shows that the max dip difference and max dip angle difference are both less than 6°, and the accuracy meets the requirement of geological analysis. Moreover, the contour plots for the two rock surface extraction results by the proposed method and the manual measuring method by geologists are shown as Figure 24, and it is not difficult to see that they are basically identical, i.e., the distributions of structural planes are consistent.

All the above experimental results demonstrate that the MFF-SLIC hybrid superpixel segmentation algorithm could achieve reliable and accurate segmentation by integrating multiple features from the TLS point cloud and digital images without any manual interactions, which provides important data materials for further geological statistics and analysis. Besides, most existing segmentation methods involves many parameters so that the quality of segmentation results heavily depends on the choice of parameters. Relatively, the whole process of the MFF-SLIC algorithm needs few parameters, which improves the automatization of segmentation.

5. Conclusions

Compared to traditional segmentation methods in pixels, superpixel segmentation methods could well consider local regional correlation, and clustering on this basis could improve segmention efficiency and accuracy. However, most existing superpixel segmentation algorithms are presented for 2D images, which only depend on the color and pixel planar position of image to define the distance measure. For geological objects, they may also be more affected by other characteristics such as the spatial location and rock surface attitude. Such segmentation methods are not always appropriate for geological objects, and especially for those special circumstances with similar image features but various spatial or other features, it is difficult to guarantee reliable and accurate results. In this paper, a new hybrid superpixel segmentation algorithm, MFF-SLIC, for rock surface extraction by integrating the TLS point cloud and digital images is proposed. In order to take full use of characteristics of Multi-Source RS data, this algorithm constructs the 1-N mapping relationship between TLS point cloud and digital images based on the registration results, defines a new distance measure model integrating multiple features including the color and planar position of pixels in digital images as well as the spatial location and local rock surface attitude of 3D points in the TLS point cloud, and introduces a hierarchical and optimized clustering strategy of superpixels which defines a regional dissimilarity as adjacent superpixels clustering criterion and combines RAG and NNG to improve clustering efficiency. The algorithm was verified by use of the columnar basalt data from Guabushan Geopark in Nanjing, China. The results demonstrate that the proposed method could achieve rock surface extraction with high precision and efficiency, the result of which would be prepared for further geological statistics and analysis.