Porosity estimation method by X-ray computed tomography

Journal of Petroleum Science and Engineering 47 (2005) 209 – 217 www.elsevier.com/locate/petrol

Porosity estimation method by X-ray computed tomography H. Tauda,T, R. Martinez-Angelesa, J.F. Parrotb, L. Hernandez-Escobedoa a

Mexican Institute of Petroleum, P.O. Box 14-805, 07730 Mexico D.F., Mexico Institute of Geography, UNAM, P.O. Box 20-850, 01000 Mexico D.F., Mexico

b

Received 15 April 2004; received in revised form 7 March 2005; accepted 29 March 2005

Abstract In X-ray computed tomography imaging, the approaches used to determine the porosity of the rock from a single computed tomography scan are based on image segmentation techniques. When these techniques are applied to the same data, different results emerge and a threshold is needed at some level of the process. Consequently, the implication is that there is an uncertainty in the porosity measurement. Because of these sensibilities, a new method, called here the grey level method, is developed avoiding the use of these techniques. Considering the computed tomography image as a surface, the volumes required in porosity estimation are obtained by means of integrating this surface with simple operations applied to the image histogram. A porosity distribution which can reflect the properties of the studied rocks is developed given the value of the estimated porosity. The method is compared to two segmentation methods and is evaluated by a conventional one. A close agreement with the conventional method is found. D 2005 Elsevier B.V. All rights reserved. Keywords: Rock; Porosity; Image analysis; X-ray tomography

1. Introduction The productivity of wells in hydrocarbon exploration and prospecting depends on petrophysical properties such as porosity, permeability, saturation and capillarity. Porosity determines the storage capacity of hydrocarbons, permeability indicates the fluid flow capacity of rock, saturation reflects the porosity

T Corresponding author. Fax: +52 5 91 75 62 77. E-mail addresses: [email protected], [email protected] (H. Taud). 0920-4105/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.petrol.2005.03.009

occupied by hydrocarbons and capillarity reveals the available hydrocarbons that can be produced. Porous media can be considered as being the union of two parts: the solid part, made up of various materials and the empty part or pore space. Knowledge of the structures of the empty space provides a better measurement of these properties. The acquisition of reservoir material is expensive. Therefore, the non-destructive method seems to be the best way to obtain petrophysical properties. The application of X-ray computed tomography (CT) imaging to porous media has been used for many years to study and understand these properties

210

H. Taud et al. / Journal of Petroleum Science and Engineering 47 (2005) 209–217

(Anderson et al., 1988; Peters and Afzal, 1992; Johns et al., 1993; Kamath et al., 1998; Montemagno and Pyrak-Nolte, 1999; Ueta et al., 2000). Some kinds of rocks such as the hydrocarbon reservoir rock contain pores, which can have dimensions varying from microns to centimetres. The porosity is then characterized by the superposition of several sizes of pore (Moctezuma-Berthier et al., 2002). The porosity of rock is the ratio between the volume of empty space and the total volume of the rock and is reported in a percentage between 0 and 1. There are several ways to estimate porosity, for example: in the petrophysics laboratory with a porosimeter by injecting mercury or argon, or in the field by means of geophysical well logs or the buildup of a pressure test. Some investigations estimate porosity using the X-ray CT imaging. Withjack (1988) and Akin et al. (2000) determined the porosity by a dual scan at the same location obtained with different fluids saturating the porous medium. Akin et al. (1996) proposed a method employing a dual scan at two energy levels. These methods and others that use a dual scan with a single energy level are described and discussed by Akin and Kovscek (2003). Porosity can also be calculated from the CT image with a single scan of a rock by detecting the pore space by image segmentation techniques. Segmentation is the first treatment applied to CT images before analysing the physical characterization. It consists of the pores spaces extraction in a given scale corresponding to the CT image resolution. This step is crucial because of the nature of the CT image and the sensibility of image segmentation techniques (Ashbridge et al., 2003). It can reduce or increase the pores space and blur or make the connection appear between them. The properties of the feature resulting from the segmented image can vary greatly with small changes in the segmentation parameters. Consequently, the implication is that there is an uncertainty in the porosity measurement and its derived properties (Ashbridge et al., 2003; Sheppard et al., 2004). To bypass the segmentation step, the use of all the values of the CT image should be considered. In fact, without using the image segmentation of X-ray CT images, Keller (1998) measured the fracture aperture of several consolidated materials and characterized the aperture distribution statistically. He used the CT

