3DMA-Rock

A Software Package for
Automated Analysis of Rock Pore Structure
in 3-D Computed Microtomography Images

W. B. Lindquist


Software Authors:   Prof. W.B. Lindquist,   Dr. S.-M. Lee,  Dr. W. Oh,  Dr. A.B. Venkatarangan,  Dr. H. Shin,  Dr. M. Prodanovic

Documentation Authors:   M. Prodanovic,   W.B. Lindquist

Table of Contents

A)   Introduction

B)   Pore Space Analysis

     B.I)    Segmentation
     B.II)   Construction and Modification of the Medial Axis
     B.III)  Throat Construction
     B.IV)  Pore Surface Construction
     B.V)   Pore-Throat Network
     B.VI)  Geometric Characterization

C)   Fluid Analysis

     C.I)   Segmentation
     C.II)   Fluid Blobs
     C.III)  Pore Saturations

D) References




A) Introduction   ( return to table of contents)   ( next section)

3DMA-Rock is a software package for analyzing the pore space in three- (and two-) dimensional X-ray computed microtomographic (XCMT) images of rock.  With appropriate images it can also be used to analyse fluid partitioning in the pore space.

This web page provides an overview of the capabilities of the package, which are continually evolving. Documentation has been made available for five released versions (12/99, 12/03, 03/04, 08/05 and 12/05). For detailed documentation developed for researchers interested in using the code, consult the links in the table below. In particular, detailed exposition of the capabilities of each version are found under "Detailed Instructions for Running 3DMA-Rock". In general the documentation for a later version over-rides the documentation for a previous version. The one exception is the documentation for the 12/99 version which contains information on the formats of data files produced by the 3DMA-Rock code that is not found in the documentation for later versions.

3DMA-Rock
System Requirements and Graphics Software
Source Code Availability
Distribution Installation and Compiling Instructions
Detailed Instructions for Running 3DMA-Rock
    Version 12/99 (3DMA General Users Manual)
    Version 12/03
    Version 03/04
    Versions 08/05 and 12/05
Release Notes and Bug Reports

The figures accompanying this html file are best viewed with the html viewer window set to full screen.




B) Pore Space Analysis   ( prev section)   ( return to table of contents)   ( next section)

The image analysis algorithms for the pore space consist of six general steps:  image segmentation, extraction and modification of the medial axis of the pore space, throat construction using the medial axis as a search path, pore surface construction via marching cubes, assembly of the pore-throat network, and geometrical characterization of the pore-throat network.  These steps are discussed below, largely by illustration, with accompanying literature references. 




B.I) Segmentation   ( prev section)   ( return to table of contents)   ( next subsection)


Fontainebleau Sandstone
Tomographic Images - slice views
7.5% porosity 13% porosity 15% porosity 22% porosity
IK Segmented Images - slice views

Fig. 1: (top) Sample regions of slices from 3-D tomographic images of a suite of Fontainebleau cores imaged at 5.7 micron resolution. (bottom) The same slice-regions after the data has been segmented using the indicator-kriging [IK] ualgorithm.

The greyscale intensity of each voxel in an XCMT image (Fig. 1 (top)) is typically an integer value (from 0-255 for 8-bit images). 3DMA-Rock provides three options for image segmentation:  simple segmentation, Mardia-Hainsworth (MH) [1988] and indicator kriging (IK) [Oh and Lindquist, 1999].  For simple segmentation, the user must enter a single threshold:  voxels having intensity above this threshold are assigned to the grain phase; remaining voxels are assigned to the background phase.  Both MH and IK methods are locally adaptive and are similar in that both utilize (an estimate of) the two-point correlation function of the image. MH is an iterative method, requiring a single threshold value which it modifies each iteration. IK is a "single pass" method requiring that a subpopulation of voxels of each phase be positively identified a priori.  This can be done by manually establishing a window of intensity values delimited by two thresholds. Intensity values above the upper threshold are assumed grain, those below the lower threshold are assumed void space. The classification of remaining voxels is estimated by indicator kriging which utilizes an estimate of a correlation function incorporating local spatial information.  These three methods are compared in Oh and Lindquist on a test images (Figs. 2, 3) and on a high resolution (laser scanning confocal microscopy) image (Fig. 4) of Berea sandstone.

