3DMA-Rock Primer
August 2005, December 2005

This page assumes you have obtained a copy of either of the 08/05 or 12/05 versions of the 3DMA-Rock code, compiled it and need to start using it. If you want to obtain a copy of the code, email Prof. Brent Lindquist at lindquis@ams.sunysb.edu.

Readers are reminded that the 3DMA General Users Manual, which accompanied the 12/99 version of 3DMA-Rock, contains information on the formats of data files produced by the 3DMA-Rock code.

This "primer" is a tutorial on the basic algorithms ('cases' in our jargon) to be run to produce an analysis of the 'pore' space and/or the distribution of resident fluid phases in a porous medium.

Table of Contents

Introduction and motivation
                    (return to table of contents)    (next section)

Given a 3D tomographic image of a porous medium, and potentially a sequence of images of the same medium with pores filled with more than one fluid, the (2D) sketch below outlines a typical 3DMA-Rock based analysis:


Pore structure characterization

Pore/grain
------------------>
segmentation

Medial axis
construction
------------->

+ some modification

Tomographic image.
(polyethylene core is light gray and pore space filled with oil is dark gray).
  Segmented Image.
Grain (polyethylene) is shown in black, and pore space in white).
  Medial axis.
Pixel color ranges from minimal burn number) to violet (maximal burn number).
                   |
           |    Throat finding +
           |    pore partitioning
          v
   
Pore/throat characterization
(distributions or correlations of properties such as pore volume, surface area, shape factor, coordination number, throat area)
<------------------
        Pore-throat network.
The medial axis is shown in red (note it is trimmed on the boundary); identified throats are shown in lilac.


Fluid characterization

Fluid
------------------>
segmentation
------------------>
Fluid characterization
(interphase surface area, blob characterization, scale fluid saturation, ...)
Tomographic image.
Polyethyelene, water and oil are light, darker and the darkest gray.
  Fluid segmented image.
Water is green and oil is red.
   



Directory set-up, download example input files
                     (previous   section)    (return to table of contents)    (next section)

- Create a directory for your data analysis. We will use the example directory /nfs/3dma_rock/Test/

  cd /nfs/3dma_rock
  mkdir Test

- Organize the following subdirectories under /nfs/3dma_rock/Test/:

  cd Test
  mkdir sw
  cd sw
  create_wdir


create_wdir is a shell script (3dma_rock/bin/create_wdir) provided with the code that creates the following list of subdirectories in the present working directory:

  tomog cases burn seg c_seg ma ma_t plots results throats

Except for checking to make sure write permissions exist and the subdirectories do not already exist, running create_wdir is equivalent to running the system command

  mkdir tomog cases burn seg c_seg ma ma_t plots results throats

Here is a summary of what will be stored in each subdirectory:

tomog   raw, tomographic files
cases   input, output files
burn   erosion number (burn) files
seg   segmented files
c_seg   cleaned segmented files
ma   medial axis files
ma_t   trimmed medial axis files
plots   graphical output, largely for quality control
results   graphical output on geometrical characterization
throats   throat and pore network files

- Use the following link to download an archived (sw-tomog.tar.gz) data set of example tomographic input files into your /nfs/3dma_rock/Test/sw/tomog/ directory.

Warning: This archived data set is 14.6 MBytes in size. The ftp transfer may take some time.

The data set contains 256 files labeled 'sw.XXX.gz' where XXX ranges from 001 to 256, containing the slices (the z-direction) of a 256x256x256 volume of a Berea sandstone core saturated with water. Each slice is an array of 256^2 numbers (unsigned chars, values ranging from 0 to 255) representing X-ray attenuation coefficients.

- Unzip and extract the tomographic files.

   gunzip sw-tomog.tar.gz
   tar -xvf sw-tomog.tar



Running the code
                      (previous section)    (return to table of contents)    (next section)

3DMA-Rock is a command-line based code (no GUI yet !!), querying the user with a scrolled menu of options. To avoid the agony of typing in answers, each functional step of 3DMA-Rock can read these answers from a pre-prepared input file with answers provided in the expected order. To facilitate the preparation of input files, 3DMA-Rock echos every menu and its user response into an output file, which can used later as an input file.

We have adopted a convention for naming the input/output file pairs associated with each functional step. The input file is labelled caseN.M.in and its associated output file is caseN.M.out where N and M are the integer answers provided to the first two menus presented by 3DMA-Rock.

The first menu is

	Input Data Options
		tomographic data (1)
		segmented data (2)
		burn data (3)
		medial axis data (4)
		throat data (5)
		pore-throat network data (6)
		fluid data (7)
	Enter choice:

and refers to the (intermediate) set of data files on which the desired analysis is to be performed.

The second menu, which also takes a numerical answer, depends on the first menu selection.

As stated, the 3DMA-Rock executable is built to read user responses either from the keyboard or from an input file (caseN.M.in), and produce an output file (caseN.M.out) containing a record of the input as well as a text record of pertinent output statistics. The output file caseN.M.out also contains information on the run time and total memory required to perform the computation. This provides critical "scale-up" information for performing runs on larger data sets. In addition, depending on the case, 3DMA-Rock outputs intermediate data files, as well as files containing requested graphs and images.

To run 3DMA-Rock with user input from the keyboard:

	3dma_linux > caseM.N.out

Inputting answers from the keyboard can be tedious. In addition, error recovery from incorrect input is poor, input must be entered all over again. Once prepared however, an output file can be converted to an input file, edited to change responses, and rerun making subsequent running relatively easy.

To reuse an output file:

	mv caseM.N.out caseM.N.in
	vi caseM.N.in                 /* edit input file, change answers */
	stripcomm_linux -b < caseM.N.in | 3dma_linux > caseM.N.out

The stripcomm facility is provided for the purpose of running from an input file. It scans the input file searching for the next occurrance of the string ": ". The remainder of the line is taken to be the next required answer and is passed to the 3DMA-Rock executable.

As the stripcomm command line is rather lengthy, we provide the shell script runcase to shorten typing. Thus

	runcase linux caseM.N
is equivalent to typing the above stripcomm command.



List of input files
                      (previous section)    (return to table of contents)    (next section)

The following is the list of the main data analysis algorithmic operations ("cases") available under this version of 3DMA-Rock. For each case, an example input file is provided for the user to run on the tomographic data set downloaded as instructed in the Directory set-up, download example input files section above. As is clear from the algorithmic descriptions, most of these cases must be run in the correct sequential order.

Input file   Function
     
    Operations on Tomographic Data
case1.1.in   slicewise raster plot of tomographic data
case1.1a.in   slicewise raster plot of tomographic data with plot merging option
case1.6.in   histogram of attenuation coefficients, dependency of porosity on simple thresholding cut-off
case1.5smpl.in   segmentation by simple thresholding
case1.5krig.in   segmentation by indicator kriging
     
    Operations on Segmented Data
case2.1ras.in   slicewise raster plot of segmented data
case2.1gmv_sub.in   Geomview format subvolume plots of the pore-grain surface as generated via the Marching-Cubes algorithm
case2.1gmv_vox.in   Geomview format subvolume plots of the voxelized pore-grain surface
case2.3.in   volume information on disconnected pores and grains
case2.4.in   "clean-up" operations on the segmented data
case2.5.in
case2.5seal.in
  erosion computation of the media axis (Lee, Kasyap, Chu [1994] algorithm)
case2.17.in   generation of a segmented file of a digitized hexagonal close packing of spheres
     
    Operations on and using Medial Axis Data
case4.10.in
case4.10seal.in
  trimming operations on the medial axis
case4.1ras.in   slice-by-slice raster plots of the medial axis
case4.1gmv.in   Geomview format volume plots of the medial axis
case4.1gmv_sub.in   Geomview format subvolume plots of the medial axis
case4.6x.in   Medial axis path tortuosity computation (in the x-direction)
case4.9.in   conversion from cluster-path to voxel-based medial axis file formats
case4.12dij.in
case4.12dij_seal.in
  computation of pore throats using the Dijkstra based algorithm [Venkatarangan 2000]