numbers directly by integrating the fractured area in a typical cross-section and obtained an adequate predictor for fracture aperture for both small and large aperture. The use of integrating the CT numbers is the basic idea of our new method, called here the grey level method, which proposes to estimate the porosity of rocks from CT imaging by a single scan. This article is organized in the following manner. First, a brief description of the X-ray CT measurement is given. An overview of CT image segmentation in geological fields and their limitations is described and illustrated by applying two segmentation methods to three samples. The main lines of the method and its application to the CT image of these samples are explained. However, it is not the purpose of this article to provide either a detailed comparison of this method with the many other methods that exist or an in-depth evaluation of various methods. Nevertheless, the results are compared with those given by two segmentation methods and evaluated by those obtained from the conventional method.

2. X-ray CT measurement The principles of X-ray CT imaging have been discussed extensively elsewhere (Keller, 1998; Ohtani et al., 2000; Van Geet et al., 2000, 2001; Ketcham and Carlson, 2001; Mogensen et al., 2001; Akin and Kovscek, 2003). Principally, X-ray attenuation measurements taken around the object from different directions produce cross-sectional images of an object. Lambert–Beer’s law is used to relate the intensity of transmitted radiation I and the intensity of incident radiation I 0 to the material described by its attenuation coefficient l, and the distance travelled through the material x: I ¼ I0 elx

ð1Þ

The signal in each point in the reconstructed images, called the CT number, is expressed in Hounsfields units and is defined as: l  lw CT ¼ 1000 ð2Þ lw where l w is the attenuation coefficient of pure water. The calibration usually used gives CT = 0 for water and CT =  1000 for air. High CT numbers correspond

H. Taud et al. / Journal of Petroleum Science and Engineering 47 (2005) 209–217 Table 1 Core characteristics Field size (cm) Core diameter (cm) Pixel size (mm2) Slices

211

Frecuency core A

50000 Core A

Core B

Core C

core B

22.6 8.84 0.20  0.20 113

15.4 9.52 0.23  0.23 77

10.4 12.96 0.47  0.47 52

core C

to high density materials. In fact, the linear attenuation coefficient is a function of both the electron density (bulk density q) and the effective atomic number, Z in the following form (Van Geet et al., 2000; Akin and Kovscek, 2003):
 l ¼ q a þ bZ 3:8 =E 3:2 ð3Þ where a is the Klein–Nishina coefficient, b is a constant and E is the X-ray energy. As an absolute interpretation that the linear attenuation coefficient is quite difficult Van Geet et al. (2000) propose a dual energy method which allows quantification of density as well as effective atomic number. Akin and Kovscek (2003) point out that for X-ray energies above 100 kV, the CT image is proportional to density and for energies well below 100 kV the CT image is proportional to effective atomic number. Measurements with X-ray CT are subject to a range of errors and image artefacts including Beam hardening, star-shaped, positioning error and machine error. These artefacts and the techniques that reduced these artefacts are discussed by Van Geet et al. (2000), Ketcham and Carlson (2001) and Akin and Kovscek (2003). The CT number of each point or pixel is a function of the average density and composition of the material in a given volume or voxel. The voxel size is equal to the pixel size by the slice thickness. In reality, the twodimensional space (2D) of the CT image represents a

25000

0 -1000

CT number 1000

3000

5000

Fig. 2. Histograms of the three cores.

