# Fracture Segmentation Method Based on Contour Evolution and Gradient Direction Consistency in Sequence of Coal Rock CT Images

Authors: Zhiwei Li; Guoying Zhang
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The coal rock exhibits obvious heterogeneity and anisotropy after a long-term geological evolution in nature. Therefore, the coal rock CT (computed tomography) image shows uneven grey scale, low contrast, and the fractures with weak boundaries. The accurate segmentation of the coal rock fracture networks is challenging. In this paper, a segmentation method of fractures based on contour evolution and gradient direction consistency is proposed to accurately segment the fracture networks in the sequence of coal rock CT images. According to the contour variation rule of the fractures in the discrete 3D (three dimensional) space formed by the sequence of CT images, the fracture contour evolution model (FCEM) is constructed and the preliminary segmentation results of fractures are obtained from FCEM. A 3D adaptive median filtering (3DAMF) and a 3D bilateral filtering (3DBF) are proposed. The high density miscellaneous point noises in the coal rock CT images are filtered by the 3DAMF. And the boundaries of fractures are enhanced by 3DBF. According to the similarity of the preliminary segmentation results of fractures and the real contours of fractures, the preliminary segmentation results of fractures are optimized based on the gradient direction consistency model (GDCM) proposed in this paper to obtain the accurate boundaries of fractures. The fracture segmentation method proposed in this paper can obtain accurate boundaries of fractures with weak boundaries, and the experimental results show that the segmentation efficiency for sequence is high and adaptability is strong.


---

## Body

## 1. Introduction

With the increase in the depth of the underground coal mining, the gas content in the coal rock increases accordingly. The disasters such as gas outbursts and gas explosions have caused huge economic losses and a large number of casualties. The fractures of the coal rock are the storage place and the migration channels of the gas. The coal mining will cause changes in the stress field of the coal rock, which in turn will affect the extension, expansion, and redistribution of the fractures in the coal rock. Among them, the fractures with weak boundaries often indicate the most likely change position of the coal rock under mining stress. The accurate segmentation of the fractures in the coal rock is conducive to obtaining the shape parameters of the fractures (such as length, area, azimuth, etc.) and then can accurately give a quantitative description of the fractures in the coal rock. This is of great significance for mastering the evolution rule of the fractures in the coal rock under mining stress, as well as the coal rock seepage calculation and the damage calculation.After many stages of the complex natural tectonic movement, the coal rock exhibits obvious heterogeneity and anisotropy, and the complex fracture networks were formed inside. Detecting the internal structure of the coal rock by industrial CT (computed tomography) is a nondestructive testing technology. Compared with other methods, CT scanning can detect the structure of fractures in the coal rock at the microscopic scale (micron or nanometre) and has the advantages of intuitive imaging, high resolution, and no limitation on the structure of the object to be detected. At the same time, a sequence of CT images contains thousands of images. The coal rock CT images are characterized by uneven grey scale and low contrast and contain high density miscellaneous point noises. Moreover, the fractures with weak boundaries are widely present in the coal rock, which makes it difficult to accurately segment the fractures in the coal rock.For many years, the researches on segmentation methods of the fractures mainly have focused on the 2D image processing techniques based on threshold, edge detection, region growth, and texture analysis [1–9]. H Oliveira et al. [1] used dynamic threshold method and entropy to segment fractures of the road image. A Landstrom et al. [3] extracted longitudinal fractures from the image of road surface using morphological methods and logistic regression models. L Liu. et al. [5] studied the segmentation of fractures in CT image based on wavelet transform and CV model. Y. Fujita and Y. Hamamoto [7] applied the eigenvalues and eigenvectors of the Hessian matrix to discriminate the shape to distinguish fractures. The above methods only utilize the 2D features of the fractures, and these methods are sensitive to noises, uneven grey scale, and low contrast. Due to the lack of more reliable information of fractures, the segmentation results of the fractures with weak boundaries have serious misjudgement and misdetection by these methods.In recent years, the segmentation methods of fractures based on 3D features have been widely concerned [10–12]. Since these methods can utilize the 3D grey scale and 3D structural features, a higher recognition rate and a lower misjudgement rate are obtained. But the corresponding 3D processing technologies of the fractures are rarely reported. MR. Jahanshahi [10] used depth information to reconstruct the image and then identify the fractures. M Voorn et al. [11] studied the fractures of CT image based on Hessian matrix multiscale planar filtering.At present, the researches on the segmentation of the fractures in the coal rock are rarely reported. In this paper, by studying the 3D grey scale and 3D structure features of fractures in sequence of the coal rock CT images, a fracture segmentation method based on contour evolution and gradient direction consistency is proposed. The fracture contour evolution model (FCEM) is proposed based on the 3D variation rule of fracture contours, and the preliminary segmentation results of fractures are obtained by FCEM. According to the characteristics of the coal rock CT images, the 3D adaptive median filtering (3DAMF) is used to reduce the image noises, and the 3D bilateral filtering (3DBF) is used to enhance the boundaries of fractures. At the same time, the gradient direction consistency model (GDCM) is proposed and is used to optimize the preliminary segmentation results of fractures, and finally the accurate segmentation results of fractures are obtained.