test image simple Mardia-Hainsworth indicator kriging

Fig. 2: Comparison of simple thresholding, MH and IK segmentation methods on a test image of overlapping circles to which uncorrelated log-normal Gaussian random noise has been applied.

test image Mardia-Hainsworth indicator kriging

Fig. 3: Comparison of simple thresholding, MH and IK segmentation methods on a test image of overlapping circles to which spatially correlated log-normal Gaussian random noise has been applied.

test image Mardia-Hainsworth indicator kriging

Fig. 4: Comparison of MH and IK segmentation methods on a high resolution (1 micron voxel size) laser scanning confocal microscopy image of Berea sandstone.




B.II) Construction and Modification of the Medial Axis   ( prev subsection)   ( return to table of contents)   ( next subsection)

The medial axis [Sirjani and Cross, 1991] of a digitized object is a 26-connected (or 6-connected) centrally-located skeleton which preserves the original topology and geometry of the object. If their are no cavities in the object, the medial axis is a network of voxel paths. Erosion-based algorithms are used to construct the medial axis; 3DMA-Rock employs the algorithm developed by Lee et al. [1994].  The medial axis is sensitive to surface features, including noise; as a result, the medial axis may contain spurious paths that are not significant descriptors of the overall geometry of the object.  Spurious paths can be trimmed; 3DMA-Rock provides a menu of options to control this trimming.  Fig. 5 shows the medial axis extracted from a small region of a 22\% porosity Berea data set (imaged at 4.93 micron voxel resolution) before and after spurious paths are trimmed. After trimming, all "apparent" dead end branches in fact exit through the boundary of the cubic region displayed.

before trimming after trimming

Fig. 5:  (left) The medial axis extracted from a sandstone image.  (right) The medial axis, after trimming all dead end (non-flow supporting) paths.  Grains are not displayed in the image as this would obscure most of the view of the medial axis. The colors represent, in rainbow scale, distance to the closest (at least two) grain voxels.

The importance of the medial axis trimming is more clearly pictured in 2D. In Fig. 6 we present the medial axis computed from a single slice of an image of vesiculated basalt, before and after trimming of dead end paths.

before trimming after trimming

Fig. 6:  (left) The 2-D medial axis extracted for a single slice of an image of vesiculated basalt.  (right) The medial axis, after trimming all dead end (non-flow supporting) paths.  Grains are displayed in black.

For more information on the medial axis algorithms, see the Online Manual to 3DMA.




B.III) Throat Construction   ( prev subsection)   ( return to table of contents)   ( next subsection)

The non-symmetrical shape of pores results in some medial axis paths lying completely with a pore, while other paths connect pores. We use a distance measure [Lindquist, 2002b], comparing the length of a path with the distance to the nearest grain surface as a first pass to eliminate pore-internal paths. For the rest, we have developed two throat finding algorithms [Venkatarangan, 2000; Shin, 2002] to determine the location of one (or more as a path may pass through several pores) minimal area cross sectional surfaces - i.e. throats. Our algorithms do not presuppose that the throat surface are planar. The throat surfaces are determined as triangulated interfaces; the set of voxels through which each triangulated throat surface passes defines a throat region.

Fig. 7: Plot of one of the triangulated throat surfaces (purple) showing the triangulated pore-grain surface (green) in the immediate vicinity




B.IV) Pore Surfaces   ( prev subsection)   ( return to table of contents)   ( next subsection)

With throat regions identified, we utilize the marching cubes algorithm [Bloomenthal, 1988; Lorensen and Cline, 1987] to determine the surface of each pore. The surface is the triangulated interface separating the pore from the set of grain/throat region voxels surrounding it. (See Fig. 8).

Fig. 8: Rotated views of the surface surrounding a single pore. The pore interior is internal to the surface displayed. View of the interior is via the throats of the pore.