volume in three-dimensional space (3D). The resolution of the CT image produced depends on the scanner used. Fine details can be detected from the best resolution. CT image resolution or spatial resolution is estimated by measuring the two nearest distinguishable objects in two-dimensional CT image. The PICKER IQ PREMIER (IQXTRA) scanner parameters and calibration treatments were set to optimise the resolution and the quality of the scanned images. Tube voltage, tube current and exposure are 130.00 kV, 125.00 mA and 500.00 mAs respectively. The image matrix is 512  512 pixels with a slice thickness of 2.00 mm. Three cores, A, B and C, are used for this study. Their characteristics are shown in Table 1 where one slice of each core is illustrated in Fig. 1. Histograms in Fig. 2 describe the CT numbers and their frequency for each core. To obtain images with different characteristics, squared images Ae (300  300 pixels), B e (290  290) and C e (130  130) are extracted from circular images cores A, B and C respectively. Observing the CT images in Fig. 1, the lowest CT number is set to black and the highest to white. The increase in the CT number is seen as an increase in the grey level from dark to bright. The horizontal

Fig. 1. CT images with CT numbers increasing from black to white. From left to right, core A (from  991 to 5289), core B (from 939 to 3402), and core C (from 839 to 7000).

212

H. Taud et al. / Journal of Petroleum Science and Engineering 47 (2005) 209–217

4000 CT number

3000 2000 1000 0 -1000 0

190

350

Pixel

Fig. 3. CT numbers along a given section.

section in the plane of a given CT image points out that the CT numbers change gradually on the margins of cavities (Fig. 3). This gradual change is due to the absorption and dispersion of the X-ray, the properties of X-ray detectors and the process of image reconstruction. Moreover, the CT number in a voxel reflects not only the density and composition in a voxel, but also those in neighbouring voxels, because the spatial resolution is larger than the voxel size. Although CT numbers in air-filled cavities should be  1000 by definition, there are cavities corresponding to different CT numbers, going from the lowest numbers to the highest one. This indicates that the cavities cannot be estimated from CT numbers of  1000. Therefore, it is often difficult to select an optimal threshold by either visual interpretation or by automatic methods in order to detect the cavities from CT images.

3. CT image segmentation The common practice in geological fields employs thresholding techniques to segment such images by applying a visual interpreted threshold or image processing approach. Many automatic thresholding techniques are efficient when the grey level histogram is bi or multimodal but not in the case of uni-modal (Parker, 1996). More sophisticated methods have been developed in the medical field because of the simplicity of patterns studied and the thorough ¨ zkan et al., knowledge of their characteristics (O 1993; Chen et al., 1998; Pardo et al., 2001; Hibbard, 2004). In porous media, the geometry of pores space is more complex and irregular, therefore the segmentation is more complicated.

The methods used to segment the CT images assume prior knowledge of the type of material studied. Ohtani et al. (2000) used three-dimensional imaging of miarolitic cavities in a granite core and discussed the relationship between the cavity distribution and the geothermal reservoir. Based on the CT images study of synthetic glass tubes, they defined the threshold by comparing the CT image and the cut surface of the rock sample. Clark and Anderson (2002) compared soils under conventional and non tillage systems using X-ray CT. To segment the image into pores and solids, they chose a threshold by an analysis of the greyscale histograms and then by identification of the form in the image. Other studies make use of CT image for porosity measurements without indicating the image segmentation technique used (Mogensen et al., 2001). Even though, in other investigations, the image processing techniques are very developed because they are applied to binary images, segmentation is guided by an external supervisor. In the thresholding method developed by Oh and Lindquist (1999), the authors utilized the spatial covariance of the image in conjunction with Indicator Kriging to determine object edges. Implementation of this method requires a prior population identification by selecting two thresholds. This method was used to segment the CT images and then to investigate 3D geometry of porous media in binary images (Lindquist and Venkatarangan, 1999; Lindquist et al., 2000). Sarti and Tubaro (2002) detected and characterized planar fractures from 3D data relative to rock samples based on the Hough transformation. Their approach works on binary data obtained by a thresholding. 3D image segmentation and analysis techniques applied to binary images are developed to measure pore size distribution in soils (Delerue et al., 1999; Delerue, 2001; Delerue and Perrier, 2002). According to the nomenclature in image processing (Gonzalez and Woods, 1992), the 3D image presentation and description would describe their technique better than 3D segmentation. Sheppard et al. (2004) present a three-stage approach for the segmentation of images of porous and composite materials obtained from Xray tomography. After removing noise and enhancing the edges, they use a combination of watershed and active contour methods. To implement this method, a prior selection of two thresholds to define the studied