## 2. Fracture Contour Evolution Model

The fractures are generated from natural evolution, and their features (such as shape, length, area, and azimuth) are random. Therefore, the structural features of the fractures are rarely mentioned in the studies of fracture segmentation. Most of the existing methods are based on the grey scale feature of the fractures. A discrete 3D spaceT(x,y,z) is formed by the sequence of coal rock CT images, where (x,y) represents the 2D coordinates of the pixel points in the single CT image and Z represents the serial number of CT image. Discriminating whether a target in the image is a fracture, one often makes a judgment by analyzing the same position in the adjacent images. That is to say, in the 3D space, the fracture contours in the adjacent CT images are very similar, and the smaller the interlayer spacing d is, the higher the contour similarity is. This characteristic of the fracture structure also clearly exists in the fractures with weak boundaries. That is, the fractures in the z-th image can be obtained by the slight shape deformation of the fractures among the adjacent CT images. Therefore, this paper proposes FCEM for obtaining the fractures in the CT image.n images of consecutive serial numbers are selected and added to the set Iz+1,Iz+2,⋯,Iz+n. The images in the set Iz+1,Iz+2,⋯,Iz+n require as much as possible including the complete fracture networks. The fracture networks in the set Iz+1,Iz+2,⋯,Iz+n are segmented and added to the set Cz+1,Cz+2,⋯,Cz+n. In this paper, the evolution rule is based on the curve evolution theory [13]. This theory mainly studies the process of curve evolution between two curves with time by the normal. During the process of curve evolution, only the point motion in its normal direction is considered.Firstly,Cz+1 is discretized into a set Ωz+1=[p0z+1,p1z+1,⋯,piz+1] [14]. At the same time, the normal lines of each point in the set Ωz+1 are available, and these normal lines are added to the set Nz+1=[N0z+1,N1z+1,⋯,Niz+1]. As shown in Figure 1, each normal line in Nz+1 and the contour curve Cz+2 intersect into the set Ωz+2=[p0z+2,p1z+2,⋯,piz+2]. In the same process, the discrete points set Ωz+n and the normal line set Nz+n on the curve Cz+n can be obtained.Figure 1
The process of contour evolutions.As shown in Figure1, each point evolves along the normal direction and reaches the fracture contour of the adjacent image. The motion time of the point is related to d. The speed v(piz+n) of each point in the set Ωz+n is different. v(piz+n) is obtained mainly by simulating the evolution process of piz+1→piz+2→,⋯,→piz+n. The evolutionary process of the points can be seen as a uniform acceleration motion. Set the speed v(piz+1) of the point piz+1 to zero. The acceleration formula for the point piz+2 is calculated as follows:(1)apiz+1,piz+2=2·Dispiz+1,piz+2-vpiz+1·dd2,where Dis(piz+1,piz+2) is the Euclidean distance between the point piz+1 and the point piz+2. The formula of the speed for the point piz+2 is as follows:(2)vpiz+2=vpiz+1+apiz+1,piz+2·d.Similarly, the accelerationa(piz+n-1,piz+n) and the speed v(piz+n) of the point piz+n are calculated. Therefore, the distances of the points in the set Ωz+n moving along the normal direction are calculated by formula (3).(3)Dispiz+n,p~iz+n+1=vpiz+n·d+12·apiz+n-1,piz+n·d2.Finally a new setΩ~z+n+1=[p~0z+n+1,p~1z+n+1,⋯,p~iz+n+1] is obtained. The points in the set Ω~z+n+1 are fitted to a new curve C~z+n+1 [15, 16]. The normal line set of the Ω~z+n+1 is N~z+n+1=[N~0z+n+1,N~1z+n+1,⋯,N~iz+n+1], and the corresponding normal direction angle set is θ~z+n+1=[θ~0z+n+1,θ~1z+n+1,⋯,θ~iz+n+1]. C~z+n+1 is very similar to the exact boundary Cz+n+1 of the fractures to be segmented, so C~z+n+1 is taken as the preliminary segmentation result of the image Iz+n+1.