B.V) Pore-Throat Network   ( prev subsection)   ( return to table of contents)   ( next subsection)

With throats located, the pore space can be divided into a network of pores separated by throat surfaces. Every pore is cross-indexed with its connecting throats and ajoining pores, and every throat is cross-indexed with the pores it connects. Note that this subdivision does not include an attempt to defince channels connecting pores** - rather two ajoining pores connect directly at the mutual throat surface. Thus the volume of the pore space is divided into pore bodies. The code allows the option to ignore all pores that touch the boundary of the imaged region, as such pores are truncated, with subsequent loss of knowledge of their true properties. Thus throats may be indicated as having only one connecting pore, and some pores may have no indicated ajoining pore through an existing throat, indicating that a truncated pore has been dropped from the cross-indexing.


**We know of no way geometrically to define where a pore ends and a connecting channel begins.




B.VI) Geometric Characterization   ( prev subsection)   ( return to table of contents)   ( next section)

As rock properties are statistical, we characterize the spatial variation of each property by plotting its distribution of values over the image. Distributions are produced for: pore coordination number, pore volume, pore surface area, and throat surface area. Using moment of inertial computations we compute a center of mass and principal directions for each pore. A diameter, passing through the center of mass in each principal direction is also computed. An effective radius for a pore can also be computed using the sphere of equivalent volume. in a similar manner, two principal diameters and an effective radius (from the circle of equivalent area) are produced for each throat surface. Distributions of the principal diameters and the effective radius are produced for the pores and throats. Cross correlations between each pair of variables may also be of interest. Currently we produce cross correlation plots for: coordination number and pore volume; throat and pore effective radii.

Fig. 9:  A summary of statistics obtained from an investigation of a suite of Fontainebleau sandstones, from 7.5% to 22\% porosity. Voxel resolution 5.687 microns.
(top left) The exponential distribution for the pore coordination number shown for the 18% porosity sample is characteristic of the entire range of porosity investigated. The inset plot shows values obtained for the exponential decay parameter C_O from the analysis of 16 images (4 at each porosity).
(top right) The channel length distribution (defined as pore center to pore center distance along the medial axis path joining neighboring pores) shown for the 18% porosity sample has an exponentially decaying tail characterisitic of the entire range of porosity investigated. The inset plot shows values obtained for the exponential decay parameter L_O from the analysis of 16 images (4 at each porosity).
(bottom left) The log-normal-like distribution for the throat surface area shown for the 18% porosity sample is characteristic of the entire range of porosity investigated. The inset plot shows values obtained for the mean and standard deviation of the throat surface area over the range of porosity investigated. Mean throat surface area increase from 7.5% to 18% porosity.
(bottom right) The log-normal distribution for the pore volume shown for the 18% porosity sample is characteristic of the entire range of porosity investigated. The inset plot shows values obtained for the mean and standard deviation of the pore volume over the range of porosity investigated. Mean pore volume decreases slightly with porosity!

Fig. 10:  A summary of statistics obtained from an investigation of a Berea core, voxel resolution 4.93 microns.
(top left) Pore volume.
(top middle) Distributions for pore principal axes diameters (D1,D2 and D3), effective radius (Reff), and diameters in an arbitrary direction (Dx). The diameter plots have been vertically offset from one another for clarity. The effective radius is plotted with the same vertical offset as D3 to show the similarity of distribution.
(top right) Pore surface area.
(bottom left) Throat surface area.
(bottom right) Distributions for throat principal axes diameters (D1 and D2), and effective radius. The diameter plots have been vertically offset from one another for clarity. The effective radius is plotted with the same vertical offset as D2 to show the similarity of distribution.




C) Fluid Analysis   ( prev section)   ( return to table of contents)   ( next section)

Using X-ray absorbing dopants (iodine and bromine) to increase X-ray contrast between fluids [Seright et al. 2002; 2003; 2004; Prodanovic et al. 2004] and careful tomographic registration to re-image the same region under a sequence of fluid injections, 3DMA-Rock, version 08/05, has been extended to examine fluid distribution in the pore space.

