|
|
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. |
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.
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.Nis equivalent to typing the above stripcomm command.
| 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 |
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).
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
|
|
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.
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.
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
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.
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.
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/
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.
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. |
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.
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.
|
| Fig. 9: Distributions resulting from running case5.2 on a larger polyethylene sample. |