case4.12wdg.in
case4.12wdg_seal.in
  computation of pore throats using the wedge based algorithm [Shin 2002]
case4.12sub_seal.in   computation of pore throats using the Dijkstra algorithm with subsectioning [Prodanovic 2005]
case4.17.in   conversion of medial axis data from padded (sealed) indicing to unpadded (unsealed) indicing
     
    Operations on Throat Data
case5.1ras.in   raster plot format slice-by-slice views of throats
case5.1gmv_sub.in   Geomview format subvolume plots of throat regions
case5.2gmv_sub.in   Geomview format subvolume plots of throat surfaces
case5.3dij.in
case5.3.in
  construction of the pore-throat network
case5.5wdg.in
case5.5mrg.in
  throat file information; recording of MA paths that do not have an identified throat
case5.6mrg.in
case5.6.in
  comparison and merging of throat data files
case5.7.in   Conversion of throat data files from padded (sealed) indicing to unpadded (unsealed) indicing
     
    Operations on Pore Network Data
case6.1dij.in   Geomview format plots of individual pore surfaces; pore surface area computation
case6.1sor.in   Geomview format plots of individual pore surfaces and fluid partitioning (within the pore)
case6.2dij.in   output of pore-throat statistics
     
    Operations on Fluid Data
case7.1.in   histogram of tomographic data in the pore space
case7.2.in   fluid segmentation (indicator kriging segmentation in the pore space)
case7.3.in   fluid saturation in pores/throats
case7.5.in
case7.5a.in
  fluid blob size (volume) distribution
case7.8ras.in   slicewise raster plots of fluid segmented files
case7.8sub.in   Geomview format subvolume plot of the water-oil interface
case7.8area.in   total fluid/fluid interface area computation
case7.9.in   slicewise comparative raster plots of segmented, tomographic and fluid segmented data
     
    Single Phase lattice Boltzmann Computations
case8.1thr10.in   lattice Boltzmann computation and visualization of the permeability for a single throat
case8.1thr.in   lattice Boltzmann computation of the permeabilities for all throats




Case 1.1:   Plotting the input (tomographic) image
                      (previous section)    (return to table of contents)    (next section)

Case 1.1 plots images of the input tomographic data.

Test example:
	To obtain slice-by-slice raster files (2-D plots of each slice in the
	stack in SUN raster file format which is viewable almost everywhere):
	
	Download the input file case1.1.in into the
	/nfs/3dma_rock/Test/sw/cases/ directory.

	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case1.1
	
	This will produce the output file case1.1.out and raster files
	sw.100.ras.gz to sw.110.ras.gz in the directory /nfs/3dma_rock/Test/tomog
	These files can be viewed by most 2D image viewing packages, e.g. "xv"
	or "gimp" (the latter comes with a RedHat Linux distribution).
	
	% cd ../tomog/
	% xv sw.*.ras.gz

	Fig. 1 shows the plot in sw.110.ras.gz.


Fig. 1:   Rasterfile image of the test example sw.110.ras.gz. The light phase represents water, the dark is the sandstone matrix.

The example input file case1.1a.in is a variation of case1.1 in which the "Merge raster files" option is activated, resulting in the merging of each subsequent set of 10 slice images on the same page in a 2 x 5 row/column arrangement.




Case 1.6:   Image intensity Histogram and porosity-threshold dependence
                      (previous section)    (return to table of contents)    (next section)

The first major algorithm in the 3DMA-Rock analysis is segmentation. Producing a distribution of the X-ray attenuation coefficients in the image is an decision making tool towards segmentation. Case 1.6 will produce a histogram plot (jgraph format) of the distribution density as well as a cumulative probability plot. The cumulative plot is interpretable as the fractional volume occupied by all voxels whose intensity lies below any X-ray attenuation coefficient value, and thus directly provides a measure of porosity versus cut-off threshold for the global (simple) segmentation method.

Test example:
	Download the input file case1.6.in 
	into the /nfs/3dma_rock/Test/sw/cases/ directory.

	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case1.6

	This will produce the output file case1.6.out and the jgraph file
	poros_thres.jgr in the directory /nfs/3dma_rock/Test/sw/results.
	The jgraph0 file can be processed and viewed as follows:

	% cd ../results/
	% jgraph < poros_thres.jgr > poros_thres.ps
	% ggv poros_thres.ps

	Fig. 2 shows the resulting plot in poros_thres.jgr.

0 The jgraph program is publically available software that is distributed with the 3DMA-Rock code. It produces postscript output which can be printed, included in documents, or viewed by a postscript interpreter (e.g. ghostscript).


Fig. 2:   X-ray attenuation coefficient (scaled 0-255) density (dashed line, right axis - not normalized to unity) and the cumulative probability (solid line, left axis) interpreted as void fraction (porosity).

As is typical with X-ray CT images of bi-phase materials having sufficient X-ray attenuation contrast, the plot is bimodal, the lower mode corresponding to the lower attenuating phase (water filled pore space in this case) and the upper mode to the higher attenuating phase (sandstone).




Case 1.5(smpl):   Simple Segmentation
                      (previous section)    (return to table of contents)    (next section)

Simple (global) thresholding interprets each voxel having attenuation coefficient value below a threshold as phase "0" and every voxel above as phase "1".

Choosing a reasonable threshold value of 80 from the histogram of case1.6, (Fig. 2) results in the following test example run.

Test example:
	Download the input file case1.5smpl.in  into the
 	  /nfs/3dma_rock/Test/sw/cases/ directory.
	
	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case1.5smpl
	
	This will produce the output file case1.5smpl.out,
	intermediate (segmented) data file sw_smpl.gz in the 
	/nfs/3dma_rock/Test/sw/seg directory, and the raster files 
	sw_smpl.XXX.ras.gz in the ../seg/ directory.
	
	The raster files provide a view of the segmentation results
	for each slice and just echo1 the output of the
	sw_smpl.gz (volume) data file.
	
	% cd ../seg/
	% xv sw_smpl.*.ras.gz

	The plot for one slice is shown in Fig. 3.

1So why print both the data files sw_smpl.gz and the raster files sw_smpl.XXX.ras.gz which, after all, contain the same information? The reason is one of disk usage. The data files, which we refer to as segmented files, contain 0/1 data and are therefore bit-packed. The raster files are byte-packed (character packed) and are 8 times larger. The raster files can be deleted (and reproduced any later time by running case2.1) and the segmented data files stored with minimized impact on disk space.

Fig. 3:   View of the test example simple segmentation results for slice 110 in the 3D image stack. Water is white; sandstone matrix is black.

Simple segmentation is extremely prone to "speckled noise" and is only employed as a fast, rough, segmentation method.




Case 1.5(krig):   Segmentation by indicator kriging
                      (previous section)    (return to table of contents)    (next section)

The indicator kriging method requires input of two thresholds, T0 and T1. T0 (T1) establishes a lower (upper) limit below (above) which any voxel can reasonably be identified as phase 0 (1). All other voxels (with coefficient values in the range [T0,T1]) have their phase determined by indicator kriging.

Test example:
	Download the input file case1.5krig.in into the
 	  /nfs/3dma_rock/Test/sw/cases/ directory.

	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case1.5krig

	This will produce the output file case1.5krig.out
	(in the cases directory), the intermediate data file
	sw.gz in the /nfs/3dma_rock/Test/sw/seg directory, and
	kriging summary raster files sw.XXX.ras.gz in the
	sw/seg directory.

	The raster files provide a visual synopsis (Fig. 4) of
	the kriging decisions for each slice.
	
	% cd ../seg/
	% xv sw.*.ras.gz

Fig. 4:   Kriging summary raster file image for slice 110 in the test example, with the attenuation coefficient histogram of the tomographic data superimposed lower right.  (upper left) The tomographic data.  (upper right) White (black) voxels are identified as phase 0 (1) through the threshold T0 (T1). Yellow (red) voxels are identified as phase 0 (1) via kriging.  (lower left) The final segmented result.  (lower right) Superimposed attenuation coefficient histogram with the indicator kriging thresholds indicated as light vertical lines.