H. Taud et al. / Journal of Petroleum Science and Engineering 47 (2005) 209–217

regions and user parameters controlling the speed functions of the propagation of regional boundaries is needed. The segmentation methods, as shown above, need a threshold at some level of the process. Moreover, even when they are applied to the same data, these methods give different results (Parker, 1996), particularly when there is only one mode in the histogram of the data. The rocks studied here have this characteristic. The first porosity estimation of these rocks, in the absence of any laboratory measurement, is given by the segmentation method called here bUser methodQ. A threshold is located according to a visual interpretation of an expert on the subject. In fact, visualization provides a mechanism for identifying spatial distribution. However, the distinction between effects produced by a visualization process and those caused by an object being visualised is decisive in the identification. Displaying the same image with a variation of grey level can lead to an erroneous interpretation. To validate the accuracy of the user method, a second method in which no user gives the threshold is required. The clustering techniques offer this alternative. The aim of clustering is to partition a data set into subsets of bsimilarQ data. Closely linked to the pattern recognition field, clustering has also been applied to image segmentation. Many approaches to clustering have been suggested (Tsu and Gonzalez, 1974). K-means algorithm is chosen here because its application requires only one parameter; the number of classes to be found. For both methods, when a threshold is identified, the porosity is determined directly from the accumulated histogram. Table 2 records the threshold obtained by both methods. It shows: close thresholds for C by both methods, close thresholds for A and B when the K-means is applied and different thresholds for A and B when the user method is applied. A difference in threshold implies a difference in porosity estimation.

213

Table 3 Porosity (%) estimated by user and K-means methods User K-means

A

B

C

6.3 7.9

10.4 5.4

3.3 3.6

Comparing the results of the porosities, which is estimated by both methods in Table 3, shows a similar porosity in core C. However, the porosity in core A is greater than in core B when the K-means is applied. On the contrary the porosity in core A is lower than in core B in the case of the user method. In other words, the two methods give contrary results for cores A and B. Furthermore, a small difference in a threshold can reduce or increase the pores space and blur or make the connection appear between them. The methods developed for the texture analysis in a grey level image, such as the surface-curvature methods (Bergbauer and Pollard, 2003) based on a view of the grey level as a three-dimensional surface, improve on the pore skeleton detection in a given CT image resolution. Combining these methods with the CT numbers (research in progress) allows extensive studies of the void/matrix geometry. The application of the mean curvature in the CT image of core A underlines the concave and convex surfaces corresponding to the white and black tone respectively in Fig. 4A. Because the concave surfaces reflect the pore skeleton, Fig. 4A reveals the presence of pores varying from small to large. The intersection of this image with the segmented image by the user method shows that significant pores can be lost (Fig. 4B). Given these limitations, a method which does not need any threshold or any parameters at all must be considered.

Table 2 Threshold by user and K-means methods User K-means

A

B

C

3000 3129

3400 3092

3300 3327

Fig. 4. Mean curvature application. Grey tone between white or concave surface and black or convex surfaces. A: Features from core A image. B: Features combined with the segmented image.

214

H. Taud et al. / Journal of Petroleum Science and Engineering 47 (2005) 209–217

4. Method In the measure theory, the volume V X of a set of points X (in three-dimensional space) is defined as a measure of the set X. The measure of X can also be expressed as a function f, which associates a number V X with the set X: VX ¼ f ðXÞ

ð4Þ

Applying this concept to the case of rock, the porosity can be defined as the ratio between the volume of empty space V E and the total volume of the rock V T: /¼

VE f ðvoidsÞ ¼ f ðentire sampleÞ VT

rmax X

ðri  rmin ÞH ðri Þ

ð6Þ

ri ¼rmin

