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 |

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

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

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

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.

Tomographic Images - slice views

7.5% porosity | 13% porosity | 15% porosity | 22% porosity |

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.

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.

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

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.

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.

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.

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.

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.

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

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.

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.