With the image segmented, all further processing steps utilize the segmented data volume file (see Appendix A for more details on the format of this file).




Cases 2.3/2.4:   Post segmentation clean up procedures
                      (previous section)    (return to table of contents)    (next section)

No segmentation procedure is perfect - some level of segmentation errors result. In some respects, the most troubling of these are small isolated phase "blobs". Isolated blobs of phase 12 can be especially troubling if medial axis analysis of phase 02 is to be performed, as the isolated phase 1 blobs are cavities in phase 0 and will cause the medial axis of phase 0 to consist of a network of linear paths and closed surfaces. The closed surfaces surround the phase 1 blobs. A mixed path/surface medial axis network is difficult to analyse. The essential feature of Case 2.4 is the ability to "delete" isolated blobs of a given phase (by converting them to the opposite phase). We use a volume threshold to decide what blobs to convert. In order to determine volume thresholds for isolated blob conversion of each phase, Case 2.3 should be run first. Case 2.3 provides both text printout and a jgraph histograms3 of the distribution of isolated blob sizes. Unless there is some interest in retaining blobs above a certain size, in practice we convert all isolated rock and all isolated pore blobs. While isolated pore spaces may be real in some rocks, they are certainly cut off from fluid flow; thus for fluid flow questions, their existence is unimportant. Isolated rock "blobs" are unphysical, except in one instance. The boundary of the imaged region may cut through a projecting rock grain, leaving the grain "tip" as a visible protrusion into the imaged volume, but losing the connectivity information of the tip. In such cases, we consider it a small error to delete the protruding tip (by converting it to pore phase). Certainly the expense of keeping it, resulting in the introduction of a surface segment in the medial axis, vastly complicates using the medial axis as a search network.


2In the output for Case 2.3, the word "pore" refers to phase 0 voxels (for the example data provided, this is in fact the water phase), while the word "grain" refers to phase 1 (for the example data provided, this is the sandstone grain phase).
3In practice the jgraph histogram is of little use; necessary information can be obtained directly from the text printout.
Test example:
	Download the input files case2.3.in and case2.4.in
	into the /nfs/3dma_rock/Test/sw/cases/ directory.

	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case2.3

	This produces the output file case2.3.out and a jgraph file
	results/disc_vol.jgr.  The distribution of volume sizes (in
	numbers of voxels) ordered by volume for both disconnected pore and
	grain volumes is printed in the case2.3.out file - as shown here.
	Note that cumulative number, cumulative volume and cumulative volume
	fraction (relative to the entire phase volume) of blobs up to the
	certain size are stated as
	well.
	
	To delete ALL but the largest blob of each phase, volume thresholds
	(e.g. 32000 for phase 0 (water) and 1000 for phase 1 (sandstone)) can
	be established simply by examining this printout. For instance,
	cumulative volume fraction of all disconnected grain volumes of up to
	1000 voxels in size is 0.000194 or less than 0.02% of the entire grain
	phase volume.  The chosen thresholds can then be entered into
	case2.4.in (In fact we chose the threshold value of 1000 for both in
	the example case2.4.in provided.)
	
	% vi case2.4.in
	% runcase linux case2.4
	
	This stores a "cleaned" segmented volume file, sw.gz, and raster
	file plots, sw.XXX.ras.gz, of the slices in the c_seg/ directory.
	
	% cd ../c_seg/
	% xv sw.*.ras.gz

A second feature in "clean-up" of the segmented image is the ability to "smooth" the void/grain interfaces. This has two identifiable goals, reducing the number of dead-end paths on the medial axis that my result from small interface irregularities, and the elimination of "1-voxel wide pore channels" running between grains that may result from segmentation/digitization artifacts. (The latter of these two goals is more critical than the former.)
In particular, if the segmented data is to be reconstructed as a pore/throat network, in case2.4 we strongly endorse either
   applying the morphological closure operation to the grain space
or
   performing the material/void boundary smoothing operation option:
     2) less than a majority of neighbors of the same type
during the clean up procedure in case2.4.

There are other options in case2.4.in that are somewhat self-explanatory. To assist in their understanding, (and for a further discussion of the two features already mentioned) a case 2.4 input file with extra commentary can be accessed here.

Analysis continues with the processed segmented volume file in the /nfs/3dma_rock/Test/sw/c_seg/ directory.



Case 2.1:   Visualization of segmented files and interphase surface area measurement
                      (previous section)    (return to table of contents)    (next section)

2D slicewise raster plots of segmented files have been encountered in the previous sections. A handy input file that will provide raster plots of cleaned-up segmented files can be downloaded from case2.1ras.in.

3D visualization routines have been included to provide plots of the pore-grain interface (surface). The pore-grain surface can be plotted as either: a surface composed of (rectangular) voxel faces (producing a plot which is a 3D analog of the 2D raster plots) in either Geomview or Inventor format; or as a triangulated surface obtained via the Marching Cubes algorithm [Lorensen and Cline, 1987]. This latter surface plot is available only in Geomview format, and outputs surface area measurement as well.

Test example:
	Download the input files case2.1gmv_sub.in, case2.1gmv_vox.in
	into the /nfs/3dma_rock/Test/sw/cases/ directory.
	
	% cd /nfs/3dma_rock/Test/sw/cases/ 
	% runcase linux case2.1gmv_sub
	% runcase linux case2.1gmv_vox
	
	This produces the output files case2.1gmv_sub.out and case2.1gmv_vox.out 
	as well as Geomview files, seg_surf_sub.list and seg_surf_vox_data.geom.gz
	(zipped file), in the plots/ directory.
 
	% cd ../plots
	% geomview seg_surf_sub.list
	% gunzip seg_surf_vox_data.geom.gz
	% geomview seg_surf_vox_data.geom

	Figs. 5 and 6 show the Geomview plots. Note that file case2.1gmv_sub.out
	also reports the measured surface area in physical units (typically
	microns2).
pore-grain surface
Fig. 5:  Geomview graphics view of the pore-grain surface, obtained by the Marching Cubes algorithm, of an [0,99]x[0,99]x[0,99] voxel subvolume of the sample.
  
contact voxel surfaces
Fig. 6:  Geomview graphics view of the pore-grain voxel surfaces of the same subvolume.



Case 2.5:   Medial axis construction
                      (previous section)    (return to table of contents)    (next section)

The processed segmented files are used to compute the medial axis of the pore space using the LKC [Lee, Kashyap, Chu; 1994] algorithm. In addition an integer distance value (referred to as a "burn" or "erosion" number) is associated with each pore voxel. This distance is the max norm distance from the pore voxel to the closest grain voxel in units of voxels. (A pore voxel sharing a face with a grain voxel has "burn" number 1. The layer of grain voxels bordering the pore space have burn number 0.)

Test example:
	Download the input file case2.5.in  into the
	/nfs/3dma_rock/Test/sw/cases/ directory.
	
	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case2.5 &
	
	This produces the output file case2.5.out,
	intermediate burn number data files sw.XXX.gz 
	in the burn/ directory and intermediate medial 
	axis data files sw.XXX.gz in the ma/ directory.

The burn number files record distance information for all pore voxels (L_infinity distance to the closest grain). All grain voxels are assigned a burn number of 0. The medial axis data files record locations of the medial axis voxels only. Complete information on the structure of these files is found in the appendices of the 3DMA General Users Manual.



Case 4.10:   Medial axis trimming
                      (previous section)    (return to table of contents)    (next section)

Any medial axis algorithm is sensitive to surface noise, resulting in dead-end paths that have no relevance to the geometrical features under investigation. These should be removed before the axis is utilized either to characterize the object or as a search path for feature location. Isolated pore spaces will produce isolated "pieces" of the medial axis which may or may not be of interest. Using terminology from graph-theory to refer to components of the medial axis, we describe case4.10 which provides the user with the ability to trim features from the medial axis. The major components of this trimming are:

Output from case4.10 is a trimmed medial axis, which can either be stored in "voxel-location" format files, or preferably in "cluster-path" format which retains all book-keeping information on (merged) clusters and remaining paths. Complete information on the structure of the cluster/path files is found in the appendices of the 3DMA General Users Manual.

Test example:
	Download the input file case4.10.in into the
	/nfs/3dma_rock/Test/sw/cases/ directory.
	
	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case4.10
	
	This produces the output file case4.10.out and
	trimmed medial axis data files cp_sngl, cp_loc and 
	cp_struct of cluster/path format in the ma_t/ directory.

Analysis continues using the trimmed medial axis files stored in the directory /nfs/3dma_rock/Test/sw/ma_t/



Case 4.1:   Medial axis plotting
                      (previous section)    (return to table of contents)    (next section)

The medial axis can be viewed either slicewise (raster file format), or in 3D (Geomview or Inventor format). Each format has advantages and disadvantages. The continuity of the medial axis is difficult to follow in slice-by-slice viewing, but the relationship of the medial axis to the grain phase is clearer. 3D viewing is preferable for seeing the continuity of the medial axis, however the simultaneous plotting of the pore-grain surface hides much of the medial axis and makes it somewhat difficult to view the medial-axis/grain surface relationship. Plotting the medial axis without the grains has the downside that it makes perspective viewing of the medial axis difficult.

We recommend plotting the medial axis in color, where the color spectrum provides information on the burn number distribution along the medial axis. Typically we use a "rainbow" color scale, small burn numbers being red, large burn number being blue/violet.

Test example:
	Download the input file case4.1ras.in  and case4.1gmv.in 
	into the /nfs/3dma_rock/Test/sw/cases/ directory.

	% cd /nfs/3dma_rock/Test/cases/
	% runcase linux case4.1ras
	% runcase linux case4.1gmv

	Case 4.1ras produces raster plots sw.XXX.ras.gz in the ma_t/
	directory. An example slice plot is shown in Fig. 7. Case 4.1gmv
	produces a Geomview format plot ma.list in the plots/ directory.
	A view of this 3D plot is shown in Fig. 8.
 
Fig. 7:  A raster plot of the (3D) medial axis voxels visible in slice 110 of the example data.   Fig. 8:  A Geomview graphics view of the 3D medial axis of the sample data set. For clarity of viewing (and to greatly reduce file size), grains are not plotted.



Case 4.6:   Geometrical tortuosity for shortest paths through the medial axis
                      (previous section)    (return to table of contents)    (next section)

Given a (trimmed) medial axis computed for the pore space inside a rectangular volume and output in MA voxel list format, the tortuosity of paths through the medial axis network can be computed. (If the medial axis has been output in cluster/path format, it can be converted using Case 4.9, "convert MA list to/from node/path formats".)

Here is a quick description of the algorithm used to compute geometrical tortuosities.

tortuosity directions
Fig. 9:  Labeling of 3D image faces in Case 4.6.

Consider the computation of tortuosities for paths connecting "entrance" face 3 to "exit" face 4. (See Fig. 9.)

For each MA voxel i on face 3 and for each MA voxel j on face 4, the algorithm will find the shortest path (if one exists) from i to j through the medial axis network. Let l_ij be the length of this shortest path, and d_ij be the Euclidean distance between voxels i and j. The geometrical tortuosity t_ij of the path is defined as t_ij = l_ij / d_ij. Geometrical tortuosities are computed for the shortest paths between all possible entrance/exit voxel pairs i, j on the MA network.

Tortuosity results are summarized in two graphs output to a single jgraph file. The first is a histogram of the geometrical tortuosities. The second histograms the distribution of the values t_i = { min_j (t_ij) | j is an exit face MA voxel}, where i ranges over all MA voxels on the extrance face.

Test example:
	Download the input files case4.9.in and case4.6x.in
	into the /nfs/3dma_rock/Test/sw/cases/ directory.

	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case4.9

	This will convert MA cluster/path files to MA voxel location format 
	needed by case4.6. The voxel location files will be "../ma_t/sw.XXX.gz"
 
	% runcase linux case4.6x

	Two geomview format plots, tort_x_abs and tort_x_all, are placed in
	the plots/ directory as well as the tortuosity summary plot which is
	put in the results/ directory. The plot tort_x_abs provides a 3D
	view of the shortest MA paths running from face (3) to (4).

	% cd ../plots
	% geomview tort_x_abs 

	A view of this plot is shown in Fig. 10. The plot file tort_x_all
	contains all of the MA paths running from (3) to (4). This plot is
	superimposed on Fig. 10.

        % cd ../results 
	% jgraph < tort_x.jgr > tort_x.ps
	% ggv tort_x.ps

	A plot of the tortuosity summary results in tort_x.jgr is presented
	in Fig. 11.
paths plot
Fig. 10:  Geomview plot of medial axis paths from face (3) to face (4). The shortest paths for each node on face (3) are in red.
  
tort stats
Fig. 11:  Histogram of geometrical tortuosities of all (upper graph) and shortest (lower graph) medial axis paths from voxels on face (3) to voxels on face (4).



Case 4.12:   Throat computaton
                      (previous section)    (return to table of contents)    (next section)

3DMA-Rock provides three throat computation algorithms. One, developed as the thesis of Dr. A. Venkatarangan, utilizes Dijkstra shortest path algorithms to locate the bounding perimeter of the throat surface. The second, developed as the thesis of Dr. H. Shin is refered to as the 'wedge finding' algorithm as it employs wedges in a cylindrical coordinate system to construct throats.  The third, a subsectioning throat algorithm is a variation of the Dijkstra based algorithm and was developed as a part of the thesis of Dr. M. Prodanovic to address higher porosity media.

Note: The Dijkstra shortest path algorithm is very computationally expensive. Runs of 30 hours or more to compute throats on a Berea image sample containing 450x450x450 voxels are common. The wedge finding algorithm provides substantial speed up in throat finding as it uses approximations to throat surfaces while searching for the minimal one on each path. The subsectioning Dijkstra based algorithm has run times that are, on average, three (!) times that of the Dijkstra algorithm and should be used only within the framework described in the Aggressive Throat Computation section.

Test example:
	Download the input files case4.12dij.in and case4.12wdg.in
	into the /nfs/3dma_rock/Test/sw/cases/ directory.

	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case4.12dij

	This produces a file sw_dij.gz4 in the ..throats/ directory
	containing information on throats found by the Dijkstra algorithm.
	On this sample, Case 4.12(dij) will run approximately
	25 minutes on a 2.0GB RAM, 1.7GHz Intel-Pentium processor.
	
	% runcase linux case4.12wdg
	
	This produces a file sw_wdg.gz4 in the ../throats/ directory
	containing information on throats found by the wedge algorithm.
	On this sample, Case 4.12(wdg) will run approximately
	7 minutes on a 2.0GB RAM, 1.7GHz Intel-Pentium processor.

4A throat file contains complete information on each throat, including it's location on the medial axis, the list of tis perimeter voxels, and its surface area. Complete information on the structure of the throat files is found in the appendices of the 3DMA General Users Manual.


Aggressive Throat Computation - medium/high porosity media
                      (previous section)    (return to table of contents)    (next section)

In moderate to high porosity (40-50%) samples, successful throat finding and pore partitioning (case 5.3) can be especially challenging. While the Berea sandstone presented in this primer does not fall into this category, we use it to illustrate a "work-around" for more complete throat finding in higher porosity media.

In analyzing 40-50% porosity polyethylene samples (Prodanovic, Lindquist and Seright) we found that the largest pore accounted for more than 80% of the entire pore space volume. Two effects were responsible for this:


5 These algorithms are bench-marked on "nice" data sets, such as sphere packs. Local geometry in real media can often exceed the parameter ranges or implicit design assumptions under which the algorithms succeed.
Test example:
	Create a directory needed to store the sealed image information

	% cd /nfs/3dma_rock/Test/
	% mkdir sw_seal
	% cd sw_seal
	% create_wdir

	In this case the sw_seal/seg and sw_seal/tomog subdirectories created
	by create_wdir will not be used. Rather the starting data will come
	from the unsealed, cleaned segmentation files in /nfs/3dma_rock/Test/sw/c_seg.

	Download the input file case2.14.in into the
	/nfs/3dma_rock/Test/sw_seal/cases/ directory.        
 
	% cd /nfs/3dma_rock/Test/sw_seal/cases/
	% runcase linux case2.14

	This produces a segmented file sw.gz in the sw_seal/c_seg/ directory
	which is padded with a layer of grain voxels on all sides of the
	imaged volume.  Note that sealed image dimensions are now
	appropriately larger (258x258x258 voxels in this case).

	Medial axis construction and trimming must be rerun on the sealed data6.
	Download the input files case2.5seal.in and case4.10seal.in into the
	/nfs/3dma_rock/Test/sw_seal/cases/ directory.
	Note: we will use the shell script "3dma_rock/bin/runcases" to run
	both of these cases consecutively. This script assumes that the input
	files begin with the letters "case".

	% runcases linux 2.5seal 4.10seal &

	The MA and the trimmed MA files are stored, respectively, in the
	sw_seal/ma/ and sw_seal/ma_t/ directories.

	Download the wedge based throat finding input file case4.12wdg_seal.in
	into the /nfs/3dma_rock/Test/sw_seal/cases/ and run it.

	% runcase linux case4.12wdg_seal

	This produces a throat file, sw_seal/throats/sw_wdg.gz which records
	the throat information for successful paths.

	Download the input files case5.5wdg.in, case4.12dij_seal.in and
	case5.6mrg.in into the /nfs/3dma_rock/Test/sw_seal/cases/ directory.

	% cd /nfs/3dma_rock/Test/sw_seal/cases/
	% runcase linux case5.5wdg

	This reads sw_wdg.gz and produces a file 'ma_no_thr_wdg' listing the
	indices of those paths for which the wedge algorithm failed (found no
	throat). This file will provide path information for running the
	Dijkstra algorithm.

	% runcase linux case4.12dij_seal &

	This runs the Dijkstra throat finding only on a subset of paths,
	and produces the throat file sw_seal/throats/sw_dij.gz

	% runcase linux case5.6mrg

	This merges the files sw_wdg.gz and sw_dij.gz, producing the merged
	throat file, sw_seal/throats/sw_mrg.gz, and an empty plot file
	sw_seal/results/cmpr.jgr. (This plot file is empty as the two throat
	files, sw_wdg and sw_dij, have no MA paths in common).
	
	Download case5.5mrg.in, case4.12sub_seal.in and case5.6.in
	into the /nfs/3dma_rock/Test/sw_seal/cases/ directory.
 
	% runcases 5.5mrg 4.12sub_seal 5.6 &

	case5.5mrg will produce a file 'ma_no_thr_mrg' listing the MA paths
	for which no throat exists in the file sw_seal/throats/sw_mrg.gz.
	Case4.12sub_seal will run the subsectioning Dijkstra algorithm
	on these paths, producing the throat file sw_seal/throats/sw_sub.gz. 
	Case5.8 will merge the throat files sw_seal/throats/sw_mrg.gz and
	sw_seal/throats/sw_sub.gz into the final throat file sw_seal/throats/sw.gz.

	Because the sealed data set is two voxels larger in each direction,
	voxel indices no longer align correctly with the unsealed data.
	The data in the sw_seal medial axis and throat files need to be
	un-padded so we can use them on the original image.

	Note: these un-padded data sets are going to be stored in the sw directory space.
	The un-padded MA data sets are put into a new subdirectory, sw/ma_t_unpad.
	The un-padded throat data is put in the sw/throat directory.

	Download case4.17.in and case5.7.in into the
	/nfs/3dma_rock/Test/sw/cases/ directory.

	% cd /nfs/3dma_rock/Test/sw/cases/
	% mkdir ../ma_t_unpad
	% runcases linux 4.17 5.7

	This produces MA cluster path files in the sw/ma_t_unpad directory and
	an unpadded throat file, sw/throats/sw.gz.

6 Because of the sealing layer of grains, the medial axis for the sealed data is different than for the unsealed data.


Case 5.1/5.2:   Throat visualization
                      (previous section)    (return to table of contents)    (next section)

Throats are plotted together with the medial axis and with or without the grains included. Throats can be viewed either slicewise (raster file format), or in 3D (Inventor format). Again the continuity of the throats are difficult to follow in slice-by-slice viewing, but the relationship of the medial axis to the grain phase is clearer. 3D viewing is preferable, with the same problem that including the pore-grain surface hides much of the figure and makes it somewhat difficult to view the throats. Plotting without the grains has the downside that it makes perspective viewing difficult.

We distinguish throat surfaces and throat barriers. A throat surfaces is a polygonal surface (zero volume) that spans the throat perimeter voxels, passing through the throat center on the medial axis. A throat barrier consists of those pores space voxels through which the throat surface passes as well as a subset of neighboring pore voxels. (The throat barrier is constructed to be a 6-connected block of pore voxels that "blocks" the channel.)

Test example:
	Download the input files case5.1ras.in, case5.1gmv_sub.in,
	case4.1gmv_sub.inand case5.2gmv_sub.in into the
	/nfs/3dma_rock/Test/sw/cases/ directory. Note that we will again
	use the command "runcases", which will run all of these cases
	consecutively.
	
	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcases linux 5.1ras 5.1gmv_sub 4.1gmv_sub 5.2gmv_sub &
	
	Case 5.1ras produces slicewise raster plots, sw.XXX.ras.gz,
	in the ../throats/ directory. The other three cases produce 3D
	plots in Geomview format. In these we plot only a subvolume of the
	data since plotting the entire data set would produce extremely "busy"
	pictures.
	Case 5.1gmv produces a plot file, ../plots/thr_barr_subvol.list, of
	the throat barriers in the subregion.
	Case 5.2gmv produces a file, ../plots/thr_surf_subvol.list, of the
	subregion throat surfaces.
	Case 4.1gmv_sub produces a file, ../plots/ma_subvol.list, of the
	subregion MA.

	The Geomview plots of the same subregion can be viewed separately,
	or merged by catenation.
	
	% cd ../plots
	% geomview thr_barr_subvol.list &
	
	% cat ma_subvol.list thr_surf_subvol.list > tmp
	% geomview tmp & 

	Catenated plots are shown in Fig. 12.
Fig. 12:   (left) A Geomview plot of the throat barriers (purple) and medial axis (rainbox colored) in a [0,99]x[0,99]x[0,99] voxel subvolume of the example data set. (right) The view of the medial axis and throat surfaces (white) in the same subvolume.



Case 5.3:   Pore-throat network construction
                      (previous section)    (return to table of contents)    (next section)

With throats constructed, the pore space is divided into pores separated by throat surfaces. Case 5.3 organizes the pores and throats into a network with identified neighbor relationships. More details on pore-throat network partitioning algorithm can be found in Prodanovic, Lindquist and Seright, 2005.

Test example:
	Download the input file case5.3dij.in into the
	/nfs/3dma_rock/Test/sw/cases/ directory.
	
	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case5.3dij
	
	This runs a pore network analysis based upon throats
	(../throats/sw_dij.gz) found only by the Diskstra-based algorithm.
	It produces the output file case5.3dij.out, throat files with
	reconstructed throat barriers, ../throats/sw_dij_cln.gz, and the
	pore-throat network data files, sw_dij_cln.th2np.gz and
	sw_dij_cln.np2th.gz, in the ../throats/ directory.

The two network data files cross reference the network: sw_dij_cln.th2np cross references throats to pores, while sw_dij_cln.np2th cross references pores to each other and to throats. Complete information on the structure of the *.np2th and *.th2np data files is found in the appendices of the 3DMA General Users Manual.

Note: If you have gone through the test example from the Aggressive Throat Computation section, you can run Case 5.3 on the set of throats determined by its more robust throat finding.