## 3. 3D Enhancement Methods of Coal Rock CT Image

The FCEM only utilizes the contour variation characteristic of the fractures in 3D space. In this paper, the resultC~z+n+1 obtained by FCEM is optimized according to the 3D grey scale feature of the fractures. At the same time, for the characteristics of the coal rock CT images with low contrast, uneven grey scale, and high density miscellaneous point noises, this paper proposed 3DAMF and 3DBF, which can reduce noises and enhance the boundaries of fractures.The coal rock is a natural body, and the coal rock CT image shows an obvious low contrast, as shown in Figure2.Figure 2
The images of grey scale distribution statistics.
(a) (b)The red line in Figure2(a) is located on a coal rock fracture with obvious fracture characteristics. The green line is located on the coal rock background and spans two fractures with weak boundaries. The curves of red and green in Figure 2(b), respectively, show grey scale value distribution of each point on the same color line in Figure 2(a). It can be seen from Figure 2(b) that the grey scale distribution of the weak boundaries is basically consistent with the coal rock background. Therefore, this paper uses the local histogram statistical method to perform contrast stretching on the coal rock CT images.The coal rock image is a heterogeneous body with high density miscellaneous point noises. According to the median filtering principle, when the size of the filtering window is small, the image details can be well protected, but the filtering effect of noises is not very good. Moreover, if the number of the noise points in the filtering window is larger than the number of the coal rock pixels in the entire window, the noises cannot be well filtered. In order to protect the fractures with weak boundaries, this paper proposes 3DAMF. This method considers all information of the pixels in the 3D neighbour space of each pixel. So the size of the filtering window becomes a cube (e.g., 3∗3∗3), and the size of the filtering window can be adaptively adjusted according to the local information. The filtering method can protect the information of weak boundaries while filtering high density miscellaneous point noises, and the principle is as follows:( 1 ) To Determine whether the median value in the filtering window is a noise point, and if so, increasing the size of window to find a nonnoise median output.( 2 ) To determine whether the point to be processed is a noise point. If not, do not process; otherwise replace the point value with the median value of step (1).For the feature of uneven grey scale in the coal rock CT images, the bilateral filtering is used to enhance the images. The bilateral filtering is a nonlinear filter composed of a spatial domain filter and a range filter. This method not only considers the distance between pixels, but also considers the similarity of grey scale in the pixel range domain. This method has characteristics of being simple, noniterative, and local. In order to improve the effect of bilateral filtering, the 3D domain information of the pixel to be processed is involved in the process of calculation, and the size of template becomes a cube (such as 3∗3∗3). The 3D range filter of this method is expressed as(4)di,j,kl,m,n=exp⁡-i-l2+j-m2+k-n22σd2,where (i,j,k) is the coordinate of the central position within the template window and (l,m,n) is the coordinates of the other positions within the template window. σd is the distance standard deviation of the Gaussian function. The 3D spatial domain filter is expressed as(5)ri,j,kl,m,n=exp⁡fi,j,k-fl,m,n2σr2,where f(x,y,z) represents the pixel value of coordinate (x,y,z) and σr is the grey standard deviation of the Gaussian function. Then the 3DBF is expressed as(6)gi,j,k=∑∑∑fi,j,k·di,j,kl,m,n·ri,j,kl,m,n∑∑∑di,j,kl,m,n·ri,j,kl,m,n.After being processed by the 3DBF,g(i,j,k) is the pixel value of the central point in template window. The 3D enhancement result of Figure 2(a) is shown in Figure 3.Figure 3
The 3D enhancement result.


