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... 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 tomog cases burn seg c_seg ma ma_t 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
results    graphical output on geometrical characterization
throats    throat and pore network files

- Use the following link to download a set of example tomographic input files into your /nfs/3dma_rock/Test/tomog/ directory.

- Unzip and extract the tomographic files.

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

The example contains 100 files labeled 'oil.XXX.gz' where XXX ranges from 001 to 100, containing 100 slices (the z-direction) of a 100x100x100 volume of a polyethylene synthetic core filled with oil. Each slice is an array of 10^4 numbers (unsigned chars, values ranging from 0 to 255) representing X-ray attenuation coefficients.




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
	special purpose data requiring pre-processing (0)
	tomographic data (1)
	segmented data (2)
	burn data (3)
	medial axis data (4)
	throat data (5)
	fluid data (6)
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.

( previous section)    ( return to table of contents)    ( next section)

List of input files   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.3.in   volume information on disconnected pores and grains
case2.4.in   "clean-up" operations on segmented data
case2.5.in   erosion computation of media axis (Lee,Hasyap, Chu [1994] algorithm)
case4.11.in   trimming operations on medial axis
case4.1ras.in   slice-by-slice rasterplots of medial axis
case4.1iv.in   volume plots of medial axis via Inventor
case4.13dij.in   computation of pore throats using Venkatarangan [2000] algorithm
case4.13ray.in   computation of pore throats using Shin [2002] algorithm
case5.1ras.in   slice-by-slice views of throats via rasterplots
case5.1iv.in   volume plots of throats via Inventor
case5.1iv_subvol.in   subvolume plots of throats via Inventor
case5.3.in   construction of pore-throat network
case5.2.in   output of pore-throat statistics




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

Case1.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/cases/ directory.

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

	- This will produce the output file case1.1.out and rasterfiles
	  oil.001.ras.gz to oil.020.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 oil.*.ras.gz
Fig. 1 shows the resulting plot in oil.002.ras.gz.

Fig. 1: The rasterfile image (magnified by a factor of 2) of the test
example oil.002.ras.gz. The dark phase represents oil, and the light
part is polyethylene core.

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. Case1.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/cases/ directory.

	% cd /nfs/3dma_rock/Test/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/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 (the polyethylene matrix in this case) and the upper mode to the higher attenuating phase (hexadecane doped with bromide).



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 130 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/cases/ directory.

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

	- This will produce the output file case1.5smpl.out,
	  intermediate data files oil_smpl.XXX.gz in the /nfs/3dma_rock/Test/seg
	  directory, and rasterfiles oil_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 oil_smpl.XXX.gz data
	  files.

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


**So why print both the data files oil_smpl.XXX.gz and the rasterfiles oil_smpl.XXX.ras.gz which, after all, contain the same information. The reason is one of disc usage. The data files, which we refer to a segmentation 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 for test example simple segmentation results for slice 2
in the 3D image stack. Oil is black; polyethylene appears white.

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




Case1.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/cases/ directory.

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

	- This will produce the output file case1.5smpl.out,
	  intermediate data files oil.XXX.gz in the /nfs/3dma_rock/Test/seg
	  directory, and kriging summary raster files oil.XXX.ras.gz also in
	  the seg directory.

	  The raster files provide a visual synopsis (Fig. 4) of the kriging
	  decisions for each slice.

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

Fig. 4: Kriging summary image (x2 magnification) for the second slice in the test example.
(left) The tomographic data. (middle) 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.
(right) The final segmented result.

