Rock Fracture Centerline Extraction based on Hessian Matrix and Steger algorithm

Abstract

The rock fracture detection by image analysis is significant for fracture measurement and assessment engineering. The paper proposes a novel image segmentation algorithm for the centerline tracing of a rock fracture based on Hessian Matrix at Multi-scales and Steger algorithm. A traditional fracture detection method, which does edge detection first, then makes image binarization, and finally performs noise removal and fracture gap linking, is difficult for images of rough rock surfaces. To overcome the problem, the new algorithm extracts the centerlines directly from a gray level image. It includes three steps: (1) Hessian Matrix and Frangi filter are adopted to enhance the curvilinear structures, then after image binarization, the spurious-fractures and noise are removed by synthesizing the area, circularity and rectangularity; (2) On the binary image, Steger algorithm is used to detect fracture centerline points, then the centerline points or segments are linked according to the gap distance and the angle differences; and (3) Based on the above centerline detection roughly, the centerline points are searched in the original image in a local window along the direction perpendicular to the normal of the centerline, then these points are linked. A number of rock fracture images have been tested, and the testing results show that compared to other traditional algorithms, the proposed algorithm can extract rock fracture centerlines accurately.

1. Introduction

In recent years, there are more and more geological damages such as mountain landslide caused by rock fractures, which results in huge economic losses and casualties [1-2]. If the relevant fractures can be detected and evaluated earlier, the disasters might be avoided or the losses would be declined. One of the fundamental steps in rock fracture image analysis is the fracture detection, which directly affects the efficiency and the accuracy of the subsequent measurements and analyses. Compared with manual measurement, the rock fracture detection by image processing is safer, more efficient and more precise [3]. At present, there is still no standard image processing algorithm for the rock fracture detection, therefore, the research for a new algorithm is significant both in theoretical and in practical. In general, rock fractures can be found out in slopes, underground tunnels, underground caverns, mines, pavements, mountains and so on, therefore, there are many literatures describing this research topic.

The rock fracture connectivity was analyzed by the fracture intersections and the advisable classification grade of rock surfaces was made based on intersections per fracture [4]. The fracture orientation in 3-D was estimated by tracing the fractures of several rock surfaces, and then the 3-D characteristics could be determined by the relevant 2-D features [5]. The gray-scale transformation and Canny edge detection were extended from 2-D to 3-D for processing 3-D seismic rock fractures, which achieved the satisfactory results [6]. The classification algorithms based on Sobel, Marr-Hildreth, Harr wavelet and texture analyses were respectively applied to classify the damage and non-damage of rock fracture structures, which verified that the algorithm, based on texture analysis, could reach to the best results [7]. The fracture detection on coal surfaces and the estimation of fluid inclusion traces were conducted but the precisions of both detection and regression were to be further improved [8]. The fractures on 3-D were enhanced by the multi-scale filter based on Hessian matrix, which gained the satisfactory results but was impacted by selected scales [9]. The rock fracture centerlines could be extracted by image processing and support vector machine [10], and the algorithm was able to extract the fracture centerline in a test, but it is not suitable for all cases.

Some algorithms were studied for rock fracture measurement rather than the detection which might be by manual [11]. A single fracture was well traced and its width was easily measured but the algorithm was only suitable for the fractures without intersection [12]. The fracture edges were detected by the algorithm based on rough set, the detection results showed that the fracture edges are not continuous and the resulted image includes a lot of noise [13]. The generalized Hilbert transform was adopted to delineate fractures, which is strongly anti-noise; however, the results contained many holes which made difficulty for rock fracture measurement and analyses [14]. Fractures were delineated by fractal geometry but the results were dissatisfactory since rock surface is too rough [15]. The wavelet transformwas utilized to detect fractures, which could locate faint fractures and involved spurious fractures [16]. A minimum spanning tree algorithm was put forward to extract fractures, the results were continuous fractures of one pixel width, but it was not suitable for fractures with intersections too [17]. The faint discontinuous fractures were linked according to the connected domain correlation; nevertheless, the reasonable parameters might be selected during the local histogram equalization and noise reduction [18]. Fractures were extracted by valley edge detection of fractional integral, which obtained the satisfied results; however, the parameter of fractional integral should be adjusted and the post processing steps might be added [19]. The pulse coupled neural network based on the minimum error principle was used to extract fractures, but the complicated algorithm was not suitable for the complexly distributed fractures [20]. It was assumed that the light intensity profile of the fracture cross section perpendicular to the centerline of a true fracture resembles a Gaussian distribution function [21]. Firstly, it used a filter to enhance the fractures with characteristics of symmetry and line-like; then, after the gray level image binarization, the fracture centerlines were extracted; finally, 1-D Savitzky-Golay smoothing/differentiation filter was adopted to detect fracture edges along the direction perpendicular to the centerline. The algorithm was valuable for simple fractures while not applicable to the fractures with intersections; meanwhile, the assumption above was not always tenable. The rock fractures could also be detected by fractional differential [22], local gray value minima [23] and color clustering and quaternion based edge detection algorithms [24]. Since rock fracture or pavement crack images are often unclear or fuzzy, some image enhancement algorithms can be used as preprocessing [25-26].

