# An Improved Algorithm in Porosity Characteristics Analysis for Rock and Soil Aggregate

Authors: Yanfang Wu; Xiao Li; Yu Wang; Bo Zhang
Journal: Discrete Dynamics in Nature and Society
 (2014)

---

## Abstract

The rock and soil aggregate (RSA) is a special inhomogeneous multiphase geomaterial. It is crucial for stability of engineering by study of RSA mesodamage characteristic. This paper aims at investigating the porosity evolution characteristics of RSA by X-ray computed tomography (CT) under uniaxial compressive loading. X-ray tomography images were used to extract defects of
RSA specimen under different compressive loading. In this paper, we proposed an improved Ostu method to calibrate the beam hardening phenomenon which is caused by X-ray. Also, based on this Ostu method, the outline of rock blocks in RSA is extracted, and the double gray level threshold of soil and rock block is obtained to ensure the reliability of the porosity calculation. We can conclude that the main reason of RSA cracking is the elasticity mismatch between rock blocks and soil, and the porosity evolution of RSA can be divided into four
typical stages. These results may be useful to reveal the mesoscopic cracking mechanism and establish mesodamage evolution equation and constitutive relation for RSA.


---

## Body

## 1. Introduction

With the development of modern science and technology, digital image processing plays an increasingly important role in medical treatment, scientific research, and military. Digital image processing uses computer algorithms to create, process, communicate, and display images. The major step of image processing is image segmentation, because its inputs and outputs are images. Some useful and attractive characteristics can be obtained by using image segmentation. Image segmentation method has an important influence on computer tomography (CT). Figure1 shows the main steps of CT processing: the first step is detecting the sample by X-ray and obtaining the projection; the second step is selecting computer algorithm to rebuild image; the third step is noise removal; the fourth step is data analysis; and the last step is deriving the result by image segmentation [1, 2].Figure 1
Image processing.For most of image segmentation methods, there are two basic properties of intensity values: discontinuity and similarity [3]. Kapur et al. proposed a method with a sound theoretical foundation and some examples are given on a number of real and artificially generated histograms [4]. White and Rohrer [5] denoted two cost-effective thresholding algorithms for use in extracting binary images of characters from documents. The first nonlinear algorithm is implemented with a minimum of hardware and is intended for many character image extraction applications, and the second is a more aggressive approach directed toward specialized, high-volume applications which justify extra complexity. Prusty and Dash [6] introduced some algorithms for image segmentation, such as iterative thresholding and region growing, implementation of Hough transform, a new image segmentation based on maximum a posteriori. From [3], we know that there are two typical algorithms of image segmentation: one is image segmentation algorithm based on edge detection and the other one is image segmentation algorithm based on thresholding. Compared with the method of edge detection, image segmentation algorithm based on thresholding has a central position of application, according to its intuitive properties, simplicity of implementation, and computational speed. In this paper, we mainly use image segmentation algorithm based on thresholding.Interest in digital image processing methods stems from two principal application areas: improvement of pictorial information for human interpretation and processing of image data for storage, transmission, and representation for autonomous machine perception. The image segmentation also has been used in rock and soil aggregate (RSA), which is formed in the quaternary period as a type of extremely inhomogeneous and loose geomaterial with a certain percentage of rock blocks, composed of rock blocks of various sizes and high strength, fine grained soil and pores [8]. The geometry of the structure of pore space plays a fundamental role in governing fluid transport through porous media. Much effort has, therefore, gone into the problem of understanding the relationship between the geometrical properties of the pore structure and the bulk flow properties of the medium. Establishing the relationship is difficult for two main reasons, the difficulty in making direct measurements on the pore space and the random nature of its structure.With X-ray tubes becoming a standard tool in geophysical laboratories, there are many researches on the three-dimensional distribution of grains and voids in natural rock over the last 10 years [8–13]. Xu et al. [8] compared the nature soil-rock and simulated-rain soil-rock mixtures. Vallejo and Mawby [9] analyzed the stability of slopes made of granular material-clay mixtures. Hirono et al. [10] used X-ray to visualize the advecting fluid image during permeability testing under atmospheric pressure. In addition to natural rock, X-ray computed tomography (CT) as a tool is also applied to analyze the property of micromechanical investigations of cement and mortar [14], geomaterials-soils [15], and the rock mesodamage propagation law [16]. Three-dimensional image analysis is used to measure internal crack growth during each load increment [17], solute transport in cracked concrete [18]. Many scholars have proposed several qualitative and quantitative description methods to “rock-like” material based on X-ray CT, using the CT to analyze the damage and fracturing feature [5, 19]. However, there are few reports about applying CT technology for RSA currently. So in this paper, we will propose an improved Ostu algorithm to obtain some characteristics of RSA CT slice, analyze the porosity extraction of RSA under three typical stress levels, and verify the reliability of the porosity calculation.The paper is organized as follows. A brief outline of the beam-hardening is presented in Section2. In Section 3, we present the greater detail on our improved Ostu algorithm. The porosity characteristic parameter statistics of RSA CT slice is shown in Section 4. Section 5 gives a brief discussion about main results.