where s 2 is the pixel size on the horizontal space size, a is the pixel size in the vertical space, r i is a grey level related to the altitude and H(r i ) is the histogram of the image with grey levels in the range [r min, r max]. The histogram is a discrete function expressed as: H ðri Þ ¼

rmax X ni where H ðri Þ ¼ 1 n ri ¼rmin

ð7Þ

where n i is the number of pixels in the image with that grey level and n is the total number of pixels in the image. From the histogram, an accumulative histogram can be defined by: Ha ðrm Þ ¼

m X i¼0

H ðri Þ

V E ¼ s2 a

rmax X

ðrmax  ri ÞH ðri Þ

ð9Þ

ri ¼rmin

V T ¼ s2 a

rmax X

ðrmax  rmin ÞH ðri Þ

ð10Þ

ri ¼rmin

ð5Þ

To illustrate the estimation of the measure f, the CT image is symbolized as a Digital Terrain Model (DTM). The grey level in the DTM image is related to the altitude or elevation terrain. In a large DTM, the volume of terrain in general terms between two altitudes r min and r max can be approximated by: V ¼ s2 a

A plot of the functions H(r i ) or H a (r i ) for all values of r i provides a global description of the appearance of an image. The total volume V T of the rectangular parallelepiped which includes the DTM data and the volume of the empty space V E corresponding to the volume of complementary solid space can be expressed as:

ð8Þ

5. CT image application and results The application of the method to a CT image is based on the following considerations: (1) The energy of the scanner is greater than 100 kV, as mentioned above the CT image is proportional to density of the material. We then assume that there is a linear dependence between the CT numbers and the density. Keller (1998) indicates a good correlation between the CT number and fracture aperture with this assumption. Furthermore, in the segmentation method, the selection of a threshold implies the same assumption. (2) The 2D CT image space is set as 3D space with (x,y) denoting 2D position and the third coordinate (z) denoting the CT number. (3) All slice images form one image. (4) Each voxel must be considered to estimate porosity. (5) The CT number varies from a minimum negative value to a maximum positive one in the range [CTmin, CTmax]. To facilitate the volume computation, the CTmin =  1000 for air is subtracted from all the CT numbers to make them positive in the range [0, CTm] where CTm = CTmax+1000. (6) By similarity between the CT number and the grey level, r i is the CT number, n i is the number of pixels in the image with that CT number, n is the total number of pixels in the CT image, and r i = 0,1,. . .,r max and r max = CTm.

H. Taud et al. / Journal of Petroleum Science and Engineering 47 (2005) 209–217

Taking these considerations into account, the porosity can be calculated from Eqs. (9) and (10) as follows: rmax X

/¼

ðrmax  ri ÞH ðri Þ

rmax X

core A

0.4

core B core C

0.2

ð11Þ 0.1

rmax H ðri Þ

ri ¼0

l X

l X

ðl  ri ÞH ðri Þ

i¼0

l

CT number

0

The porosity / j in each slice image shows the distribution of the pore space in the core data and can be calculated from Eq. (11) by replacing H(r i ) with the histogram of slice image j H j (r i ). When there are a few isolated points with a very high CT number (r max = CTm is very large), the ratio between the volume of empty space V E and the total volume V T is disproportionate. Furthermore, r max needs to be large enough to enable a good assessment of the pore space. To find an adequate r max, the porosity distribution /(l) is calculated varying l from 0 to CTm and computing the corresponding volumes:

l X

¼

ðl  ri ÞH ðri Þ

i¼0

H ðri Þ

ð12Þ

lHa ðl Þ

i¼0

The plot of this distribution against the CT number, for the couple (A, Ae) data, shows that /(l) changes, increasing the CT number (Fig. 5). The /(l) curve decreases until a given position, corresponding to the minimum in the curve, and then, /(l) increases almost equally reflecting the effect of the few high isolated points. We assume that this position corresponds to the limit between the empty and solid spaces

0.5

Porosity

0.3

ri ¼0

/ðl Þ ¼

0.5

215

Porosity

A Ae