The paper proposes a novel algorithm. It firstly obtains a mask (template) for fracture enhancement, then makes centerline extraction roughly, and finally does fracture delineation accurately. To gain a fracture mask, Hessian Matrix and Frangi filter are firstly adopted to enhance the curvilinear structures, then after image binarization, the spurious-fractures and noises are removed by synthesizing the area, circularity and rectangularity. For the second step, for extracting a fracture centerline, according to the mask above, Steger algorithm is used on the original image to detect and link the centerline points, and then the centerline segments are linked based on the distances and angle differences. In the last step, in order to detect fracture edges accurately, the rough result is improved in the original image, which including the noise removal, centerline points searching and centerline gap linking.

 

2. Image Enhancement on Hessian Matrix at Multi-scales

For a 2-D matrix I : Ω → R,Ω ∈ R2, the definition of Hessian matrix corresponding to I at point p(x,y) is:

Hessian matrix is a real symmetric matrix. If a 2-D image is represented by I, the eigenvalues λ1 and λ2 of Hessian matrix can be used for constructing a filter to enhance the curvilinear structures. The parameters λ1 and λ2 are calculated as:

Where, K and Q are calculated by the formula

The eigenvectors e1 and e2 corresponding to eigenvalues λ1 and λ2 are important shape feature parameters. For a 2-D image, the maximum and minimum of a surface curvature are given by λ1 and λ2, meanwhile, the corresponding directions are determined by e1 and e2. As shown in Fig.1, a dark synthetic line is on the bright background; furthermore, the intensity profile of the synthetic cross section perpendicular to the centerline resembles a Gaussian distribution function. Assume |λ1|≤|λ2|, Fig.1 shows e1 and e2 at the center point. Corresponding to λ1 which is close to 0, the surface is not bending along the tangent plane; while corresponding to λ2 which has a large absolute value, the surface is bending along the normal plane.

Fig. 1.2-D synthetic line where the intensity profile of the cross section perpendicular to the centerline resembles a Gaussian distribution function

For a given signal I(x,y), its linear scale-space representation is a family of signals I(x,y;σ). The definition of I(x,y;σ) is:

Where, g(x,y;σ) is a Gaussian kernel, ∗ represents the convolution operation, σ is the scale parameter which is the deviation of Gaussian kernel. I(x,y;0) = I(x,y). indicates that when σ = 0, linear scale-space representation of I(x,y) is I(x,y) itself. As σ larger, I(x,y;σ) is the result by smoothing I(x,y) with a larger Gaussian kernel.

Fig. 2 lists the relationship between eigenvalues of Hessian matrix and shape structures. For a 3-D tubular structure, if a voxel belongs to the tubular, the eigenvalues of Hessian matrix at the point satisfies |λ1|≈0, |λ1|<<|λ2|,|λ1|<<|λ3|,λ2 ≈ λ3. For a 2-D curve structure, if a pixel belongs to the curve, the eigenvalues of Hessian matrix at the point meets |λ1|≈0, |λ1|<<|λ2|.

Fig. 2.Relationship between eigenvalues of Hessian matrix and shape structures (H = high, L = low, +/- represent the sign of the eigenvalue, the eigenvalues are ordered as: |λ1|≤|λ2|≤|λ3|)

According to Fig. 2, Frangi [27] proposed a method to extract 3D tubular structures. The algorithm is based on the following probability function:

Where, α,β,c are the linear filter parameters to control the sensitivity of RA,RB,S, 0 <νo(σ) < 1, RA = |λ2|╱|λ3|, . The probability function for 2-D curve structure is as follows:

In light of scale-space theory, the output of Frangi filter reaches the maximum when the curve width w matches the scale σ. The multi-scale enhancement algorithm based on Hessian matrix adopts various σ to compute νo and chooses the maximum as the final output f(x,y) = max[νo(x,y;σ)]. Frangi filter and its improved version were widely used for blood vessel enhancement [28].