## 2. Experimental Device and Beam-Hardening

A high-resolution X-ray CT facility is used to detect the RSA at the Institute of High Energy Physics, Chinese Academy of Sciences. The specifications of CT machine are as follows: the maximum penetration thickness of X-ray is 50 mm Fe, the spatial resolution is 4 Lp/mm, and the minimum focus-object-distance is 0.13 mm. Figure2 shows the 450 KV X-ray computed tomography machine [7].Figure 2
450 kV X-ray computed tomography machine [7].Image a traditional chest X-ray, obtained by placing the subject against an X-ray sensitive plate and “illuminating” the individual with an X-ray beam in the form of a cone. The X-ray plate produces an image whose intensity at a point is proportional to the X-ray energy impinging on that point after it has passed through the subject. This image is the 2D equivalent of the projections. The method now is referred to commonly as the Radon transformation. One of the questions which has been raised in the literature is the effect of heterogeneity of the X-ray beam.Since all substances attenuate low-energy X-rays more strongly than high-energy ones, primarily because of photoelectric absorption, a heterogeneous beam traversing an absorbing medium becomes proportionately richer in high-energy photons and, hence, more penetrating or “harder.” Hence, the attenuation produced by a given material, defined as the negative logarithm of the ratio of transmitted to incident intensities, is not strictly proportional to its thickness. In other words, Beer’s law does not apply to a heterogeneous beam. The result is that the projections are distorted, a projection being a plot of attenuation against beam position as the beam scans linearly across the object.The problem can be approached by resolving the total beam intensityI into its energy spectrum:
(1)I=I0e-μx,
where I0 is the incident energy flux density per unit energy interval, μ is the linear attenuation coefficient for photons of energy, and x is the thickness of absorber.Traditional reconstruction algorithms assume that the X-ray is monochromatic, which means the energy obeys Beer’s law, while in fact X-ray is polychromatic either in industrial or in medical CT, as Figure3 shows. So only polychromatic projection data have been obtained in actual CT system. Once this kind of data is used to reconstruct the image directly, beam-hardening artifacts appear in the reconstructed image, as Figure 4 shows, which contains 512×512 pixels with gray levels from 1 to 255. The main phenomenon for the beam-hardening is that the gray center is smaller than the edge.Figure 3
Relationship between attenuation coefficient and thickness of the absorber.Beam-hardening artifacts. (a) Reconstructed CT image. (b) Distribution of image intensity.
(a)
(b)It is well known that beam-hardening is relative to the material, and the detecting parameters of CT are different with varying materials. Hence, it is difficult to afford a high quality using the same algorithm to avoid the beam-hardening. In the following section, we will introduce our improved Ostu’s method for image thresholding.


## 3. Improved Ostu’s Method