0.4

0

1000

2000

3000

4000

5000

6000

7000

Fig. 6. Porosity distribution for the three data A, B, and C.

indicating the adequate CT number for measuring the real porosity. The porosity estimation corresponds then to the minimum of the /(l) distribution. The difference between /(l) in the two data belonging to the same core for a lower CT number can be related to the core border effect that is simply the presence of more empty space in one data than in the other. In spite of this difference, the curves /(l) converge in the same porosity estimation. The same happens with the other couples. Moreover, the /(l) distribution indicates the nature of the studied cores. Cores A and B have the some geological characteristics in contrast to core C. In fact, the 3 cores come from the same lithology unit and are from the same age and origin (Cretaceous–Tertiary). The difference is assigned to the diagenesis process. Cores A and B are characterized principally by vugular structures and core C by fractures. Fig. 6 shows that cores A and B generally have the same distribution, and that core C is different. The distribution /(l) of core A is greater than the distribution of core B inside the CT number range [350, 1600]. The reverse happens in the range [1600, 3900]. In spite of this difference and the disparity between the values of CTm in cores A and B (isolated point at 4402 and 6289 respectively), their distribution converges at a similar limit point

0.3 0.2 0.1 CT number

0 0

1000

2000

3000

4000

5000

Fig. 5. Porosity distribution for data A and Ae.

6000

Table 4 Porosity (%) estimated by the grey level and conventional methods Grey level Conventional

A

B

C

7.9 9.4

7.3 8.9

3.1 5.1

216

H. Taud et al. / Journal of Petroleum Science and Engineering 47 (2005) 209–217

25

User

Porosity

method in that the porosity of cores A and B are similar (Table 4).

K-means

20

Proposed 15 10

6. Conclusion

5 0 1

21

41

61

81

101 Slice

Fig. 7. Porosity measurement in each slice of core A by the 3 methods.

(CT number c 3930) given a similar estimation of the porosity in both cores. Comparing the results of the porosities estimated by the proposed method (Table 4) and by the segmentation methods (Table 3), the porosity in core C is almost the same in all three methods. The result from the grey level method is comparable to the one obtained by K-means and coincides in the fact that the porosity in core B is not greater than the porosity in core A. Furthermore, this result shows that the porosity in core A and B is very similar. The porosity / j in each slice image versus the number of slices in core A is illustrated in Fig. 7. The curve obtained from the segmentation methods is more accentuated than the one obtained by the proposed method. This accentuation is related to the effect of classifying the CT numbers into just two groups: belonging or not belonging to the empty space. In contrast, this accentuation does not occur with the proposed method because it takes into account each voxel present in the image. Determination of the uncertainty related to an image-based measurement system depends on the propagation of the uncertainty for each step forming this system. It can be done by investigating models and methods for evaluating it (De Santo et al., 2000, 2004; Vieira and Paciornik, 2001; Ballester et al., 2002). This article does not propose to investigate these models and methods, nevertheless, in order to have a range of uncertainty and to have a reliable evaluation, the porosities estimated by the grey level method are compared and evaluated with those given by the conventional method (Helium injection at 400 PCI). A close agreement (about F2% porosity percent) is found. The results from the grey level method coincide with the conventional

A method has been developed to estimate the porosity of rock from single X-ray CT imaging. The porosity is calculated directly from the CT numbers without using segmentation techniques and without using any user-defined parameters. From a novel distribution /(l) which can indicate the properties of the studied rocks, the minimum of the curve gives the value of the estimated porosity. The method is applied to 3 samples chosen to represent difference in porosity. In fact, 2 segmentation methods by threshold provide different porosities but give a contrary result for 2 samples. The distribution /(l) corresponding to these 2 samples shows a similar value of their porosities, and the results agree with those obtained by the conventional method. The method evaluation is made by comparing the results with those obtained by the conventional method. A close agreement with the conventional method is found. In terms of operations applied to the image histogram, the algorithm gives a simple method for estimating the porosity.