The curvilinear structures are obviously enhanced by Frangi filter; however; some spurious fractures are also removed in this operation. After image binarization, the fractures have the large areas and obvious line-like structures while the spurious fractures have small areas and obsolete line-like structures. Hence, for the fracture with no intersection, spurious fractures can be removed directly by synthesizing the areas, circularity and rectangularity; for the fracture with intersections, they should be extracted before the spurious fracture removal.

Area is a metric of object size and is an important geometry feature. It plays a significant role in image analysis and is often used as a threshold for removing false objects. The common method of area calculation is accumulating the pixel number in the interior domain:

Where, f(x,y) is a binary image, f(x,y) = 1 represents for objects, f(x,y) = 0 represents for the image background.

Circularity is a shape feature widely applied in image analysis. The definition commonly used is:

Where, P and S are perimeter and area of an object respectively. When an object is of a circulated shape, C = 4π. The more circular the object is, the smaller C becomes; the more slender the object is, the larger C turns.

Rectangularity is another important shape feature which is usually applied to distinguish the slender object from the spheroidal/quadrate object. One of the parameters named length to width ratio is defined as:

Where, WMER and LMER are the width and length of the minimum circumscribed rectangle of an object, respectively.

 

3. Fracture Centerline Extraction on Gray Level

3.1 Feature Point Detection of Fracture Centerline

Hessian matrix was applied to track the vessel centerlines in Coronary Arteriograms according to the curvilinear structures [29]. The centerlines in other objects can also be extracted based on Hessian matrix [30].

Steger [31] drew a conclusion by studying the scale-space properties of 1-D curve. The second order derivative reaches to the extreme as , where, w is the width of 1-D curve and σ is the deviation of the Gaussian kernel. The conclusion can be extended to 2-D ridge/valley which can be modeled as s(t). Suppose n(t) is the normal of s(t). If a point belongs to the feature point of ridge/valley, the first order derivative at this point along n(t) equals to 0, meanwhile, the second derivative at this point reaches to the extreme.

The partial differentials of Gaussian kernel are defined as:

Where, gσ(x) and gσ(y) are 1-D Gaussian kernels with deviation σ. Firstly, image I(x,y) is convoluted with these kernels to obtain the partial differentials of I(x,y) rx(x,y),ry(x,y),rxx(x,y),rxy(x,y),ryy(x,y):

Secondly, Hessian matrix is constructed as:

Thirdly, the eigenvalue λ with maximal absolute value and its corresponding eigenvector (nx,ny) are calculated. (nx,ny) satisfies ║(nx,ny)║2 = 1, and gives the normal direction. The zero-crossing position of the first derivative along (nx,ny) can be determined by the Taylor polynomial interpolation:

Where,

If , the pixel is the feature point of the centerline, furthermore, if the second derivative reaches to the extreme, the pixel is considered as a salient point of the centerline.

3.2 Feature Point Linking in Fracture Centerline

Three kinds of feature point information of the fracture centerline can be obtained as: the direction (nx,ny) = (cos α,sin α), the intensity is the second derivative at the point, and the sub-pixel position (px,py) is the zero-crossing position of the first derivative along (nx,ny). According to these kinds of information, the linking algorithm suggested by Steger [31] is as follows:

There are many curves in a fracture image, each of which is corresponding to a maximum of the second derivative extremes. These curves are descending ordered by their corresponding maximum. For every element/curve in the descending ordered array, the points linking process is: firstly, the algorithm locates the feature point at which the second derivative extreme reaches to the maximum, and it considers the point as the current feature point; secondly, the algorithm searches the candidates and chooses the next feature point according to n⊥ of the current feature point; finally, the process is executed until there is no feature point which satisfies the constraint. The feature point at which the second derivative extreme reaches to the maximum is usually not an endpoint of the curve; therefore, the process above can only complete the feature point linking on one side of the feature points at which the second derivative extreme reaches to the maximum, therefore, to complete the feature point linking on the other side, the same process is operated along the direction -n⊥ of the feature points at which the second derivative extreme reaches to the maximum. The curves are linked in order until a maximum of the second derivative extreme of a certain curve is lower than the high threshol high_threshold.