In order to export our improved Ostu’s method, we firstly introduce the traditional Ostu method [20]. Suppose an image is a 2D grayscale intensity function and contains N pixels with gray levels from 1 to L. The number of pixels with gray level i is denoted fi, giving a probability of gray level i in an image of
(2)pi=fiN,p1(k)=∑i=1k‍pi.In the case of bilevel thresholding of an image, the pixels are divided into two classes,C1 with gray levels [1,2,…,t] and C2 with gray levels [t+1,…,L]. Then, the gray level probability distributions for the two classes are
(3)C1:p1ω1(t),p2ω1(t),…,ptω1(t),C2:pt+1ω2(t),pt+2ω2(t),…,pLω2(t),
where
(4)ω1(t)=∑i=1t‍pi,ω2(t)=∑i=t+1L‍pi.
Also, the means for classes C1 and C2 are
(5)μ1=∑i=1tipiω1(t),μ2=∑i=t+1Lipiω2(t).
Let μT be the mean intensity for the whole image. It is easy to show that
(6)ω1μ1+ω2μ2=μT,ω1+ω2=1.Using discriminant analysis, Ostu defined the between-class variance of the thresholded image [15]:
(7)σB2=ω1(μ1-μT)2+ω2(μ2-μT)2.For bilevel thresholding, Ostu verified that the optimal thresholdt* is chosen so that the between-class variance σB2 is maximized; that is,
(8)t*=Argmax⁡{σB2},1≤t≤L.For general image, the traditional Ostu method can be easily used to obtain thresholding of this image. But for CT image with the phenomenon of beam-hardening, it is hard to afford the thresholding image. So we need to introduce the improved Ostu’s method to solve this problem. Compared with the traditional Ostu’s method, the main change of the improved Ostu’s method is that the thresholding method can be extended to dual thresholding, which is combined with Poisson distribution. Improved dual thresholding generalizes to the expression:(9)P(fj)=(∑i=1LIpi)fje(-∑i=1LIpi)fj!,
where P(fj) is subject to the Poisson probability distribution. The same as with (7), we can give the following equation:
(10)σB2=ϖ1(μ1-μT)2+ϖ2(μ2-μT)2+ϖ3(μ3-μT)2,
where
(11)ϖ1=∑i=1k1pi,ϖ2=∑i=k1+1k2pi,ϖ3=∑i=k2+1Lpi,(12)μ1=1ϖ1∑i=1k1ipi,μ2=1ϖ2∑i=k1+1k2ipi,μ3=1ϖ3∑i=k2+1Lipi.As in (5) and (6), we can obtain that the following relationships hold:
(13)ϖ1μ1+ϖ2μ2+ϖ3μ3=μT,ϖ1+ϖ2+ϖ3=1.We see that theϖ and μ terms, there for σB2, are functions of k1 and k2. The two optimum threshold values, k1* and k2*, are the values that maximize σB2(k1,k2). In other words, as in the single-threshold case discussed in (6), we find the optimum thresholds by finding
(14)σB2(k1*,k2*)=max⁡0<k1<k2<L-1σB2(k1,k2).The procedure starts by selecting the first value ofk1. Next, k2 is incremented through all its values greater than k1 but less than L-1 (i.e., k2=k1+1,…,L-2). Then, k1 is incremented to its next value and k2 is incremented again through all its values greater than k1; this procedure is repeated until k1=L-3. The result of this process is a 2D array, σB2(k1,k2), and the last step is to look for the maximum value in this array. The values of k1 and k2 corresponding to that maximum are the optimum thresholds, k1* and k2*. The thresholded image is then given by
(15)g(x,y)={a,f(x,y)≤k1*,b,k1*≤f(x,y)≤k2*,c,f(x,y)≥k2*,
where a, b, and c are any three valid intensity values. As shown in Figure 5, the red circle in the left image represents the smaller thresholding, as the k1* in (15), and the red circle in the right image represents the larger thresholding, as the k2* in (15).Dual thresholding.
(a)
(b)The improved between-class variance is summarized in the following algorithm, wheref(x,y) is input image:(1)
compute an CT image;(2)
compute the Poisson probability distribution, using (9);(3)
compute the cumulative sums ofϖi, for k=1,2,…,L, using (11);(4)
compute the cumulative means ofμi, for k=1,2,…,L, using (12);(5)
compute the global intensity mean ofμT, using (13);(6)
compute the between-class variance,σB2(k1,k2), for k=1,2,…,L, using (14);(7)
obtain the improved Ostu threshold, as the value ofk1* and k2*, for σB2(k1*,k2*) which is maximum.For evaluation of the performance of our proposed method, the RSA specimen is detected by CT method, which is mentioned in the previous section. According to the geotechnical test, technical manuals, and soil specimen preparation standard BS1377 (1990), RSA is cylinder shaped (50 mm diameter and 100 mm height) and the diameter of the blocks should not be greater than 10 mm for the test [7].The following figures can illustrate our improved algorithm by CT slice of RSA. Figure6 shows the image segmentation of CT slice under three stress levels with different algorithms. The first column of Figure 6 (also see Figure 6  (1), (5), and (9)) shows the same CT slice under different stress strain. The second column of Figure 6 (also see Figure 6  (2), (6), and (10)) shows the results obtained by Ostu algorithm; it is easy to see that image segmentation has failed. The reason for the failure can be traced to the fact that there exists the beam-hardening of the boundary between object and background of RSA. Comparing with the second column of Figure 6, in the third column, we add to the method of smoothing histogram under the basis of Ostu algorithm. For the effect with the method of smoothing histogram which can be shown in Figure 7, Figure 7(a) shows the histogram of a CT image, and Figure 7(b) shows the smoothing histogram of a CT image. Comparing with Figures 7(a) and 7(b), it is also easy to see that the improvement in the shape of the histogram due to smoothing is evident. Through adding to the method of smoothing historgam, we want the thresholding of smoothing image to be successful, but in the third column of Figure 6 (also see Figure 6  (3), (7), and (11)), it is easy to see that image segmentation has also failed and ring artifact also exits. The fourth column of Figure 6 (also see Figure 6  (4), (8), (12)) shows the results obtained by our improved Ostu algorithm; we can see that the boundary of rock can be recognized obviously. So we can conclude that the segmentation is highly successful, and our improved Ostu algorithm has an important influence on statistics porosity characteristic parameter.Figure 6
Porosity extraction of RSA under three typical stress levels using improved Ostu’s method (first line is 0 MPa, Second line is 3.87 MPa, Third line is 6.33 MPa).Histogram and smoothing histogram of CT slice.
(a)
(b)


## 4. The Discussion and Analysis of Porosity Characteristic Parameter Statistics

The excellent quality of this thresholding segmentation with RSA image is that it can be applied to the changing stress, and the ideal segmentation result is to get a complete rock contour edge. In the following section, we will give the porosity statistical of the thresholding segmentation with RSA image, which includes each CT slice and the whole RSA specimen under five strain-stresses. The CT slice-porosity-stress strain with different stress-strain is shown in Figure8.Figure 8
Porosity variation along the 50 slices of the RSA specimen (rock percentage 40%) under five successive stages of compressive loading.From Figure8, we can see that when the strain-stress is 0.7 MPa, the porosity of all CT slices is in decline, which shows the RSA is in crack closure and soil consolidation stage; this phenomenon continues to 2.74 MPa. When the stress reaches 2.74 MPa, the porosity of every layer goes down to a minimum; as the load continues to increase, the porosity of the sample increases suddenly, and the cracks in sample gradually appears, which shows the specimens are in linear elastic deformation stage. When the strain-stress reaches 3.87 MPa, RSA will be in the rock-soil jointed cracks with initiation stage and the stable crack growth stage, and porosity increases quickly. When the strain-stress reaches 6.33 Mpa, the RSA is in the unstable propagation stage; there are a large number of cracks inside the sample. More details are shown in Table 1.Table 1
The porosity calculation for typical CT slice under different stress level.
Axial stress
Slice-1
Slice-5
Slice-13
Slice-17
Slice-25
Slice-31
(MPa)
(%)
(%)
(%)
(%)
(%)
(%)
0
9.78445
8.66302
8.48999
8.66955
11.93590
10.21865
0.7
5.95077
4.03601
3.96500
4.38860
14.84958
9.21467
2.74
7.20932
4.54531
4.82852
4.91912
8.52500
7.24196
3.87
23.29622
9.83196
7.36521
9.81701
22.90364
18.02289
5.26
24.50009
9.93480
7.81982
10.60299
24.95878
19.60627
6.33
27.07269
10.14701
8.83923
12.39614
29.69426
23.39253It is necessary to analyze the specimens including all CT slices as shown in Figure9, which shows the changes of porosity under different stages. Every stage of the mesoscopic damage extension situations is as follows:(1)
1-2 stage is corresponding to crack closure and soil consolidation stage, and the porosity decreases from 9.55% to 6.55%;(2)
2-3 stage is corresponding to linear elastic deformation stage, and the porosity is 6.32%. After the closure of the microcrack, the stress-strain curve exhibits a linear elastic relationship;(3)
3-4 stage is corresponding to rock-soil jointed cracks initiation stage, and the porosity is 14.29%. The growth of these cracks has been shown to occur around the rock blocks. The opening of cracks with faces around the block surface is detected as a departure from the linear lateral and volumetric strain behavior;(4)
4–6 stage is corresponding to stable crack growth stage, and the porosity is 14.29%; local deformation is in the rapid development and then destruction.Figure 9
Relationship between RSA porosity and axial stress.


## 5. Conclusion

With the help of X-ray computed tomography, we discuss the internal structure of the rock mixture specimens under loading. CT technology can not only interpret accurately details of the data inside the object, but also give the quantitative data of radiation absorbed. In this paper, in order to afford the thresholding image with the phenomenon of beam-hardening, an improved Ostu method has been developed. Using our improved Ostu method, we can obtain image segmentation of the rock in RSA. Compared with other methods, our improved Ostu method can help us to understand and research the cracking characteristics of RSA under uniaxial compression test real-time. We conclude that the statistical method of image segmentation algorithm based on mathematical morphology can extract the crack porosity and other characteristics. With increasing of strain-stress, the porosity of RSA increases first, and then it is reduced. In addition, there exist four stages of the mesoscopic damage extension situations during the process. Thus, for beam-hardening CT image, our new algorithm is a significant improvement. It is also worth mentioning that our algorithm can be applied to obtain the crack in RSA after segmentation of rock.



---

*Source: 798235-2014-08-27.xml*