Acknowledgements We thank Anna Reid for improving the English grammar and style of our manuscript, the reviewers for their helpful criticism and suggestions and the Mexican Petroleum Institute (IMP) for its support.

References Akin, S., Kovscek, A.R., 2003. Computed tomography in petroleum engineering research. In: Mees, F., Swennen, R., Van Geet, M., Jacobs, P. (Eds.), Application of X-ray Computed Tomography in the Geosciences. Special Publication - Geological Society of London, vol. 215, pp. 23 – 38. Akin, S., Demiral, B., Okandan, E., 1996. A novel method of porosity measurement utilizing computerized tomography. In Situ 20, 347 – 365. Akin, S., Schembre, J.M., Bhat, S.K., Kovscek, A.R., 2000. Spontaneous imbibition characteristics of diatomite. J. Pet. Sci. Eng. 25 (3 and 4), 149 – 165.

H. Taud et al. / Journal of Petroleum Science and Engineering 47 (2005) 209–217 Anderson, S.H., Gantzer, C.J., Boone, J.M., 1988. Rapid non destructive bulk density and soil water content determination by computer tomography. Soil Sci. Soc. Am. J. 52, 35 – 40. Ashbridge, D.A., Thorne, M.S., Rivers, M.L., Muccino, J.C., O’Day, P.A., 2003. Image optimization and analysis of synchrotron X-ray computed microtomography (CAT) data. Comput. Geosci. 29, 823 – 836. Ballester, M.A.G., Zisserman, A.P., Brady, M., 2002. Estimation of the partial volume effect in MRI. Med. Image Anal. 6, 389 – 405. Bergbauer, S., Pollard, D.D., 2003. How to calculate normal curvatures of sampled geological surfaces. J. Struct. Geol. 25, 277 – 289. Chen, Ch.W., Luo, J., Parker, K.J., 1998. Image segmentation via adaptive-mean clustering and knowledge-based morphological operations with biomedical applications. IEEE Trans. Image Process. 7 (12), 1673 – 1683. Clark, J.G., Anderson, S.H., 2002. Computed tomographic measurement of macroporosity in chisel-disk and no-tillage seedbeds. Soil Tillage Res. 64 (1 and 2), 101 – 111. Delerue, J.F., 2001. Segmentation 3D, Application a` l’extraction de Re´seaux de Pores et a` la Caracte´risation Hydrodynamique des Sols. The`se de doctorat, Universite´ Paris XI, no d’ordre 6574, 194 pages. Delerue, J.F., Perrier, E., 2002. DXSoil, a library for 3D image analysis in soil science. Comput. Geosci. 28 (9), 1041 – 1050. Delerue, J.-F., Perrier, E., Yu, Z.Y., Velde, B., 1999. New algorithms in 3D image analysis and their application to the measurement of a spatialized pore size distribution in soils. Phys. Chem. Earth, Part A Solid Earth Geod. 24 (7), 639 – 644. De Santo, M., Liguori, C., Pietrosanto, A., 2000. Uncertainty characterization in image-based measurements: a preliminary discussion. IEEE Trans. Instrum. Meas. 49 (5), 1101 – 1107. De Santo, M., Liguori, C., Paolillo, A., Pietrosanto, A., 2004. Standard uncertainty evaluation in image-based measurements. Measurement 36, 347 – 358. Gonzalez, R.C., Woods, R.E., 1992. Digital Image Processing. Addison-Wesley, New York. Hibbard, L.S., 2004. Region segmentation using information divergence measures. Med. Image Anal. 8, 233 – 244. Johns, R.A., Steude, J.D., Castanier, L.M., Roberts, P.V., 1993. Non destructive measurements of fracture aperture in crystalline rock cores using X-ray computed tomography. J. Geophys. Res. 98, 1889 – 1900. Kamath, J., Xu, B., Lee, S.H., Yortsos, Y.C., 1998. Use of pore network models to interpret laboratory experiments on vugular rocks. J. Pet. Sci. Eng. 20, 109 – 115. Keller, A., 1998. High resolution, non-destructive measurement and characterization of fracture apertures. Int. J. Rock Mech. Min. Sci. 35 (8), 1037 – 1050. Ketcham, R.A., Carlson, W.D., 2001. Acquisition, optimization and interpretation of X-ray computed tomographic imagery: applications to geosciences. Comput. Geosci. 27 (4), 381 – 400. Lindquist, W.B., Venkatarangan, A., 1999. Investigating 3D geometry of porous media from high resolution images. Phys. Chem. Earth, Part A Solid Earth Geod. 24 (7), 593 – 599.

217

Lindquist, W.B., Venkatarangan, A., Dunsmuir, J., Wong, T.F., 2000. Pore and throat size distributions measured from synchrotron X-ray tomographic images of Fontainebleau Sandstones. J. Geophys. Res., B [Solid Earth Planets] 105 (9), 21509 – 21527. Moctezuma-Berthier, A., Vizika, O., Adler, P.M., 2002. Macroscopic conductivity of vugular porous media. Transp. Porous Media 49, 313 – 332. Mogensen, K., Stenby, E.H., Zhou, D., 2001. Studies of water flooding in low-permeable chalk by use of X-ray CT scanning. J. Pet. Sci. Eng. 32 (1), 1 – 10. Montemagno, C.D., Pyrak-Nolte, L.J., 1999. Fracture network versus single fractures: measurement of fracture geometry with X-ray tomography. Phys. Chem. Earth, Part A Solid Earth Geod. 24 (7), 575 – 579. Oh, W., Lindquist, W.B., 1999. Image thresholding by Indicator Kriging. IEEE Trans. Pattern Anal. Mach. Intell. 21, 590 – 602. Ohtani, T., Nakashima, Y., Muraoka, H., 2000. Three-dimensional miarolitic cavity distribution in the Kakkonda granite from borehole WD-1a using X-ray computerized tomography. Eng. Geol. 56 (1 and 2), 1 – 9. ¨ zkan, M., Dawant, B.M., Maciunas, R.J., 1993. Neural-networkO based segmentation of multi-modal medical images: a comparative and prospective study. IEEE Trans. Med. Imag. 12 (3), 534 – 544. Pardo, X.M., Carreira, M.J., Mosquera, A., Cabello, D., 2001. A snake for CT image segmentation integrating region and edge information. Image Vis. Comput. 19, 461 – 475. Parker, J.R., 1996. Algorithms for Image Processing and Computer Visions. Wiley Computer Publishing, New York. Peters, E.J., Afzal, N., 1992. Characterization of heterogeneities in permeable media with computed tomography imaging. J. Pet. Sci. Eng. 7, 283 – 296. Sarti, A., Tubaro, S., 2002. Detection and characterisation of planar fractures using a 3D Hough transform. Signal Process. 82 (9), 1269 – 1282. Sheppard, A.P., Sok, R.M., Averdunk, H., 2004. Techniques for image enhancement and segmentation of tomographic images of porous materials. Physica A 339, 145 – 151. Tsu, J.T., Gonzalez, R.C., 1974. Pattern Recognition Principles. Addison-Wesley, New York. Ueta, K., Tani, K., Kato, T., 2000. Computerized X-ray tomography analysis of three-dimensional fault geometries in basementinduced wrench faulting. Eng. Geol. 56 (1 and 2), 197 – 210. Van Geet, M., Swennen, R., Wevers, M., 2000. Quantitative analysis of reservoir rocks by microfocus X-ray computerised tomography. Sediment. Geol. 132 (1 and 2), 25 – 36. Van Geet, M., Swennen, R., Wevers, M., 2001. Towards 3-D petrography: application of microfocus computer tomography in geological science. Comput. Geosci. 27 (9), 1091 – 1099. Vieira, P.R.M., Paciornik, S., 2001. Uncertainty evaluation of metallographic measurements by image analysis and thermodynamic modeling. Mater. Charact. 47, 219 – 226. Withjack, E.M., 1988. Computed tomography for rock-property determination and fluid-flow visualization. Soc. Pet. Eng. Form. Eval. 3, 696 – 704.