The scheme of searching and choosing a feature point is as follows. (nx,ny) = (cos α,sin α) is accurate enough to simply consider three pixels in the eight neighborhood of the current feature point as the candidates. For example, if the position is (qx,qy) and the direction is α ∈ (−π/8,π/8], (qx + 1,qy - 1), (qx + 1,qy) and (qx + 1,qy + 1) are candidates. The distances dis = ║p1 - p2║2 and angle differences dif = |α1 - α2|,dif ∈ [0,π/2] are calculated for every candidate and the point with minimal dis + dif will be considered as the linking point if its corresponding second derivative extreme is higher than the low_threshold.

3.3 Feature Segment Linking in Centerline

However, as depicted above, many centerline segments of a fracture can be obtained. We can refine the segments: for every segment, there are at most two points in its eight neighborhood. For the convenience of further measurement, the segments are linked based on the thresholds of the distance and angle differences. The process is:

Firstly, for every endpoint in image I, the position (x,y) and the direction angle are marked. The angle calculation process is: set the endpoint (x,y) as the start point; backtrack n points on the segment; calculate the mean of coordinates , where, → represents rotation. It is notified that if the other endpoint is encountered during the backtracking process, all points on the segment are used for mean computation. Secondly, for every endpoint pcurrent of (x,y), the other endpoints are searched in the range of [x - δ,x + δ]×[y - δ,y + δ]; if there exist other endpoints, the distances dis and angle differences dif are calculated from the current endpoint pcurrent to all the other endpoints. Finally, among all the other endpoints of which both dis and dif satisfy the setting thresholds of the distance and angle difference, the endpoint pselected with minimal dis + dif is selected and pcurrent is connected to pselected by a line.

3.4 Improvement of Fracture Mask on Centerline Extraction

The idea of taking advantage of the mask for centerline extraction originated from [32], which made use of the retinal vessel mask to extract centerline.

In the process of the feature point detection, for every pixel in an original image, firstly, the convolutions of the pixel with different partial differentials of Gaussian kernel are implemented; secondly, Hessian matrix is constructed for the point; finally, the point is judged whether it is the feature point or not. It is a time-consuming process. In this study, the feature point of the centerline is only detected and linked on the original image in a range of mask, so the algorithm speed is intensively accelerated while the algorithm accuracy is not decreased.

In the process of feature point detection, the mainly used parameter is σ, which is set as the maximal scale parameter σmax in the multi-scale enhancement algorithm. In the process of feature point linking, the involved parameters are high_threshold and low_threshold. The fractures on the original image in a range of mask are to be extracted as much as we can, so the high_threshold is set as the value which is slightly greater than 0. As we tested five different values 0.05, 0.08, 0.1, 0.12, 0.15 for more than 65 rock fracture images, we set high_threshold = 0.01 and low_threshold = 0 for all fracture images in this study. Therefore, compared with the algorithm operated on the whole original image, Steger algorithm implemented in the range of mask omits the tedious parameter setting procedure.

 

4. Experiments and Analyses

To verify the availability of the new algorithm, more than 100 rock fracture images and about 50 pavement crack images were tested, and the several typical images in this paper are presented as follows. For all the images, the algorithm studied is compared with some widely used algorithms, such as Canny detection [33], auto thresholding [34], Valley edge detection [21,25], Minimum spanning tree (MST) [27], Clustering analysis [35], and Fuzzy C-Means (FCM) based segmentation algorithms [36]. For all the images, as tested, in Frangi filter, the threshold parameters are set as high_threshold = 0.01 and low_threshold = 0 as discussed in the above section.

Fig. 3(a) presents a simple pavement crack image, where, the pavement surface is rough as usual, the illumination is even, and only one crack clearly shows up.

Fig. 3.Procedure of centerline extraction on pavement for a simple crack image: (a) Original crack image; (b) Multi-scale enhancement; and (c) Centerline point detection and linking.

First, the crack image is enhanced by Frangi filter, which obtains the satisfactory result by choosing only one scale, the crack is strengthen as shown in Fig. 3(b); and then after image binarization by Otsu thresholding [34], the noise are filtered away by simply using an area threshold (20 pixels), finally the centerlines are extracted by Steger algorithm, as shown in Fig. 3(c).

In the comparison between our new algorithm and conventional algorithms on this image, the conventional algorithms involve three different threholdings (gray level similarity), two edge detection (discontinuity) algorithms, and a graph based algorithm (MST).