Test example using Aggressive Throat Computation data:

	Download the input file case5.3.in into the
	/nfs/3dma_rock/Test/sw/cases/ directory.
	
	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case5.3
	
	This runs a pore network analysis based upon throats
	(../throats/sw.gz) found by running the aggressive throat computation
	in the sealed image. (It uses the MA data in ../ma_t_unpad.) It
	produces the output file case5.3.out, the throat file with
	reconstructed throat barriers ../throats/sw_cln.gz and the
	corresponding pore-throat network data files, sw_cln.th2np.gz and
	sw_cln.np2th.gz, in the ../throats/ directory.
The difference in the resulting network obtained in using only throats found by the Dijkstra algorithm (case5.3dij.out) and those using the aggressive throat finding approach (case5.3.out) are summarized in the following table.

number
of pores
number
of throats
size of
largest pore
% pore volume (pv)
in interior pores
case5.3dij.out 431 1135 286514 voxels
(9.7% pv)
45.66
case5.3.out 619 1751 160344 voxels
(5.4% pv)
47



Case 6.1:   Surface plots of selected pores
                      (previous section)    (return to table of contents)    (next section)

After construction of the pore-throat network, one can visualize a selected set of pores, compute their area and/or gather statistics on the surface area of the selected set of pores. Each pore is delimited by grain bodies and by the throats that connect it to neighboring pores. Case 6.1 plots the bounding (grain + throat) surface for each pore of a selected set. The output is (set of) 3D Geomview plots.

Test example:
	Download the input file case6.1dij.in into the
	/nfs/3dma_rock/Test/sw/cases/ directory.
	
	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case6.1dij
	
	This produces the output file case6.1dij.out, pore surface area
	statistics in the plot file ../results/pore_area_dij.jgr,
	and the Geomview format files pore_49.list, pore_49_THR.list,
	pore_49_MA.list and pore_49_PORE.list in the ../plots/ directory. 
	File pore_49.list contains data on the pore-grain surface for
	pore 49; pore_49_THR.list contains data on the throat barriers for
	this pore; pore_49_MA.list contains data on the MA passing through
	this pore; pore_49_PORE.list contains data on all the voxels
	belonging to the pore.

	To view the area statistics file:
	
	% cd ../results
	% jgraph pore_area_dij.jgr > pore_area_dij.ps
	% ggv pore_surf.ps
	
	To view a 3D plot of the surface of throat 49:

	% cd ../plots
	% geomview pore_49.list pore_49_THR.list pore_49_MA.list &
	
	In the Geomview main window, click on 'pore_49.list" in the 'Target'
	list and then set the material to be transparent in the drop down
	menu (Inspect -> Material). Plots of the data on pore 49 are
	shown in Fig. 13 with and without transparency of the pore-grain
	surface.
	
Fig. 13:   (Left) Pore-grain surface (brown) of a pore (identified as pore 49) and its 4 adjacent throat barriers (voxels in white with black edges). (Right) A view of the same pore plotted with the pore-grain surface transparent. The medial axis (red) and branch cluster (blue) voxels inside the pore are now visible.



Case 6.2:   Pore-throat network statistics plots
                      (previous section)    (return to table of contents)    (next section)

This case produces a jgraph format output file containing the following distributions measured from either the throat or pore-throat network files.

Test example:
	Download the input file case6.2dij.in  into the
	/nfs/3dma_rock/Test/sw/cases/ directory.
	
	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case6.2dij
	
	This will produce the output file case6.2dij.out, and the jgraph
	summary statistics file ../results/np_thr_net.jgr.
	
	% cd results
	% jgraph < np_thr_net.jgr > np_thr_net.ps
	% ggv np_thr_net.ps
	
	A plot of np_thr_net.jgr is presented in Fig. 14.
Fig. 14:   Distribution plots produced by Case 6.2dij



Fluid analysis:   Directory setup and example data download for multiple fluid phase studies
                      (previous section)    (return to table of contents)    (next section)

3DMA-Rock has been extended to analyze two or more fluids in the pore space, with particular emphasis on comparison of fluid distribution in the same region at two different times. Fluid displacements and X-ray microtomography (CT) have to be performed without removing the rock core from the CT imaging mount to allow for direct comparison of imaged regions.

In our work to date, the first image consists of a single fluid phase in the pore space. Subsequent images of the same region contain either two liquid phases, a liquid phase and a gas, or two liquid phases and a gas phase in the pore space. In the case of two liquid phases, one liquid phase is typically doped with an X-ray attenuating additive (e.g. iodine or bromine) in order to distinguish it from the other liquid phase. This usually reduces the contrast between the doped liquid phase and the grain phase. Thus the first image (containing only a single fluid in the pore space) is used to establish where the grain phase is. The grain locations from the first image are carried over to subsuquent images. (Hence an assumption of strict alignment of the images in the sequence.)

In this primer we use the sw data set above, which is water saturated, as the first image. For a second image, we use the data set 'sor', which contains water and oil in the same pore space. The oil phase is at residual conditions in the sor data set.

- Create the subdirectory /nfs/3dma-rock/Test/sor.
- Create a subset of the usual 3DMA-Rock subdirectories within Test/sor/ as well as extra fld_seg subdirectories to hold segmented data on the individual fluids. (In this case we have two fluids, oil and water.)

  % cd /nfs/3dma-rock/Test/
  % mkdir sor/
  % cd sor
  % mkdir tomog cases results plots fld_seg fld_seg/ras fld_seg/h2o fld_seg/oil

- Use the following link to download the archived (sor-tomog.tar.gz) data set of tomographic input files into the /nfs/3dma_rock/Test/sor/tomog/ directory.

Warning: This archived data set is 14.4 MBytes in size. The ftp transfer may take some time.

The example contains 256 files labeled 'sor.XXX.gz' where XXX ranges from 001 to 256, containing 256 slices (the z-direction) of a 256x256x256 volume of a Berea sandstone core which contains both water and oil. Each slice is an array of 256^2 numbers (unsigned chars, values ranging from 0 to 255) representing X-ray attenuation coefficients. The oil used in flooding the Berea core was doped; in plots the oil will therefore appear darker than water (but not as dark as the Berea grains).

- Unzip and extract the tomographic files.

  % gunzip sor-tomog.tar.gz
  % tar -xvf sor-tomog.tar



Case 7.1:   Pore space histogram
                      (previous section)    (return to table of contents)    (next section)

If the "sw" and "sor" images are carefully aligned, the grain phase (e.g. Berea) voxels occupy the same position in both images. The segmented files from the "sw" image can therefore be used to identify the grain space of the "sor" image. The remaining voxels in the "sor" image are then pore voxels containing either oil or water phase.

In the example data set provided, the oil phase is doped with X-ray attenuating bromine, and has X-ray attenuation coefficients higher than water. In order to classify voxels in the pore space as oil or water, the indicator kriging algorithm is utilized. To generate the T0, T1 threshold values required by the kriging algorithm, Case 7.1 uses the sw segmentated data files to isolate the pore space in the sor data and compute the distribution of attenuation coefficients only in the sor data pore space.

Note: Case 7.1 differs from Case 1.6 in that Case 7.1 does not include any attenuation coefficient data from grain phase voxels!

Test example:
	Download the input file case7.1.in  into the
	/nfs/3dma_rock/Test/sor/cases/ directory.
	
	% cd /nfs/3dma_rock/Test/sor/cases/
	% runcase linux case7.1
	
	This produces the output file case7.1.out and the jgraph file
	sor/results/pore_hist.jgr.  The jgraph file is processed and viewed
	as usual.
	
	% cd ../results/
	% jgraph pore_hist.jgr > pore_hist.ps
	% ggv pore_hist.ps

	The resulting plot is shown in Fig. 15.

Fig. 15:  Attenuation coefficient distribution in the pore space of sample sor. Water and oil peaks are around attenuation coefficients 60 and 75 respectively; the peaks are not well separated.



Case 7.2:   Fluid segmentation via Indicator Kriging
                      (previous section)    (return to table of contents)    (next section)

