Discrete Dynamics in Nature and Society

Volume 2014 (2014), Article ID 798235, 7 pages

http://dx.doi.org/10.1155/2014/798235

Research Article

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

Yanfang Wu,¹ Xiao Li,¹ Yu Wang,1,2 and Bo Zhang1,2

¹ Key Laboratory of Shale Gas and Geomechanics, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China

²Graduate University of Chinese Academy of Sciences, Beijing 100049, China

Received 1 July 2014; Accepted 24 July 2014; Published 27 August 2014

Academic Editor: Li Li

Copyright © 2014 Yanfang Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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). Figure 1 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 Section 2. 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. Figure 2 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 intensity into its energy spectrum:

where is the incident energy flux density per unit energy interval, is the linear attenuation coefficient for photons of energy, and 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 Figure 3 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 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.

(a)

(b)

Figure 4: Beam-hardening artifacts. (a) Reconstructed CT image. (b) Distribution of image intensity.

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 pixels with gray levels from 1 to . The number of pixels with gray level is denoted , giving a probability of gray level in an image of

In the case of bilevel thresholding of an image, the pixels are divided into two classes, with gray levels and with gray levels . Then, the gray level probability distributions for the two classes are

where

Also, the means for classes and are

Let be the mean intensity for the whole image. It is easy to show that

Using discriminant analysis, Ostu defined the between-class variance of the thresholded image [15]:

For bilevel thresholding, Ostu verified that the optimal threshold is chosen so that the between-class variance is maximized; that is,

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:

where is subject to the Poisson probability distribution. The same as with (7), we can give the following equation:

where

As in (5) and (6), we can obtain that the following relationships hold:

We see that the and terms, there for , are functions of and . The two optimum threshold values, and , are the values that maximize . In other words, as in the single-threshold case discussed in (6), we find the optimum thresholds by finding

The procedure starts by selecting the first value of . Next, is incremented through all its values greater than but less than (i.e., ). Then, is incremented to its next value and is incremented again through all its values greater than ; this procedure is repeated until . The result of this process is a 2D array, , and the last step is to look for the maximum value in this array. The values of and corresponding to that maximum are the optimum thresholds, and . The thresholded image is then given by

where , , and are any three valid intensity values. As shown in Figure 5, the red circle in the left image represents the smaller thresholding, as the in (15), and the red circle in the right image represents the larger thresholding, as the in (15).

(a)

(b)

Figure 5: Dual thresholding.

The improved between-class variance is summarized in the following algorithm, where is input image:

(1)compute an CT image;(2)compute the Poisson probability distribution, using (9);(3)compute the cumulative sums of , for , using (11);(4)compute the cumulative means of , for , using (12);(5)compute the global intensity mean of , using (13);(6)compute the between-class variance, , for , using (14);(7)obtain the improved Ostu threshold, as the value of and , for 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. Figure 6 shows the image segmentation of CT slice under three stress levels with different algorithms. The first column of Figure 6 (also see Figure 6  , , and ) shows the same CT slice under different stress strain. The second column of Figure 6 (also see Figure 6  , , and ) 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  , , and ), 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  , , ) 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).

(a)

(b)

Figure 7: Histogram and smoothing histogram of CT slice.

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 Figure 8.

Figure 8: Porosity variation along the 50 slices of the RSA specimen (rock percentage 40%) under five successive stages of compressive loading.

From Figure 8, 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.39253

It is necessary to analyze the specimens including all CT slices as shown in Figure 9, 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 41227901) and the Strategic Priority Research Program of the Chinese Academy of Sciences (Grants no. XDB10030000, XDB10030300, and XDB10050400).

References

W. H. Oldendorf, “Isolated flying spot detection of radiodensity discontinuities displaying the internal structural pattern of a complex object,” IRE Transactions on Bio-Medical Electronics, vol. 8, pp. 68–72, 1961.

G. T. Herman, Image Reconstruction from Projection, The Fundamentals of Computerized Tomography, Academic Press, New York, NY, USA, 1980.