Fig. 11:  Example from the analysis of brine and oil distribution in the pore space of Berea at residual fluid conditions. A water based gel is injected at one point, to investigate brine and oil distributions before and after gel placement.




C.I) Fluid Segmentation   ( prev subsection)   ( return to table of contents)   ( next subsection)

With a single fluid occupying the pore space, the segmentation described in B.I will differentiate the pore space and grain matrix. If the same region is subsequently imaged so that individual voxels align, when there are two fluids (having sufficient X-ray attenuation contrast) in the pore space the grain space can be identified from the single fluid picture, and the segmentation utilized in B.I can now be used to identify the fluid phase of each pore voxel.

Fig. 12:  (a) A subregion of a single slice of a tomographic image stack of Berea sandstone with brine in the pore space. (b) The same subregion with brine and oil (hexadecane) in the pore space at residual oil conditions. (c) Pore/rock and brine/oil segmentation results for this region.




C.II) Fluid Blobs   ( prev subsection)   ( return to table of contents)   ( next subsection)

With fluid phases identified, a simple grass-fire technique is used to locate all disconnected blobs of each phase. (Consistency in a digitized image requires on phase be treated as 26-connected, the other as 6-connected. We treat the wetting phase as 6-connected.) The size distributions for the blobs of each phase are illustrated in Fig. 13.

Fig. 13:  Example from the analysis of brine and oil distribution in Berea sandstone from the sequence of floods illustrated in Fig. 11 on an initially water saturated rock consisting of (from top to bottom): oil injection to residual water (S_wr); water injection to residual oil (S_or); water-based gel injection, shut-in and curing at high temperature to complete polymerization (gel); oil injection to residual water (gelS_wr); and water injection to residual oil (gelS_or). With 4.93 micron voxel resolution, in S_wr the residual water phase appears as a continuous spectrum of disconnected water blobs up to size 0.02 mm^3. The oil phase consists of a single, large (80% pore occupancy) oil blob and a spectrum of smaller disconnected blobs below 0.002 mm^3 volume. Analagous statements can be made for the fluid partitioning at the other residual conditions. On the middle (gel) plot, the pore volume distribution has been superimposed to demonstrate that the continuum spectrum of small blob sizes is always bounded by the size of the largest pore, with one exception. In the gel flooding stage (S_or to gel), the injected gel mobilizes upstream residual oil, into the imaged region. This modilized oil appears to merge into blobs of size exceeding the larges pores. Full discussion of the results can be found in Prodanivic et al. [2004].




C.III) Pore Saturations   ( prev subsection)   ( return to table of contents)   ( next subsection)

Full analysis (B.I - B.VI) of the image of the single phase saturated sample produces the entire pore-throat network. With throat locations known, fluid contents can be determined on a pore-by-pore bases in the subsequent images of two-phase occupancy at residual conditions. Fig. 14 displays results for the flooding study discussed in Figs. 11 and 13.

Fig. 14:  Distribution of oil saturations (at residual water conditions) and water saturations (at residual oil conditions) as a function of pore volume at the residual fluid conditions investigated. For each volume range, we present a box plot summarizing the distribution of saturations measured in all pores having volume in the specified range.




D) References   ( prev section)   ( return to table of contents)

J. Bloomenthal.  Polygonization of implicit surfaces.  IEEE Comput. Graph. App. 5 (1988) 341.

T. C. Lee, R. L. Kashyap, and C. N. Chu.  Building skeleton models via {3-D} medial surface/axis thinning algorithms.  CVGIP:  Graph. Models Image Process. 56 (1994) 462-478.

W.B. Lindquist, S.M. Lee, D. Coker, K. Jones and P. Spanne.  Medial axis analysis of void structure in three-dimensional tomographic images of porous media.  J. Geophys. Res. 101B (1996) 8297.  ftp://ftp.ams.sunysb.edu/pub/papers/1995/susb95_01.ps.gz

W.B. Lindquist and A. Venkatarangan.  Investigating 3D geometry of porous media from high resolution images.  Phys. Chem. Earth (A) 25 (1999) 593.  ftp://ftp.ams.sunysb.edu/pub/papers/1998/susb98_01.ps.gz