Fluid segmentation options in Case 7.2 include simple thresholding and indicator kriging, both adapted to work only in the pore space. As shown above, simple thresholding produces segmented images of poorer quality; indicator kriging is preferred. The user input threshold values T0 and T1 for indicator kriging can be determined from the distribution obtained in Case 7.1.

Case 7.2 uses the segmented sw data to identify the pore space in the sor data and produces two segmented files for two fluid phases7 (oil and water in this case). Each segmented file is bitpacked, with the voxels of one fluid phase labeled with 0 and all other voxels labeled 1.


7 Thus to segment three fluid phases, Case 7.2 must be run twice, with two separate pairs of values for T0 and T1. The first pair of values separates fluid 1 from the other two. The second pair of values separates fluid 3 from the other two.
Test example:
	Download the input file case7.2.in
	into the /nfs/3dma_rock/Test/sor/cases/ directory.

	% cd /nfs/3dma_rock/Test/sor/cases/
	% runcase linux case7.2

	This will produce the output file case7.2.out, segmented volume
	files for both water (fld_seg/h2o/sor.gz) and oil (fld_seg/oil/sor.gz)
	as well as raster plots of the fluids for each slice in the files
	fld_seg/ras/sor.XXX.ras.gz.

	% cd ../fld_seg/ras
	% xv sor.*.ras* &

	An example raster plot is not shown here. Note that in the images,
	the water phase is green and the oil phase is red.



Case7.3:   Pore/throat fluid saturation analysis
                      (previous section)    (return to table of contents)    (next section)

We are interested in water saturation behaviour in individual pores and throats, specifically how saturation may be related to pore and throat size. Case 7.3 uses the pore/throat network obtained from the sw data set to identify individual pores and throats in the sor data and the sor fluid segmentation files to locate water and oil voxels. Case 7.3 outputs total water saturation in the output file (note that oil and water saturations add to 1.0) as well as a summary plot of saturation statistics organized by pore and throat size. For a throat, the saturation is computed for the set of voxels comprising the throat barrier.

Test example:
	Download the input file case7.3.in into the
	/nfs/3dma_rock/Test/sor/cases/ directory.

	% cd /nfs/3dma_rock/Test/sor/cases/
	% runcase linux case7.3

	This produces the output file case7.3.out and the summary plot file
	../results/fld_sat.jgr. To view the graph

        % cd ../results
        % jgraph fld_sat.jgr > fld_sat.ps
        % ggv fld_sat.ps

	The plots in fld_sat.jgr are shown in Fig. 16.
Fig. 16:  Distibution plots produced by Case 7.3. The first (top) and the fourth (bottom) scatterplots show water saturation vs. volume for pores and throats respectively. In addition, for each volume range, a box plot summarizes the distribution of measured saturations for all pores in that volume range. The second graph shows the fraction of all pores (dashed line) and throats (solid line) in each size range that have water saturations exceeding 0.5. For both pores (solid lines) and throats (dashed lines) the third graph shows two curves. For each respective geometric region, the upper curve of each pair plots the fraction of total pore (throat) volume occupied by all pores (throats) in each size range. The lower curve shows the fraction of total pore (throat) volume that is occupied by water in all pores (throats) of the specified size range. Pore values are given on the left axis, throat values on the right.



Case7.5:   Fluid blob distribution
                      (previous section)    (return to table of contents)    (next section)

Case 7.5 analyzes the sor fluid segmentation files to summarize the size distribution of disconnected blobs of both fluids in the pore space. The wetting fluid phase (water in the case of Berea sandstone) is treated as 26-connected; the non-wetting phase (oil in Berea) as 6-connected. This case is analogous to the action of Case 2.3 for pore/grain phases.

Test example:
	Download the input file case7.5.in into the
	/nfs/3dma_rock/Test/sor/cases/ directory.

	% cd /nfs/3dma_rock/Test/sor/cases/
	% runcase linux case7.5

	This produces the output file case7.5.out and the summary plot
	file ../results/fld_blob.jgr.  To view the graph

        % cd ../results
        % jgraph fld_blob.jgr > fld_blob.ps
        % ggv fld_blob.ps

	The plots in fld_blob.jgr are shown in Fig. 17.
Fig. 17:  Fraction of total pore volume occupied by water (green) and oil (red) by disconnected fluid blobs in each size range. Note that both x- and y- axis are plotted in log scale. (Note that in the jgraph file fld_blob.jgr, these two curves are plotted in separate graphs.)

A variation of the input file (case7.5a.in) is available if one wants to plot the distribution of pore volumes on the same ordinate and abcissa scales for comparison with the fluid blob distributions.



Case 7.8/7.9:   Fluid visualization
                      (previous section)    (return to table of contents)    (next section)

Cases 7.8 uses the sw pore/grain segmentation files and the sor fluid segmentation files to provides three different visualization options:

As these cases are straightforward to run and are similar to previous visualization cases discussed above, input files are provided without further "Test example" commentary. An example of the output from the Geomview plot (3D) capability (case7.8sub.in) is displayed in Fig. 18.
water-oil interface
Fig. 18:   Geomview plot of water-oil interface surface for a randomly chosen subvolume of size [0,86]x[0,62]x[21,153] produced by Case7.8sub.in.

For each image slice, Case 7.9 organizes a raster file showing side-by-side comparison of:   the sor tomographic image; the sor tomographic image with the grain phase (from the sw segmentation data) overlayed in a single color; and the results of the fluid segmentation with each fluid colored separately. (See Fig. 19.) The middle image can be used to visually check the grain/pore alignments between the first and subsequent CT images. It also helps to focus the eyes only on the pore space to visually check the accuracy of the fluid segmentation shown in the third image.

Test example:
	Download the input file case7.9.in into the
	/nfs/3dma_rock/Test/sor/cases/ directory.

	% cd /nfs/3dma_rock/Test/sor/cases/
	% runcase linux case7.9

	This will produce the output file case7.9.out, and comparative
	raster plot files ../plots/compare.XXX.ras.gz.

	To view the raster plot files with 'xv':

        % cd ../plots
        % xv compare* &

	An example slice is shown in Fig. 19.

comparative raster plot
Fig. 19:   Case 7.9 segmentation summary for slice 110:   (left) the tomographic data;   (center) tomographic data in the pore space view (grain phase colored yellow) - oil is darker than water;   (right) fluid segmentation results - grain, water, oil.

Note that the fluid segmented files (which have the voxels of one fluid marked 0 and all other voxels marked 1) can be also plotted using Case 2.1, which will then plot only this particular fluid or its surface.



Case 6.1:   Pore scale fluid visualization
                      (previous section)    (return to table of contents)    (next section)

Case 6.1 discussed earlier also has the capability to produce 3D plots of the fluid partitioning in individual pores (aside from pore surface plot/area computation). It uses the sw pore/grain segmented, throat and pore-throat network data files, and the sor fluid segmented data files.

Test example:
	Download the input file case6.1sor.in into the
	/nfs/3dma_rock/Test/sw/cases/ directory.
	
	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case6.1sor
	
	This will produce the output file case6.1dij.out, statistics on
	the pore surfaces in the file ../results/pore_area_dij.jgr and the
	following files in the ../plots/ directory:
	The Geomview plots can be viewed in any catenation combination. For
	example:

	% cd ../plots
	% geomview pore_49.list pore_49_OILintfc.list pore_49_THR.list &

	Click on 'pore_49.list' in the 'Target' list in the Geomview main
	menu and set material to be transparent in the drop down menu
	(Inspect->Material). Click on 'pore_49_OILintfc.list' in the
	'Target' list and then choose Inspect->Appearance in the drop
	down menu. Change the color of the faces to be red (i.e. RGB 1.0 0 0).
	The result should be similar to Fig.20 (left). When rotating (on
	Tools menu) make sure you click on 'World' in the 'Target' list.

	% geomview pore_49.list pore_49_OIL.list pore_49_THR.list &

	With the pore/grain surface set to transparent and the oil voxels set
	to red, the result is Fig. 20 (right).