R. C. Gonzalez and R. E. Woods, Digital Image Processing, Publishing House of Electronics Industry, 3rd edition, 2010.

J. N. Kapur, P. K. Sahoo, and A. K. C. Wong, “A new method for gray-level picture thresholding using the entropy of the histogram.,” Computer Vision, Graphics, and Image Processing, vol. 29, no. 3, pp. 273–285, 1985.

J. M. White and G. D. Rohrer, “Image thresholding for optical character recognition and other applications requiring character image extraction,” IBM Journal of Research and Development, vol. 27, no. 4, pp. 400–411, 1983.

S. J. Prusty and R. Dash, Study of image segmentation schemes and their applications to microscopic images [Bachelor's thesis], NIT, Rourkela, India, 2009.

Y. Wang, X. Li, B. Zhang, and Y. F. Wu, “Meso-damage cracking characteristics analysis for rock and soil aggregate with CT test,” Science China Technological Sciences, vol. 57, pp. 1–11, 2014.

W. J. Xu, R. L. Hu, and R. J. Tan, “Some geomechanical properties of soil-rock mixtures in the Hutiao Gorge area, China,” Geotechnique, vol. 57, no. 3, pp. 255–264, 2007.

L. E. Vallejo and R. Mawby, “Porosity influence on the shear strength of granular material-clay mixtures,” Engineering Geology, vol. 58, no. 2, pp. 125–136, 2000.

T. Hirono, M. Takahashi, and S. Nakashima, “In situ visualization of fluid flow image within deformed rock by X-ray CT,” Engineering Geology, vol. 70, no. 1-2, pp. 37–46, 2003.

R. A. Ketcham, “Three-dimensional grain fabric measurements using high-resolution X-ray computed tomography,” Journal of Structural Geology, vol. 27, no. 7, pp. 1217–1228, 2005.

T. Ohtani, T. Nakano, Y. Nakashima, and H. Muraoka, “Three-dimensional shape analysis of miarolitic cavities and enclaves in the Kakkonda granite by X-ray computed tomography,” Journal of Structural Geology, vol. 23, no. 11, pp. 1741–1751, 2001.

A. Vervoort, M. Wevers, R. Swennen, et al., “Recent advantages of X-ray CT and its applications for rock material,” in X-ray CT for Geomaterials, J. Otani and Y. Obara, Eds., pp. 79–91, A.A. Balkema, Kumamoto, Japan, 2003.

E. N. Landis, “X-ray tomography as a tool for micromechanical investigations of cement and mortar,” in Advances in X-Ray Tomography for Geomaterials, J. Desrues, G. Viggiani, and P. Be'suelle, Eds., pp. 79–93, Antony Rowe Ltd, Chippenham, UK, 2006.

J. Otani and Y. Obara, “X-ray CT for geomaterials-soils, concrete, rocks,” in Proceedings of the International Workshop on X-Ray CT for Geomaterials, pp. 227–287, A.A. Balkema Publishers, Kumamoto, Japan, November 2003.

X. Ge, J. Ren, Y. Pu, W. Ma, and Y. Zhu, “Real-in time CT test of the rock meso-damage propagation law,” Science in China, Series B: Chemistry, vol. 44, no. 3, pp. 328–336, 2001.

E. N. Landis and E. N. Nagy, “Three-dimensional work of fracture for mortar in compression,” Engineering Fracture Mechanics, vol. 65, no. 2-3, pp. 223–234, 2000.

I. S. Darma, T. Sugiyama, and M. A. B. Promentilla, “Application of X-ray CT to study diffusivity in cracked concrete through the observation of tracer transport,” Journal of Advanced Concrete Technology, vol. 11, pp. 223–234, 2013.

N. Otsu, “A threshold selection method from gray-level histogram,” IEEE Transactions on System Man Cybernetics, vol. 9, no. 1, pp. 62–66, 1979.

P.-S. Liao, T.-S. Chen, and P.-C. Chung, “A fast algorithm for multilevel thresholding,” Journal of Information Science and Engineering, vol. 17, no. 5, pp. 713–727, 2001.