3DMA-Rock Primer

This page assumes you have obtained a copy 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.

The 3DMA General Users Manual, while a more complete reference, might be a bit overwhelming, and has not been updated for new features. Here are examples of the basic routines ('cases' in our jargon) to be run to produce an analysis of the 'pore' space in a bi-phase medium.



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

- Create a directory for your data analysis. We will use an 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
  mkdir tomog cases burn seg c_seg ma ma_t plots results throats

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

tomog       raw, tomographic files
cases       input, output files
burn        storage of erosion numbers
seg         segmented files
c_seg       storage of 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 100 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



Execution
                      ( 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 profided 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

As mentioned, 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)

Input file   Function
case1.1.in   slicewise rasterplot of tomographic data
case1.1a.in   slicewise rasterplot of tomographic data with 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
case2.1ras.in   slicewise rasterplot of segmnented data
case2.1gmv_sub.in   subvolume plots of triangulated surface separating pore-grain surface via Geomview
case2.1gmv_vox.in   subvolume plots of contact surface between pore and grain voxels via Geomview
case2.3.in   volume information on disconnected pores and grains
case2.4.in   "clean-up" operations on segmented da
case2.5.in   erosion computation of media axis (Lee, Kasyap, Chu [1994] algorithm)
case4.10.in   trimming operations on medial axis
case4.1ras.in   slice-by-slice rasterplots of medial axis
case4.1gmv.in   volume plots of medial axis via Geomview
case4.1gmv_sub.in   subvolume plots of medial axis via Geomview
case4.6x.in   Medial Axis tortuosity computation (in x-direction)
case4.9.in   conversion from cluster-path to voxel Medial Axis file format
case4.12dij.in   computation of pore throats using Venkatarangan [2000] algorithm
case4.12wdg.in   computation of pore throats using Shin [2002] algorithm
case5.1ras.in   slice-by-slice views of throats via rasterplots
case5.1gmv_sub.in   subvolume plots of throat regions via Geomview
case5.2gmv_sub.in   subvolume plots of throat surfaces via Geomview
case5.3dij.in   construction of pore-throat network
case6.1.in   individual pore surface plots via Geomview; pore surface area computation
case6.2dij.in   output of pore-throat statistics
case7.1.in   histogram of tomographic data in pore space
case7.2.in   fluid segmentation (Indicator Kriging segmentation in pore space)
case7.3.in   fluid saturation in pores/throats
case7.5.in   fluid blob size (volume) distribution
case7.7.in   slicewise rasterplot of fluid segmented files
case7.8sub.in   subvolume plot of water-oil interface surface via Geomview
case7.9.in   slicewise comparative rasterplot of segmented, tomographic and fluid segmented data



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 rasterfiles (2-D plots of each slice in the
	stack which are 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 rasterfiles
	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 resulting 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 jgraph file can be processed and viewed as follows:

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

The jgraph program is publically available software that is distributed along with 3dma code. It produces postscript output which can be printed, included in documents, or viewed by a postscript interpreter (e.g. ghostscript). Fig. 2 shows the resulting plots obtained in the above example output file poros_thres.ps.

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

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 (Berea rock core).



Case 1.5smpl:   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 rasterfiles 
          sw_smpl.XXX.ras.gz also in the ../seg/ directory.

	  The rasterfile provide a graphical of the segmentation 
	  results for each slice and just echo** the output of 
	  the sw_smpl.gz (volume) data file.

	% cd ../seg/
	% xv sw_smpl.*.ras.gz

**So why print both the data files sw_smpl.gz and the rasterfiles 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.5krig:   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 also 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 image for the 110th slice in the test example, with the addition of image histogram.   (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 sample histogram with  specifically marked indicator kriging thresholds (cut-off values).

With the image segmented, all further processing steps utilize the segmented data volume file. Complete information on the structure of these files is found in the appendices of the 3DMA General Users Manual.



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 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, case2.3 should be run first. Case 2.3 provides both text printout and jgraph histogram 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 grain "tip". In such cases, we consider it 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 will produce the output file case2.3.out and a jgraph
	  file disc_vol.jgr in the /nfs/3dma_rock/Test/sw/results directory.
	  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.
	
	  To delete ALL blobs but the largest of each phase, volume thresholds
	  (e.g. 32000 for phase 0 (water) and 1000 for phase 1 (Berea) can be
	  established simply by examining this printout. 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 these in the /nfs/3dma_rock/Test/sw/c_seg/
	  directory.
	
	% cd ../c_seg/
	% xv sw.*.ras.gz

**In the output for Case 2.3, the word "pore" refers to phase 0 voxels (for the example data provided, these are in fact water phase), while the work "grain" refers to phase 1 (for the example data provided, this is the sandstone grain phase).

There are other options in case2.4.in that are somewhat self-explanatory. To assist in their understanding, a case 2.4 input file with extra commentary can be accessed here.

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



Case 2.1:   Visualization of segmented files
                      ( previous section)    ( return to table of contents)    ( next section)

2D slice-wise rasterplots of segmented files have been encountered in the previous sections. A handy input file that will provide rasterplots 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 rasterplots) 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.

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 will produce the output files case2.1gmv_sub.out, case2.1gmv_vox.out 
	  and Geomview files in seg_surf_sub.list and seg_surf_vox_data.geom.gz
	  (zipped file) in plots/. See Figs. 5 and 6. Note that
	  case2.1gmv_sub.out contains the area (in units of voxel length -
	  typically microns) of the surface.

        % cd ../plots
        % geomview seg_surf_sub.list
        % gunzip seg_surf_vox_data.geom.gz
        % geomview seg_surf_vox_data.geom
pore-grain surface
Fig. 5:  Geomview graphics view of pore-grain surface of an [0,99]x[0,99]x[0,99] voxel subvolume of the sample obtained by the Marching Cubes Algorithm
  
contact voxel surfaces
Fig. 6:  Geomview graphics view of pore-grain contact voxel surfaces of an [0,99]x[0,99]x[0,99] voxel 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 will produce 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. 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 media axis, case4.11 provides the user with the ability to trim features from the medial axis. The major components of this trimming are:

Output from case4.11 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 will produce the output file case4.11.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 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)

Medial axis can be viewed either slicewise (rasterfile 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 graphic plot ma.list in the plots/ directory.
	  A view of this 3D plot is shown in Fig. 8.
Fig. 7:  A rasterplot 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. 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)

Assume a (trimmed) medial axis has been computed for the pore space inside a rectangular volume and output in MA voxel list format. (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.6x.in and case4.9.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 path list format 
          needed by case4.6. The resulting files will be "../ma_t/sw.XXX.gz"
 
	% runcase linux case4.6x
        % cd ../results 
	% jgraph < tort_x.jgr > tort_x.ps
	% ggv tort_x.ps

	- A plot of tort_x.ps is presented in Fig. 11.

	% geomview ../plots/tort_x_abs 

        - This will plot the shortest paths from face (3) to (4).
        '../plots/tort_x_all' contains all of the paths (the two are
	superimposed in Fig. 10).
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 two throat computation algorithms. One, developed as the thesis of Dr. A. Venkatarangan, utilizes Dijstra 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.

Note: The Dijkstra shortest path algorithm is the most computationally expensive algorithm in 3DMA-Rock. 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.

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.gz** in the ..throats/ directory containing
	  information on throats found. This case on the given sample will run
	  approximately 6 hours on a 2.0GB RAM, 1.7GHz Intel-Pentium processor.

	% runcase linux case4.12wdg

	- This produces a file sw_wdg.gz** in the ../throats/ directory containing
	  information on throats found.

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



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 view either slicewise (rasterfile 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 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 regions. Throat surfaces are polygonal surfaces that span throat center on medial axis and throat perimeter voxels which belong to grain phase. Throat barrier (region) consists of actual pores space voxels through which the throat surface passes.

Test example:
	- download the input files case5.1ras.in , case5.1gmv_sub.in ,
	case4.1gmv_sub.in and case5.2gmv_sub.in into the
	/nfs/3dma_rock/Test/sw/cases/ directory. Note that we will use
	the command "runcases", which will consecutively run all of these
	cases.

	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcases linux 5.1ras 5.1gmv_sub 4.1gmv_sub 5.2gmv_sub

	- Case 5.1ras produces raster plots sw.XXX.ras.gz in the ../throats/
	  directory. Case 5.1gmv produces a Geomview graphic plot
	  'thr_barr_subvol.in" of throat regions in the plots/ directory.
	  We plot only a subvolume of the original sample as plotting the
	  entire sample would produce an extremely "busy" picture. 
          Case 5.2gmv will produce an equivalent plot of throat surfaces in
	  the output file "../plots/thr_surf_subvol.list". As there is no
	  option for plotting medial axis in case5.2, we plot MA with
	  Case 4.1gmv_sub which produces "../plots/ma_subvol.list".

        % cd ../plots
        % geomview thr_barr_subvol.list &

        - The following will concatenate the graphic files and ensure proper
	superposition when viewing in Geomview: 

        % cat ma_subvol.list thr_surf_subvol.list > tmp
        % geomview tmp & 
Fig. 12:  (left) A Geomview 3D plot of the throat regions (purple) and medial axis (rainbox colored) in a [0,99]x[0,99]x[0,99] voxel subvolume of the example data set. A throat region is defined as the set of voxels through which the computed throat surface passes. (right) The eqivalent view of the throat surfaces (white).



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 explicitly specified neighboring relationships.

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 will produce the output file case5.3dij.out, and intermediate
	  data files sw_dij.th2np.gz and sw_dij.np2th.gz in the throat/
	  directory.

The two files cross reference the network: sw_dij.th2np cross references throats to pores, while sw_dij.np2th cross references pores to each other and to throats.



Case 6.1:   Individual pore surface plot
                      ( previous section)    ( return to table of contents)    ( next section)

After the pore-throat network construction, we might want to visualize individual pore, compute its area and/or gather statistics on
areas of all pores. Each pore is delimited by grain and by throats that connect it to the neighboring pores. Case 6.1 will plot closed surface of specified pore(s). The output is 3D Geomview plot.

Test example:
	- download the input file case6.1.in into the
	/nfs/3dma_rock/Test/sw/cases/ directory.

	% cd /nfs/3dma_rock/Test/sw/cases/
	% runcase linux case6.1.in

	- This will produce the output file case6.1.out, statistics on
	the pore surfaces in the file pore_surf.jgr in the ../results/
	directory and files pore_surf_51.list and pore_surf_51_THR.list
	in the ../plots/ directory. 
 
        % cd ../results
        % jgraph pore_surf.jgr > pore_surf.ps
        % ggv pore_surf.ps
        
        % cd ../plots
        % geomview pore_surf_51.list pore_surf_51_THR.list
Fig. 13:  Different rotational views of the pore (identified as pore 51) and its 4 adjacent throat regions (voxels in white,voxels sizes scaled to improve visibility). An isolated group of white voxels does not belong to pore 51 but merely resides in its bounding box. The subvolume of the original volume that contains the pore is [135,157]x[96,119]x[0,31] voxels.



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

This case produces various distributions measured from either the throat or pore-throat network files. The output file is in jgraph format This file contains:

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 file
	  np_thr_net.jgr in the ../results/ directory.

	% cd results
	% jgraph < np_thr_net.jgr > np_thr_net.ps
	% ggv np_thr_net.ps

	- A plot of np_thr_net.ps is presented in Fig. 14.
pt-net statistics
Fig. 14:  Distributions resulting from running Case 6.2dij



Multiphase 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 emaphsis on comparison of fluid distribution in the same region at two different times. Fluid flow and X-ray scanning has to be performed without removing the rock core from the imaging mount to allow for direct comparison of imaged regions. One fluid is typically doped with an X-ray attenuating additive (e.g. iodine or bromine) in order to distinguish the fluids in the image.

Create a subdirectory   Test/sor/ in /nfs/3dma-rock/Test/
Create the usual 3DMA-Rock subdirectories within Test/sor/

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

All directories will have the same function as before with the exception of ../fld_seg/ which will be used to store fluid segmented files.

- Use the following link to download an archived (sor-tomog.tar.gz) data set of example tomographic input files into your /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 100 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. Oil used in flooding Berea core was doped and will therefore in images appear darker than water (but not as dark as Berea).

- 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 computes a histogram of the distribution of attenuation coefficients only in the pore space. Note Case 7.1 differs from Case 1.6 in that it does not histogram 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 will produce the output file case7.1.out and the jgraph file
	  pore_hist.jgr in the directory /nfs/3dma_rock/Test/sor/results
	  The jgraph file can be processed and viewed as follows:

	% cd ../results/
	% jgraph < pore_hist.jgr > pore_hist.ps
	% ggv pore_hist.ps
Fig. 15:  Pore space histogram. Water and oil peaks are around attenuation coefficients 60 and 75 respectively - populations of water and oil are not well separated.



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

Fluid segmentation options include simple thresholding and indicator kriging, both adapted to work only in the pore space. As shown above, simple thresholding, gives segmented images of poor quality and the indicator kriging method is prefered. The user input threshold values T0 and T1 for indicator kriging can be determined from histogram obtained in Case 7.1.

Case 7.2 produces oil and water fluid segmented files - bitpacked files that contain 0 in place of particular fluid voxels and 1 in place of all other voxels.

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 and oil(w_sor.gz and o_sor.gz) in the ../fld_seg/
          directory as well as rasterplot of fluids for each slice - 
          sor.XXX.ras.gz in the ../fld_seg/ras/ directory.

        - To view rasterplot files, say with 'xv'    , do the following:
          Note that water phase is green and oil phase is red.

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



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

We are interested in water saturations behaviour in individual pores and throats, specifically related to pore/throat size. Case 7.3 will output individual water saturation in the output file (note that oil and water saturations add to 1.0) as well as a couple of plots. If one is analyzing a sequence of images, these plots may point out how fluid moves in pores/throats of certain size as well as if there is any difference in behaviour in pore and throat regions.

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 will produce the output file case7.3.out and fld_sat.jgr in 
          the ../results directory. To view the graph:

        % cd ../results
        % jgraph fld_sat.jgr > fld_sat.ps
        % ggv fld_sat.ps
Fig. 16:  Distibution plots produced by Case 7.3. The first (uppermost) and the fourth (lowermost) 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 will summarize the size distribution of disconnected blobs of both fluids in the pore space. (The "water" phase is treated as 26-connected; the "oil" phase as 6-connected.) This case is analagous to the action of Case 2.3 for pore/grain phases. For a time series of images, disconnected blob analysis shows bulk fluid movement.

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 will produce the output file case7.5.out and fld_blob.jgr in 
          the ../results directory. To view the graph:

        % cd ../results
        % jgraph fld_blob.jgr > fld_blob.ps
        % ggv fld_blob.ps
Fig. 17:  Fraction of total pore volume occupied by water (upper graph) and oil (lower graph) by disconnected fluid blobs in each size range. Blob size is expressed as volume; blob volumes are normalized to the total pore volume.



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

Cases 7.7-9 provides various options of visualizing the fluid segmentation in the pore space.

Case 7.7 (case7.7.in) provides slicewise rasterplots of the contents of the fluid segmented files.

Case 7.8 plots a triangulated contact surface between two fluids (case7.8sub.in). (See Fig. 19.)

For each image slice, Case 7.9 organizes a rasterfile showing side-by-side comparison of: the tomographic image; the tomographic image with the one phase overlayed in a single color; and the results of the fluid segmentation with each fluid colored separately. (See Fig. 18.)

Finally, note that fluid segmented files are as any other segmented files that have particular fluid marked as 0 and everything else as 1. Therefore, in order to plot only particular fluid or it's surface, one can use options in Case 2.1.

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
	  rasterplot files compare.XXX.ras.gz in the plots/ directory

        - To view the rasterplot files (with, for example the viewer 'xv'):

        % cd ../plots
        % xv compare* &

comparative rasterplot
Fig. 18:  Case 7.9 fluid 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 - water, oil.

water-oil interface
Fig. 19:  Geomview plot of water-oil interface surface for a randomly chosed subvolume of size [0,86]x[0,62]x[[21,153] as produced from executing Case 7.8sub.in



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


Authors: Masa Prodanovic, Brent Lindquist