W.B. Lindquist, A. Venkatarangan, J. Dunsmuir and T.-f. Wong.  Pore and throat size distributions measured from sychrotron X-ray tomographic images of Fontainebleau sandstones.  J. Geophys. Res. 105B (2000) 21508.  ftp://ftp.ams.sunysb.edu/pub/papers/1999/susb99_13.pdf

W.B. Lindquist.  Network flow model studies and 3D pore structure.  Contemp. Math. 295 (2002a) 355.  ftp://ftp.ams.sunysb.edu/pub/papers/2001/susb01_14.pdf

W.B. Lindquist.  Quantitative analysis of three dimensional X-ray tomographic images.  In Developments in X-ray Tomography III, U. Bonse (ed.), Proceedings of SPIE 4503 103.  SPIE, Bellingham, WA, 2002b.  ftp://ftp.ams.sunysb.edu/pub/papers/2001/susb01_06.pdf

W.E. Lorensen and H.E. Cline.  Marching cubes: a high resolution 3-D surface construction.  ACM Comput. Graph. 21 (1987) 163.

K.V. Mardia and T.J. Hainsworth.  A spatial thresholding method for image segmentation.  IEEE Trans. Pattern Anal. Machine Intell. 10 (1988) 919-927.

W. Oh and W. B. Lindquist.  Image thresholding by indicator kriging.  IEEE Trans. Pattern Anal. Mach. Intell. 21 (1999) 590-602.  ftp://ftp.ams.sunysb.edu/pub/papers/1998/susb98_02.ps.gz

M. Prodanovic, W.B. Lindquist, and R.S. Seright.  3D microtomographic study of fluid displacement in rock cores.  To be presented at the Computational Methods in Water Resources XV Conference, Chapel Hill, NC, June 13--17.  ftp://ftp.ams.sunysb.edu/pub/papers/2004/susb04_01.pdf

R.S. Seright, J.-T. Liang, W.B. Lindquist, and J.H. Dunsmuir.  Characterizing disproportionate permeability reduction using synchrotron X-ray computed microtomography.  SPE Reserv. Eval. Eng. 5 (2002) 355-364.  ftp://ftp.ams.sunysb.edu/pub/papers/2002/susb02_02.pdf

R.S. Seright, J. Liang, W.B. Lindquist and J.H. Dunsmuir.  Use of X-ray computed microtomography to understand why gels reduce permeabilityto water more than to oil.  J. Petroleum Sci. Eng. 39 (2003) 217-230.  ftp://ftp.ams.sunysb.edu/pub/papers/2002/susb02_03.pdf

R.S. Seright, M. Prodanovic, and W.B. Lindquist.  X-ray computed microtomography studies of disproportionate permeability reduction.  SPE paper #89393, to presented at the 14th SPE/DOE Symposium on Improved Oil Recovery, Tulsa, OK, Apr. 17-21, 2004.  ftp://ftp.ams.sunysb.edu/pub/papers/2003/susb03_19.doc

H. K. Shin.  A Throat Construction Algorithm for Medial Axis Analysis of 3D Images of Vesiculated Basalts.  PhD Thesis, Department of Applied Mathematics and Statistics, Stony Brook University. 2002.

A. Sirjani and G. R. Cross.  On representation of a shape's skeleton.  Pattern Recogn. Lett. 12 (1991) 149-154.

S.R. Song, K.W. Jones, W.B. Lindquist, B.A. Dowd and D.L. Sahagian.  Synchrotron X-ray computed microtomography: studies on vesiculated basaltic rocks.  Bull. Volcanol. 63 2001 253-263.

P. Spanne, J.F. Thovert, C.J. Jacquin, W.B. Lindquist, K.W. Jones and P.M. Adler.  Synchrotron Computed Microtomography of Porous Media: Topology and Transports.  Phys. Rev. Lett. 73 (1994) 2,001-2,004.

A.B. Venkatarangan.    PhD Thesis, Department of Applied Mathematics and Statistics, Stony Brook University, 2000.