For comparing to Fig. 3(c), as shown in Fig. 4(a), Even the crack is extracted roughly, Otsu thresholding results in a lot of noise which makes crack tracing hard in post processing; Iterative thresholding and Entropy thresholding [34] can produce less noise as shown in Fig. 4(b-c), but they create some gaps in the crack, which can also make crack delineation difficult; Canny edge detector gives crack boundary location well, but many spurious edges and noise appear in the resulted image (Fig. 4(d)); Valley edge detector can extract the centerline of the crack, the noise cannot be suppressed (Fig. 4(e)); MST detects the most parts of the crack, but it produces some spurious cracks; and the new algorithm can make a good detection result as shown in Fig. 3(c). By the comparison, we conclude that for such a crack image, the new algorithm can achieve more excellent result than others. In the aspect of anti-noise, the new algorithm overcomes the background noise affections, which cannot be completely filtered out by the conventional algorithms.

Fig. 4.Comparison between six conventional algorithms on image in Fig. 3(a): (a) Otsu thresholding; (b) Iterative thresholding; (c) Entropy thresholding; (d) Canny edge detection; (e) Valley edge detection; and (f) Graph based operation (MST).

It is very different to the image in Fig. 3(a), the image in Fig. 5(a) is a complicated multiple rock fracture image, where, several fractures intersected each other, in addition to the rough surface, it has 3D geometry and shadow properties. Iterative thresholding separates the image into different blocks, where, black color presents fractures and spurious fractures (Fig. 5(b)), which is an over-segmentation result; Canny edge detector extracts main fractures roughly, but the resulted image includes many spurious fractures and noise, and the detected fractures are not continuous (Fig. 5(c)); Valley edge detector extracts the centerlines of the fractures, but there are still some noise and gaps (Fig. 5(d)); FCM result is the similar to that by Iterative thresholding, and it consists of some fault blocks; both Clustering analysis and MST detect the most parts of fractures (Fig. 5(f-g)), the results are better than that by FCM, but they also produce some spurious fractures; and the new algorithm can make a good detection result as shown in Fig. 5(h).

Fig. 5.Comparison between new algorithm and six conventional algorithms for a complicated rock fracture image: (a) Original rock fracture image; (b) Iterative thresholding; (c) Canny edge detection; (d) Valley edge detection; (e) FCM operation; (f) Cluster analysis operation; (g) MST operation; and (h) New algorithm operation.

Fig. 6 lists four complex fracture images and their corresponding processing results. The first row presents four original images, and the second row is for the above images’ corresponding centerline extraction results.

Fig. 6.Four complex fracture image processing results: (a), (b), (c) and (d) are for four original fracture images; (e), (f), (g) and (h) present the corresponding centerline extraction results of the above four fracture images.

As shown in Fig. 6(a), it is a pavement image including a crack network, and the cracks cut the image into 12 blocks with 10 junctions. The extraction result is that the boundaries of the blocks are traced well and the 10 junction points are detected; the only shortage is the block boundaries need to extend 5-10% longer for the weak edge parts, which can be resolved by a post function.

The image in Fig. 6(b) presents for rock fractures, it includes a fracture network with shadows and noise, and the fractures are thick. The result is that the shadows and noise are not detected as fractures; and the centerlines are well located in the fractures. The remaining two rock fracture images (Fig.6(c-d)) are more complex, where, surfaces are very rough, illuminations are uneven; in addition to the noise, and there are characteristics of 3D geometry and shadows, which make fracture tracing harder, even though, the new algorithm satisfactorily extracts the rock fracture centerlines.

As above experiments, compared to other widely used image segmentation algorithms, the new algorithm can accurately and well extract centerlines of rock fractures even in a noised image, which is significance for rock fracture measurement and analysis in different rock engineering applications.

 

5. Conclusion

According to the rock fracture characteristics, the paper proposes a novel method for the fracture centerline extraction based on Frangi filter and Steger algorithm. Firstly, Frangi filter is used to obtain the fracture mask. Secondly, based on the mask, Steger algorithm and the segments linking algorithm are applied to extract the fracture centerlines roughly. Finally, to do the centerline tracing accurately, the centerline points are searched along the direction and linked based on the original image information. The paper chooses the pavement crack and rock fracture images for testing. The results verify that the centerline extraction by the new algorithm is more accurate than traditional algorithms.

The centerline extraction on the original fracture images is based on Steger algorithm, and the results are more precise in the studied images. Since the centerline extraction is directly from a gray level image, without image binarization step and any binary image operation step, the whole procedure of centerline extraction is more efficient.

The limitation of the proposed method is that there are some leakages of centerlines on intersections in some cases, which will be studied further.