With the image segmented, all further processing steps utilize the segmented data files (stored in the seg/ directory. 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. Case2.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, there 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/cases/ directory.

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

	- This will produce the output file case2.3.out and a jgraph file
	  oil_dv.jgr in the /nfs/3dma_rock/Test/results directory. The
	  distribution of volume sizes (in numbers of voxels) ordered by volume
	  for both disconnected pore and polyethylene blobs is printed in the
	  case2.3.out file:
	  
	  disconnected pore** volumes
	  total volume 555336 (disconnected 0.000122428 connected 0.999878)
	  size    number
	  1       4
	  4       1
	  28      1
	  32      1
	  555268  1

	  disconnected grain** volumes
	  total volume 444664 (disconnected 0.000206888 connected 0.999793)
	  size    number
	  1       12
	  2       1
	  78      1
	  444572  1
	  
	  To delete all blobs but the largest of each phase, volume
	  thresholds (e.g. 40 for phase 0 (polyethylene) and 80 for phase 1
	  (hexadecane) can be established simply examining this printout.
	  The chosen thresholds can then be entered into case2.4.in
	  (In fact the threshold value of 50 is entered for both in the
	  example case2.4.in provided.)

	% vi case2.4.in
	% runcase linux case2.4

	- This stores new "cleaned" segmented files oil.XXX.gz and raster file
	plots oil.XXX.ras.gz of these in the /nfs/3dma_rock/Test/c_seg/ directory.

	% cd ../c_seg/
	% xv oil.*.ras.gz

**In the output of case2.3, the word "pore" refers to phase 0 voxels (which are in fact polyethylene), while the work "grain" refers to phase 1 (which is in fact hexadecane in the void space).

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

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




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

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

	- This will produce the output file case2.5.out,
	  intermediate burn number data files oil.XXX.gz in the burn/ directory
	  and intermediate medial axis data files oil.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.11: 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.11.in 
	  into the /nfs/3dma_rock/Test/cases/ directory.

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

	- 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 formate in
	  the ma_t/ directory.

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




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

Medial axis can be view either slicewise (rasterfile format), or in 3D (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.1iv.in 
	  into the /nfs/3dma_rock/Test/cases/ directory.

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

	- Case4.1ras produces raster plots oil.XXX.ras.gz in the ma_t/
	  directory. An example slice plot is shown in Fig. 5. Case4.1iv
	  produces an inventor graphic plot oil.iv.gz in the ma_t/ directory.
	  The uncompressed oil.iv file will only be viewable with an inventor
	  graphics interpreter package. A view of this 3D plot is shown in
	  Fig. 6.

Fig. 5: A rasterplot (2x magnification) of the (3D) medial
axis voxels visible in slice 20 of the example data.
 
Fig. 6: An inventor graphics view of the 3D medial axis of the
example data set. Grains are not plotted.



Case 4.13: 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 'ray-tracing' algorithm from the original motivation behind the algorithm's development.

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 ray 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.13dij.in and case4.13ray.in
	  into the /nfs/3dma_rock/Test/cases/ directory.

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

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

	% runcase linux case4.13ray

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

The throat file** contains complete information on each throat including information on it's location on the medial axis, perimeter, 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:  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.

Test example:
	- download the input files case5.1ras.in , case5.1iv.in  and
	  case5.1iv_subvol.in  into the /nfs/3dma_rock/Test/cases/ directory.

	% cd /nfs/3dma_rock/Test/cases/
	% runcase linux case5.1ras
	% runcase linux case5.1iv
	% runcase linux case5.1iv_subvol

	- Case4.1ras produces raster plots oil.XXX.ras.gz in the throat/
	  directory. Case5.1iv produces an inventor graphic plot oil.iv.gz
	  in the throat/ directory.  The uncompressed oil.iv file will only be
	  viewable with an inventor graphics interpreter package. A view of this
	  3D plot is shown in Fig. 7. Case5.1iv_subvol produces an inventor
	  graphic plot oil_p2.iv.gz in the throat/ directory. It shows a few
	  throats only in a selected subvolume. A view of this region is shown
	  in Fig. 8.

Fig. 7: An inventor graphics view of the throat regions (purple) and medial
axis (rainbox colored) in the example data set. The throat region is defined
as the set of voxels through which the computed throat surface passes.
 
Fig. 8: Three inventor views of the throat regions (purple) and medial axis (rainbox colored)
in a subvolume [36,48]x[50,67]x[59,78] of the example data set. The view is of a coordination
number 5 pore (coordination number 5 branch cluster) with a throat crossing each medial axis path
leaving the pore.




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. Case5.3 organizesthe pores and throats into a network with explicitly specified neighboring relationships.

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

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

	- This will produce the output file case5.3s.out, and intermediate data
	  files oil.th2np.gz and oil.np2th.gz in the throat directory.
The two files cross reference the network: oil.th2np cross references the throats to the pores, will oil.np2th cross references the pores each other and to the throats.



Case 5.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 case5.2.in 
	  into the /nfs/3dma_rock/Test/cases/ directory.

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

	- This will produce the output file case5.2.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** similar to that of np_thr_net.ps is presented in Fig. 9.


**The example data set used to produce np_thr_net.jgr contains too few pores to produce meaningful statistics. The data in Fig. 9 is a plot of from a larger sample of polyethylene.

Fig. 9: Distributions resulting from running case5.2 on a larger polyethylene sample.



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

Authors: Masa Prodanovic, Brent Lindquist