Research Article  Open Access
An Improved Algorithm in Porosity Characteristics Analysis for Rock and Soil Aggregate
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 Xray computed tomography (CT) under uniaxial compressive loading. Xray 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 Xray. 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 Xray 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].
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 costeffective 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, highvolume 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 Xray tubes becoming a standard tool in geophysical laboratories, there are many researches on the threedimensional distribution of grains and voids in natural rock over the last 10 years [8–13]. Xu et al. [8] compared the nature soilrock and simulatedrain soilrock mixtures. Vallejo and Mawby [9] analyzed the stability of slopes made of granular materialclay mixtures. Hirono et al. [10] used Xray to visualize the advecting fluid image during permeability testing under atmospheric pressure. In addition to natural rock, Xray computed tomography (CT) as a tool is also applied to analyze the property of micromechanical investigations of cement and mortar [14], geomaterialssoils [15], and the rock mesodamage propagation law [16]. Threedimensional 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 “rocklike” material based on Xray 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 beamhardening 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 BeamHardening
A highresolution Xray 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 Xray is 50 mm Fe, the spatial resolution is 4 Lp/mm, and the minimum focusobjectdistance is 0.13 mm. Figure 2 shows the 450 KV Xray computed tomography machine [7].
Image a traditional chest Xray, obtained by placing the subject against an Xray sensitive plate and “illuminating” the individual with an Xray beam in the form of a cone. The Xray plate produces an image whose intensity at a point is proportional to the Xray 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 Xray beam.
Since all substances attenuate lowenergy Xrays more strongly than highenergy ones, primarily because of photoelectric absorption, a heterogeneous beam traversing an absorbing medium becomes proportionately richer in highenergy 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 Xray is monochromatic, which means the energy obeys Beer’s law, while in fact Xray 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, beamhardening 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 beamhardening is that the gray center is smaller than the edge.
(a)
(b)
It is well known that beamhardening 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 beamhardening. 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 betweenclass variance of the thresholded image [15]:
For bilevel thresholding, Ostu verified that the optimal threshold is chosen so that the betweenclass 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 beamhardening, 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 singlethreshold 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)
The improved betweenclass 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 betweenclass 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 beamhardening 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.
(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 strainstresses. The CT sliceporositystress strain with different stressstrain is shown in Figure 8.
From Figure 8, we can see that when the strainstress 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 strainstress reaches 3.87 MPa, RSA will be in the rocksoil jointed cracks with initiation stage and the stable crack growth stage, and porosity increases quickly. When the strainstress 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.

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)12 stage is corresponding to crack closure and soil consolidation stage, and the porosity decreases from 9.55% to 6.55%;(2)23 stage is corresponding to linear elastic deformation stage, and the porosity is 6.32%. After the closure of the microcrack, the stressstrain curve exhibits a linear elastic relationship;(3)34 stage is corresponding to rocksoil 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.
5. Conclusion
With the help of Xray 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 beamhardening, 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 realtime. 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 strainstress, 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 beamhardening 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 BioMedical Electronics, vol. 8, pp. 68–72, 1961. View at: Publisher Site  Google Scholar
 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 graylevel picture thresholding using the entropy of the histogram.,” Computer Vision, Graphics, and Image Processing, vol. 29, no. 3, pp. 273–285, 1985. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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, “Mesodamage cracking characteristics analysis for rock and soil aggregate with CT test,” Science China Technological Sciences, vol. 57, pp. 1–11, 2014. View at: Google Scholar
 W. J. Xu, R. L. Hu, and R. J. Tan, “Some geomechanical properties of soilrock mixtures in the Hutiao Gorge area, China,” Geotechnique, vol. 57, no. 3, pp. 255–264, 2007. View at: Publisher Site  Google Scholar
 L. E. Vallejo and R. Mawby, “Porosity influence on the shear strength of granular materialclay mixtures,” Engineering Geology, vol. 58, no. 2, pp. 125–136, 2000. View at: Publisher Site  Google Scholar
 T. Hirono, M. Takahashi, and S. Nakashima, “In situ visualization of fluid flow image within deformed rock by Xray CT,” Engineering Geology, vol. 70, no. 12, pp. 37–46, 2003. View at: Publisher Site  Google Scholar
 R. A. Ketcham, “Threedimensional grain fabric measurements using highresolution Xray computed tomography,” Journal of Structural Geology, vol. 27, no. 7, pp. 1217–1228, 2005. View at: Publisher Site  Google Scholar
 T. Ohtani, T. Nakano, Y. Nakashima, and H. Muraoka, “Threedimensional shape analysis of miarolitic cavities and enclaves in the Kakkonda granite by Xray computed tomography,” Journal of Structural Geology, vol. 23, no. 11, pp. 1741–1751, 2001. View at: Publisher Site  Google Scholar
 A. Vervoort, M. Wevers, R. Swennen et al., “Recent advantages of Xray CT and its applications for rock material,” in Xray CT for Geomaterials, J. Otani and Y. Obara, Eds., pp. 79–91, A.A. Balkema, Kumamoto, Japan, 2003. View at: Google Scholar
 E. N. Landis, “Xray tomography as a tool for micromechanical investigations of cement and mortar,” in Advances in XRay Tomography for Geomaterials, J. Desrues, G. Viggiani, and P. Be'suelle, Eds., pp. 79–93, Antony Rowe Ltd, Chippenham, UK, 2006. View at: Google Scholar
 J. Otani and Y. Obara, “Xray CT for geomaterialssoils, concrete, rocks,” in Proceedings of the International Workshop on XRay CT for Geomaterials, pp. 227–287, A.A. Balkema Publishers, Kumamoto, Japan, November 2003. View at: Google Scholar
 X. Ge, J. Ren, Y. Pu, W. Ma, and Y. Zhu, “Realin time CT test of the rock mesodamage propagation law,” Science in China, Series B: Chemistry, vol. 44, no. 3, pp. 328–336, 2001. View at: Publisher Site  Google Scholar
 E. N. Landis and E. N. Nagy, “Threedimensional work of fracture for mortar in compression,” Engineering Fracture Mechanics, vol. 65, no. 23, pp. 223–234, 2000. View at: Publisher Site  Google Scholar
 I. S. Darma, T. Sugiyama, and M. A. B. Promentilla, “Application of Xray CT to study diffusivity in cracked concrete through the observation of tracer transport,” Journal of Advanced Concrete Technology, vol. 11, pp. 223–234, 2013. View at: Google Scholar
 N. Otsu, “A threshold selection method from graylevel histogram,” IEEE Transactions on System Man Cybernetics, vol. 9, no. 1, pp. 62–66, 1979. View at: Publisher Site  Google Scholar
 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. View at: Google Scholar
Copyright
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.