## 4. Gradient Direction Consistency Model

According to the 3D grey scale feature of the fractures, the points in the setΩ~z+n+1 continue to evolve into a new point set Ωz+n+1=[p0z+n+1,p1z+n+1,⋯,piz+n+1]. The point p~iz+n+1 in the set Ω~z+n+1 will always be evolved into the point piz+n+1 in Ωz+n+1. According to the curve evolution theory, the point piz+n+1 is located on the normal line of the point p~iz+n+1. Since the points piz+n+1 and p~iz+n+1 are, respectively, located on two adjacent curves with similar shapes, the gradient direction of the point piz+n+1 and the normal direction of the point p~iz+n+1 are substantially identical. A segmentation method of fractures based on GDCM is proposed. This method finds a closest point piz+n+1 on the normal line of the point p~iz+n+1, and the gradient direction of the point piz+n+1 is similar to the normal direction of the point p~iz+n+1. In the fractures with weak boundaries, although the gradient value is small, the characteristic of gradient direction consistency is unchanged. Therefore, this method is effective for fractures with weak boundaries.There is a setΩiz+n+1=[pjΔd⋯p2Δd,pΔd,p0,p-Δd,p-2Δd,⋯p-jΔd] on the normal line N~iz+n+1, where the point p0 and the point p~iz+n+1 are the same point, and Δd is the size of step when the discrete points are taken on the normal line N~iz+n+1, as shown in Figure 4.Figure 4
The process of gradient direction consistency.Find a point inΩiz+n+1 that is as close as possible to the point p~iz+n+1, and the gradient direction of the point is close to the direction θ~iz+n+1 of the normal N~iz+n+1. In this paper, the gradient values and gradient directions of the points in the set Ωiz+n+1 are calculated using the Sobel operator, and the convolution kernels are(7)-10+1-20+2-10+1-1-2-1000+1+2+1.For the pixel informationI of the point domain to be detected, I is convoluted with two convolution kernels, respectively, to obtain gradient values Gx and Gy in the horizontal and vertical directions as follows(8)Gx=-10+1-20+2-10+1∗IGy=-1-2-1000+1+2+1∗I.The gradient valueG of this point is calculated by Gx and Gy:(9)G=Gx2+Gy2.The gradient directionθ of this point can also be calculated.(10)θ=arctan⁡GyGx.The mathematical expression of GDCM is(11)EGDCM=λ1·θ-θ~iz+n+1+λ2·Disp±jΔd,p~iz+n+1+λ3·1G,where λ1, λ2, and λ3 are the weights of each influencing factor, and Dis(pi′z+n+1,p~iz+n+1) is the distance from the point p±jΔd to the point p~iz+n+1. When EGDCM is the smallest, the point p±jΔd is the point that the point p~iz+n+1 is finally reaching along the evolutionary direction. This point p±jΔd is added into the set Ωz+n+1. Finally, the points in the set Ωz+n+1 are fitted to form a curve Cz+n+1 which is the final segmentation result of the fracture.


## 5. Experimental Results and Analysis

