An Improved Algorithm in Porosity Characteristics Analysis for Rock and Soil

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. 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 themain reason of RSA cracking is the elasticitymismatch 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.


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 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 threedimensional distribution of grains and voids in natural rock over the last 10 years [8][9][10][11][12][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].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 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.

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 focusobject-distance is 0.13 mm. Figure 2 shows the 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 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  0 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 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 512 × 512 pixels with gray levels from 1 to 255.The main phenomenon for the beamhardening is that the gray center is smaller than the edge.
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 beamhardening.In the following section, we will introduce our improved Ostu's method for image thresholding.

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,  1 with gray levels [1, 2, . . ., ] and  2 with gray levels [ + 1, . . ., ].Then, the gray level probability distributions for the two classes are where Also, the means for classes  1 and  2 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 between-class variance  2  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  2  , are functions of  1 and  2 .The two optimum threshold values,  * 1 and  * 2 , are the values that maximize  2  ( 1 ,  2 ).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  1 .Next,  2 is incremented through all its values greater than  1 but less than  − 1 (i.e.,  2 =  1 + 1, . . .,  − 2).Then,  1 is incremented to its next value and  2 is incremented again through all its values greater than  1 ; this procedure is repeated until  1 =  − 3. The result of this process is a 2D array,  2  ( 1 ,  2 ), and the last step is to look for the maximum value in this array.The values of  1 and  2 corresponding to that maximum are the optimum thresholds,  * 1 and  * 2 .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  * 1 in (15), and the red circle in the right image represents the larger thresholding, as the  * 2 in (15).
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 (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  , 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.

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

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

Figure 3 :
Figure 3: Relationship between attenuation coefficient and thickness of the absorber.

Figure 6 :
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).

Figure 7 :
Figure 7: Histogram and smoothing histogram of CT slice.

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

Figure 9 :
Figure 9: Relationship between RSA porosity and axial stress.

Table 1 :
The porosity calculation for typical CT slice under different stress level.