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.

• Introduction and motivation
• Directory set-up, download example input files
• Quick look-up:
  • Running the code
  • List of input files
• Operations performed on tomographic data
  • Case 1.1: Plotting the grey scale image
  • Case 1.6: Image intensity histogram, porosity-threshold dependence
  • Case 1.5: Segmentation
    • Case 1.5(smpl): Simple segmentation
    • Case 1.5(krig): Indicator kriging
• Operations performed on segmented data
  • Cases 2.3/2.4: Post-segmentation clean-up procedures
  • Case2.1: Visualization of segmented files and interphase surface area measurement
  • Case 2.5: Burning and generating the medial axis
• Operations performed on and using medial axis data
  • Case 4.10: Medial axis trimming (clean-up)
  • Case 4.1: Medial axis plotting
  • Case 4.6: Medial axis path geometric tortuosity
  • Case 4.12: Throat computaton
  • Aggressive throat computation for medium-high porosity media (NEW)
• Operations performed on throat data
  • Cases 5.1/5.2: Throat visualization
  • Case 5.3: Pore-throat network construction
• Operations performed on pore network data
  • Case 6.1: Individual pore surface area calculation and plots
  • Case 6.2: Pore-throat network statistics plots
• Fluid analysis: two or more pore space fluids
  • Directory setup and example data download
  • Case 7.1: Pore space histogram
  • Case 7.2: Segmentation of fluids in the pore space
  • Case 7.5: Fluid blob size distribution
  • Case 7.8/7.9: Fluid visualization
  • Case 7.3: Pore-by-pore water saturation analysis
  • Case 6.1: Pore scale fluid visualization (NEW)
    
• Lattice Boltzmann method flow simulation (NEW)
  • Case8.1: 3D, single phase, throat permeability computations
Additional topics
• Case 2.17: Create a segmented file of a hexagonal close packing of spheres (NEW)
• Appendix : Segmented volume file format
• Further reading

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

- 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

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

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

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.

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

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.

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.

The indicator kriging method requires input of two thresholds, T₀ and T₁. T₀ (T₁) 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 [T₀,T₁]) 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).

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 1² can be especially troubling if medial axis analysis of phase 0² 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 histograms³ 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.

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.

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

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.

Fig. 6:  Geomview graphics view of the pore-grain voxel surfaces of the same subvolume.

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.

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:

• remove isolated components corresponding to isolated pores. These come in two types, isolated paths, and single isolated medial voxels (paths of length 1 voxel). A path length threshold is used to determine what isolated paths to keep, a distance measure (burn number) threshold is used for isolated voxels.
• remove dead end (branch-leaf) paths. Decision parameters are based upon path length and whether a path exits through the boundary. For pore studies, we recommend removing all branch-leaf paths except those exiting the boundary.
• needle eye-paths are paths which loop, i.e. return to the same branch point in the network from which they left. Unless it is known that there are closed loop paths in the pore space, all occurances of needle eyes are probably small path segments produced by local segmentation effects and should be removed.
• Reduce "surface remants": The medial axis consists of a network of paths that join at "branch points". In a digitized medium, a "branch point" may in fact not be a single voxel, but a set of voxels, which we refer to as branch clusters. Ideally the set of voxels forming the branch cluster has genus 1 (one object, no holes or embedded cavities). In fact, it is possible to produce branch clusters that have holes. We refer to such clusters as "surface remnants", a phrase which is somewhat of a misnomer. 3DMA-Rock can attempt to simplify the structure ot branch clusters, reducing them to genus 1 objects if possible (not guaranteed).
• Merge close clusters: one goal of the application of medial axis to the study of pore structure is to identify those medial axis paths that connect separate pores, as opposed to paths which lie entirely in a single pore. The "merge close cluster" algorithm attempts to identify paths that remain inside a single pore, essentially by performing one of a number of test variants that compare the path length to the distance "from the path to the grain surface". If a path is judged to lie in a single pore, it is "removed" from the path list by "merging" it and the two branch clusters it joins into a "super cluster". No actual change in medial axis voxels occurs, this is merely a book-keeping device.

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/

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.

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.

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

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.

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.

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

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.

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:

• Boundary effects. Throat algorithms can often not find a throat in a channel cut by the image boundary, due to lack of information beyond the imaged region. To avoid this issue one can "seal" the image volume by padding it with an external layer of grain voxels and perform throat finding in the sealed image.
• Too few throats found by each of the throat algorithms. Each throat finding algorithm has a certain failure rate⁵. For every throat not found, two pore bodies are merged. Even a relatively small failure rate can generate large pore bodies, especially if the failures predominate in a single region. We exploit the fact that the "failure regimes" are different for the different algorithms. Let P denote the set of all medial axis paths. The following analysis sequence employs all available throat finding algorithms to maximize throat finding while controlling the computational time to do so.


• Compute throats by the fast wedge algorithm. Let P_W denote the set of paths with successfully identified throats by this algorithm.
• Run the Dijkstra throat algorithm on the set P/P_W, i.e. paths for which the wedge algorithm found no throat. Let P_WD denote the set of paths that have had throats identified in these two steps.
• Run the subsectioning algorithm on the set P/P_WD of medial axis paths.
• Merge all the throats found in these three steps into one throat file.

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.

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.

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.

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.

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

• Pore body volume distribution.
• Throat surface area distribution.
• Effective pore radii distribution. The effective radius is the radius of a sphere of equivalent volume.
• Effective throat radii distribution. The effective radius is the radius of the circle of equivalent area.
• Joint pore body/pore body size scatterplot. For each pore, the average of the effective radii of its neighboring pore bodies is computed and plotted vs. the effective radius of the pore at hand.
• Joint throat/pore body size scatterplot. For each throat, the average effective radius of the two adjoining pore bodies (if they both have been found) is computed and plotted vs. the effective throat radius.
• Pore body coordination number/pore body size scatterplot. The coordination number of a pore is the number of neighboring pore bodies, and is plotted vs. the effective pore radius.
• Effective throat/pore radius ratio distribution. For each throat, the average of two (different) adjoining pore bodies (if they exist) is computed. This plot gives the distribution of the ratio of the effective throat radius to the computed average pore body radius.

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

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

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.

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 phases⁷ (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.

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.

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.

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.

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

• Rasterfile plot (2D) - provides slicewise raster plots of the contents of the fluid segmented files (case7.8ras.in).
• Geomview plot (3D) - plots the contact surface between the two fluids (using the Marching Cubes algorithm) and calculates the area of this interfluid surface. (case7.8sub.in).
• Total fluid/fluid interface area computation only - This option avoids producing a (potentially very big) Geomview file and only computes the fluid interfacial area (for the entire image) (case7.8area.in).

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.

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

• pore_49.list (the pore/grain surface (via Marching Cubes)),
• pore_49_THR.list (the throat barriers for the pore),
• pore_49_OIL.list (oil voxels within pore),
• pore_49_OILintfc.list (surface between the oil voxels and all others (via Marching Cubes)),
• pore_49_H2O.list (water voxels within the pore)
• pore_49_H2Ointfc.list (surface between the water voxels and all others (via Marching Cubes)),
• pore_49_intfc.list (surface between the water and oil voxels within the pore (via Marching Cubes))
• pore_49_PORE.list (all voxels within the pore).

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

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

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.

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, b_nxyz, of unsigned chars used excluding the last, and the number, n_resid, of remaining bits stored in the the final unsigned char. Note 1≤n_resid ≤ 8. The total number of unsigne chars used to store the segmented volume is therefore b_nxyz +1, plus a 28 byte (7 integer) header.

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

• endian (binary character). 'e' designates that the file was written on little endian hardware; 'E' designates that the file was written on big endian hardware
• n_x (binary integer). number of columns of voxels in each slice
• n_y (binary integer). number of rows of voxels in each slice
• z_s (binary integer). identification number for the first slice stored in the file (convention is to start from 1)
• z_e (binary integer). identification number for the last slice stored in the file
• b_nxyz (binary integer). number of full unsigned chars needed to store the volume
• n_resid (binary integer). number of remaining bits stored in the final unsigned char
• An array of b_nxyz + 1 unsigned chars containing the bit-packed segmentation information for the volume