In order to verify the effectiveness and superiority of the proposed segmentation method, the X-Ray (ACTIS300-320/225) industrial CT detection system of the State Key Laboratory of Coal Resources and Safe Mining of China University of Mining and Technology (Beijing) was used to carry out CT real-time scanning experiments on different stages of coal rock triaxial loading. The coal samples used in the experiment were taken from the Liu-pan-shui mining area in Guizhou Province. In this experiment, CT images of the coal rock under different loading conditions were obtained, and the CT scan scale was 50 microns. A sequence of CT images with relatively developed fractures was selected, and this sequence contained 876 CT images of the size680×680.As shown in Figure2(a), the coal rock image contains an unrelated peripheral cavity. The coal rock images were processed by the Gaussian blur method and binarization method, and then the search method of contours is used to determine the coal rock region. Before the segmentation of the fractures, the peripheral cavity is first removed, so as to avoid the influence of unrelated information of the coal rock CT image on the segmentation of fractures.Figure5 is a set of two CT images with the peripheral cavity removed and located at different locations in the coal rock column.Figure 5
The coal rock CT images after cavity removed.
(a) (b)The sequence of the coal rock images was processed by the traditional image processing methods. The median filtering [17] and bilateral filtering [18] methods were used to denoise and enhance the images, and the adaptive threshold method [19] was used for binarization. At the same time, morphological methods were used to reduce the connection of noises and fractures. The processed results are shown in Figure 6.Figure 6
The results of the traditional image processing method.
(a) (b)The processed results of Figure6 contain a large number of nonfracture noises, and the loss and discontinuity of the fractures in the results of segmentation are serious. Due to the use of the morphological methods, the thickness of the fractures obtained does not match the actual fractures.In this paper, the method in [11] was used to segment the fractures in the coal rock CT images. This method based on Hessian matrix multiscale planar filtering is used to segment the fractures from a sequence of CT images, and the main parameters used are as shown in Table 1. The results are shown in Figure 7.Table 1
Parameter settings.
Parameter a v g m a t c o n s t h r e s h m a x m a t s m i n s m a x s s t e p Value 43 31 87 2 5 1Figure 7
The results of multiscale Hessian fracture filtering.
(a) (b)As described in [11], the multiscale Hessian fracture filtering has a better effect on the segmentation of fractures in homogeneous materials, and the noises have a great influence on the segmentation of fractures with weak boundaries. Therefore, the fractures with weak boundaries in the rectangular mark area of Figure 7 are not recognized or incompletely segmented. At the same time, since the intersection position of multiple fractures is not a planar feature, this method cannot segment the fractures in the circular mark area in Figure 7(b). The Gaussian scale and the grey scale threshold in this method have a large influence on the noises, and the less the noise that is retained, the more the weak boundaries are lost. So the appropriate values of these parameters are difficult.The processing results obtained by the method of this paper are shown in Figure8, and the main parameters used are as shown in Table 2.Table 2
Parameter settings.
Parameter d ( u m ) λ 1 λ 2 λ 3 σ d σ r Δ d Value 50 0.6 0.2 0.2 5 10 3Figure 8
The results of the method in this paper.
(a) (b)In the experiment, firstly, five continuous CT images of the complete fracture networks were selected from the sequence, and the corresponding fracture networks were obtained by manual segmentation. The images on both sides of the five serial numbers were iteratively processed using the method proposed in this paper. Therefore all results of segmentation are based on five a priori fracture networks. By the principle of FCEM and GDCM, the segmentation results of fractures are only the extension or contraction of the a priori fracture networks, so the fractures segmented do not appear at other positions of the image. That is to say, a small amount of noises in the results of segmentation will only be presented by the burrs in boundaries of the fractures. At the same time, a fracture that suddenly appears in an image area other than the a priori fracture networks cannot be segmented. Therefore, the a prior fracture networks are required to include all the fractures in the sequence as much as possible, or the sequence is divided into multiple subsequences with similar fracture networks for processing. Because FCEM and GDCM have strong applicability to fractures with weak boundaries and 3DAMF and 3DBF have strong ability to denoise and enhance images, the accurate results of segmentation were obtained in Figure8.


## 6. Conclusions

The segmentation of fractures in the coal rock CT images is very important for other related basic researches. In this paper, a segmentation method of fractures based on contour evolution and gradient direction consistency is proposed, which provides an effective method to segment fractures from the sequence of coal rock CT images. By studying the contour change rule of fractures in 3D space, FCEM is proposed, and the preliminary segmentation results of fractures are obtained. Aiming at the characteristics of uneven grey scale, low contrast, and high density miscellaneous point noises in the coal rock CT images, 3DAMF and 3DBF methods are proposed to improve image quality. Finally, the GDCM based on the grey scale feature of the fractures is proposed to optimize the preliminary segmentation results, and the accurate segmentation results of fractures are obtained. The proposed method has strong antinoise ability and exhibits strong segmentation ability for fractures with weak boundaries. For a fracture that may suddenly appear, the sequence needs to be divided into multiple subsequences with similar fracture networks for processing.



---

*Source: 2980747-2019-02-04.xml*