Fig. 20.   (Left) Residual oil fluid surface (red) within pore 49 (transparent brown). Throat voxels attached to the pore are shown in light grey.   (Right) Actual residual fluid voxels are shown in red. (cf. Fig. 13.)



Case 8.1:   Lattice Boltzmann single phase flow simulation of throat permeabilities
                      (previous section)    (return to table of contents)    (next section>

3DMA-Rock has the capability of simulating single phase flow via both the lattice Bhatagnar-Gross-Krook (LBGK) and multiple-relaxation-time (MRT) models. Both models compute 3D single phase flow in various geometries ( flow in a rectangular channel, flow in the pore space of a segmented image of a porous medium, flow through a single pore or through a single throat extracted from an imageimages) utilizing the D3Q19 lattice. With somewhat reduced capabilities it can also compute 2D flow utilizing the D2Q9 lattice. For more details refer to the thesis of Dr. M. Prodanovic and references therein.

Note that the code is developed for a single processor. While it can be used for calculation of a rock sample absolute permeability, the simulation might take days even for a 1283 volume.

We have developed the code primarily for calculation of individual throat (as captured in a digitized rock image) permeabilities; this information is usually needed by network flow simulators.

Test example 1 (permeability calculation and visualization for a single throat):

	Download the input file case8.1thr10.in into the
	/nfs/3dma_rock/Test/sw/cases/ directory.

	(optional) In the 3DMA-Rock code file 'src/latt_boltz/latt_boltz.h'
		set the variable 'DBG_PLOT' to 1 and recompile the code. Create
		the subdirectory sw/dbg_plots.

		% cd /nfs/3dma_rock/Test/sw
		% mkdir dbg_plots
	
	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case8.1thr &
	
	This will produce the output file case8.1thr.out and a Geomview
	visualization of individual voxel velocities in ../plots/veloc_thr10.list.
	Calculated permeabilities can be found in the output file.
	
	% cd ../plots
	% geomview ../plots/veloc_thr10.list &

	See Fig. 21 b,d for views of the velocity field through this throat.
	
	(optional) If the DBG_PLOT option was turned on (above), this
		will cause a number of debugging Geomview files to be output
		in a separate directory.

		% cd ../cases/dbg_plots
		% geomview pdat_thr10_THR.list nodes_thr10.list &
	
		In this example, these debugging plots show the voxels in the
		throat channel that are involved in the flow computation.
		Fig 21 a,c shows views of the voxels for throat 10.

a) b)
  
c) d)

Fig. 21. a)   Computational domain visualization (available as a debugging option) for throat 10. Original barrier voxels are in violet, the sides of inlet boundary voxels are in bright green and the sides of outlet boundary voxels are bright blue. b) Velocity vectors for all computational domain voxels. The rainbow coloring reveals the spectrum of absolute velocity values (the smallest in red, the largest in violet). c) and d) are rotated views of a) and b).
Test example 2 (permeability calculation only for all throats in the pore-throat network):

	Download the input file case8.1thr.in into the
	/nfs/3dma_rock/Test/sw/cases/ directory.
	
	(optional) in the 3DMA-Rock code file 'src/latt_boltz/latt_boltz.h',
		set the variable 'DBG_PLOT' back to 0 and recompile the code.
		(Leaving the DBG_PLOT option set to 1 will create way too many
		output files!)
	
	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case8.1thr &
	
	This will produce the output file case8.1thr.out. Extensive
	output is printed out for each throat and crucial information
	for all throats is assembled into a table at the end of the
	file.



Case2.17: Create a segmented image of a hexagonal close packing of spheres
                      (previous section)    (return to table of contents)    (next section)

3DMA-Rock has the capability of creating a segmented data set of a digitized hexagonal close packings of spheres. Once a segmented file is created, one can proceed with porous structure analysis as described above - medial axis, throat finding....

Test example
	Create a new directory to store the sphere pack data and
	analysis subdirectories

	% cd /nfs/3dma_rock/Test/
	% mkdir hcp
	% cd hcp
	% create_wdir
	
	Download the input file case2.17.in into the
	/nfs/3dma_rock/Test/sw/cases/ directory.
	
	% cd cases
	% runcase linux case2.17
	
	This creates a segmented file, seg/hcp.gz, of 1283 voxels
	with spheres of diameter 32 voxels closely packed on a hexagonal
	grid. Furthermore, case2.17.out provides information on both
	analytical spheres and their digital models, including
	analytically predicted throats. Inventor plots of analytical and
	digitized spheres are produced in seg/hcp_anal.iv.gz and
	seg/hcp_dig.iv.gz. These can be viewed by:
	
	% cd ../seg
	% gunzip hcp_anal.iv.gz hcp_dig.iv.gz
	% ivview hcp_anal.iv
	% ivview hcp_dig.iv

	The contents of these plots are shown in Fig. 22.

Fig. 22.   Inventor plots of a hexagonal close packing of spheres (left) and their digitized surfaces (right) within a 1283 voxel volume.




Appendix A: Segmented Volume File Format
                      (previous section)    (return to table of contents)    (next section)

One segmented file is printed for the entire data volume. In this file, each voxel is labeled either 1 or 0. The segmented information is therefore stored in bit-packed format, the values for each set of 8 successive voxels are stored in an unsigned char (byte).

The data is stored slice-wise (z-direction). In general the number of voxels in a slice may not be divisible by 8 and the last unsigned char needed for bit packing may not be full. The header for the bit packed file therefore records the number, bnxyz, of unsigned chars used excluding the last, and the number, nresid, of remaining bits stored in the the final unsigned char. Note 1≤nresid ≤ 8. The total number of unsigne chars used to store the segmented volume is therefore bnxyz +1, plus a 28 byte (7 integer) header.

The format for the entire segmented file, including the header, is:




Further Reading
                      (previous section)    (return to table of contents)

SUNYSB-AMS-05-04:   3D image-based characterization of fluid displacement in a Berea core, Prodanovic, Seright and Lindquist.
SUNYSB-AMS-05-13:   Porous structure and fluid partitioning in polyethylene cores from 3D X-ray microtomographic imaging, Prodanovic, Seright and Lindquist.
SUNYSB-AMS-05-15:   Computed microtomography studies of fluid partitioning in drainage and imbibition before and after gel placement: disproportionate permeability reduction, Seright, Prodanovic and Lindquist.
SUNYSB-AMS-03-18:   Analysis of the vesicular structure of basalts, Shin, Lindquist, Sahagian and Song.
SUNYSB-AMS-03-01:   Synchrotron X-ray computed microtomography (CMT) studies of vesiculated basaltic rocks, Song, Jones, Lindquist, Dowd and Sahagian.
SUNYSB-AMS-02-03:   Use of X-ray computed microtomography to understand why gels reduce permeability to water more than to oil, Seright, Liang, Lindquist and Dunsmuir.
SUNYSB-AMS-02-02:   Characterizing disproportionate permeability reduction using synchrotron X-ray computed microtomography, Seright, Liang, Lindquist and Dunsmuir.
SUNYSB-AMS-01-14:   Network flow model studies and 3D pore structure, Lindquist.
SUNYSB-AMS-01-06:   Quantitative analysis of three dimensional X-ray tomographic images, Lindquist.
SUNYSB-AMS-99-20:   3DMA General Users Manual, Lindquist.
SUNYSB-AMS-99-13:   Pore and throat size distributions measured from synchrotron X-ray tomographic images of Fontainebleau sandstones, Lindquist, Venkatarangan, Dunsmuir and Wong.
SUNYSB-AMS-98-02:   Image thresholding by indicator iriging, Oh and Lindquist.
SUNYSB-AMS-98-01:   Invesigating 3D geometry of porous media from high resolution images, Lindquist and Venkatarangan.
SUNYSB-AMS-95-01:   Medial axis analysis of three dimensional tomographic images of drill cores, Lindquist, Lee, Coker, Jones and Spanne.


Authors: Masa Prodanovic, Brent Lindquist