A Software Package for
Automated Analysis of Fiber Structure
in 3-D Computed Microtomography Images
| Software Authors: | Prof. W.B. Lindquist, Dr. H. Yang, M. Oh, Dr. S.-M. Lee, Dr. W. Oh, Dr. A.B. Venkatarangan, Dr. H. Shin, M. Prodanovic |
| Documentation Author: | W.B. Lindquist |
|
| Fig. 1: Three-dimensional network (medial axis) structure of an unwoven fiber mat obtained from analysis of a computed tomographic image by 3DMA-Fiber. Green fibers are thicker than red. |
Table of Contents
A) Introduction
B) Segmentation
C) Construction and modification of the medial axis
D) Specific surface area computation
E) Fiber reconstruction
F) Fiber statistics
G) Void space tortuosity
H) References
A) Introduction
( return to table of contents) ( next section)
3DMA-Fiber is a software package for analyzing the fiber network in three- (and two-) dimensional X-ray computed microtomographic (XCMT) images of fiber mats.
The figures accompanying this html file are best viewed with the html viewer window set to full screen.
The image analysis algorithms for the fiber network consist of six general steps: image segmentation, extraction and modification of the medial axis of the fiber network, fiber surface construction via marching cubes, fiber reconstruction using the medial axis as a search path, fiber statistics, and computation of void space tortuosity. These steps are discussed below, largely by illustration. They are also discussed in the 3DMA-Fiber primer.
B) Segmentation
( previous section) ( return to table of contents) ( next section)
The greyscale intensity of each voxel is typically an integer value (from 0-255 for 8-bit images). Fig. 2 shows views of a single slice from a 3D image stack of an unwoven fiber mat. 3DMA-Fiber provides three options for image segmentation: simple (global-threshold) segmentation, Mardia-Hainsworth (MH) [1988] and indicator kriging (IK) [Oh and Lindquist, 1999]. For simple segmentation, the user must enter a single threshold: voxels having intensity above this threshold are assigned to the fiber phase; remaining voxels are assigned to the void phase. Both MH and IK methods are locally adaptive and are similar in that both utilize (an estimate of) the two-point correlation function of the image. MH is an iterative method, requiring a single threshold value which it modifies each iteration. IK is a "single pass" method requiring that a subpopulation of voxels of each phase be positively identified a priori. This can be done by manually establishing a window of intensity values delimited by two thresholds. Intensity values above the upper threshold are assumed fiber, those below the lower threshold are assumed void space. The classification of remaining voxels is estimated by indicator kriging which utilizes an estimate of a correlation function incorporating local spatial information. These three methods are compared in Oh and Lindquist on a test images involving berea sandstone. The simple thresholding and IK methods are compared in Fig. 3.
| ||
| Fig. 2: Views of a single slice from a tomographic data stack of the image of an unwoven fiber (left) without and (right) with the use of a fiducial polygon mask to remove data of no interest. |
| ||
| Fig. 3: View of the same single slice as in Fig. 2 under (left) simple thresholding and (right) indicator kriging based tresholding. Note the susceptibility of simple thresolding to "speckled" noise. |
C) Construction and modification of the medial axis
( previous section) ( return to table of contents) ( next section)
The medial axis [Sirjani and Cross, 1991] of a digitized object is a 26-connected (or 6-connected) centrally-located skeleton which preserves the original topology and geometry of the object. If their are no cavities in the object, the medial axis is a network of voxel paths. Erosion-based algorithms are used to construct the medial axis; 3DMA-Fiber employs the algorithm developed by Lee et al. [1994]. The medial axis is sensitive to surface features, including noise; as a result, the medial axis may contain spurious paths that are not significant descriptors of the overall geometry of the object. Spurious paths can be trimmed; 3DMA-Fiber provides a menu of options to control this trimming. Fig. 4 illustrates trimming on a medial axis obtained from a two-dimensional data set consisting of a single slice. Notice how the trimming is used to remove the medial axis associated with short segments.
| ||
| Fig. 4: View of the medial axis obtained from a two dimensional data set consisting of a single slice (left) before and (right) after trimming. Trimming is used to remove the medial axis associated with short length fiber segments. |
Fig. 1 shows the entire medial axis extracted the 3D unwoven fiber image after medial axis trimming. More detailed information on trimming can be found in the 3DMA-Fiber primer and the historical documents: the Online Manual to 3DMA, and the 3DMA General Users Manual, describing the underlying code upon which 3DMA-Fiber is constructed.
D) Specific surface area computation
( previous section) ( return to table of contents) ( next section)
We utilize the marching cubes algorithm [Bloomenthal, 1988; Lorensen and Cline, 1987] to determine the surface of each fiber, from which the specific surface area (surface area per unit volume) is determinted. Fig. 5 shows the resulting fiber surface computed in a small subvolume of the unwoven fiber mesh.
|
| Fig. 5: The fiber surface computed using the marching cubes algorithm for a small subvolume of the 3D image of an unwoven fiber mesh. |
E) Fiber reconstruction
( previous section) ( return to table of contents) ( next section)
Fiber reconstruction consists of connecting medial axis paths into long fibers. If we refer to a "branch cluster" as a set of (1 or more) voxels at which paths cross each other, the main problem to fiber reconstruction would appear to be pairing paths that "lie 180 degrees from each other" (allowing for the possbility of some bending) on each side of the branch cluster. However, as indicated in Fig. 6, the medial axis local to the crossing (bonding) site of two fibers may consist of two close coordination number 3 clusters rather than a single cluster of coordination number 4. (This occurs because the bond between two fiber occurs at their surface, whereas the medial axis traces fiber midlines.) Thus, close cluster must be merged before tracing across clusters can be completed.
|
| Fig. 6: A medial axis resolving itself into two coordination number three clusters rather than a single coordination number 4 cluster at the bonding site of two fibers. |
When two fibers bond at an acute angle of incidence, the medial axis of the digitized sample can form a complex group of branch_clusters at the bonding site as illustrated in Fig. 7 (left). Cluster merging does not simpligy the complex bonding structure (which we refer to, somewhat incorrectly, as a "surface remnant"). A "surface remnant reduction" algorithm is provided to simplify such structures as indicated in Fig. 7 (right).
| ||
| Fig. 7: (Left) View of the complex merged branch cluster formed in the medial axis of the digitized image when fibers bond at acute angle of incidence. (Right) The complex branch cluster ("surface remnant") after reduction. |
The two publications Three-dimensional image analysis of fibrous material and A geometric analysis of 3D fiber networks from high resolution images discuss the fiber reconstruction in more detail.
F) Fiber statistics
( previous section) ( return to table of contents) ( next section)
As fiber properties are statistical, we characterize the spatial variation of any characterizing parameter by plotting its distribution of values determined from th3e image. Distributions are produced for: free fiber segment length, fiber length, bonds per fiber, bonds per fiber length, fiber crossing angles and fiber curl. Distributions for four of these parameters, determined for the unwoven fiber mat are shown in Fig. 8.
| ||||||
| Fig. 8: Distributions for (top left) fiber length; (top right) bonds per fiber; (bottom left) fiber crossing angle; and (bottom right) fiber curl determined from the 3DMA-Fiber analysis a 3D image of an unwover fiber mat. |
G) Void space tortuosity
( previous section) ( return to table of contents) ( next section)
3DMA-Fiber has the ability to compute the tortuosity through the void space. This is illustrated with test data in the in the 3DMA-Fiber primer. Here is a quick description of the algorithm used to compute, and results for the computation of, geometrical tortuosities of the void space.
With a trimmed medial axis computed for the void space inside a rectangular
volume and output in MA voxel list format, tortuosity distributions for paths
between any two faces of the rectangular volume can be computed.
|
| 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 (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_34 be the Euclidean distance between faces 3 and 4. The geometrical tortuosity t_ij of the path is defined as t_ij = l_ij/d_34. 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 ways. The first is the geometrical tortuosity distribution for all entrance/exit voxel pairs:
{t_ij | i is an entrance (face 3) MA voxel, j is an exit (face 4) MA voxel},
The second is the geometrical tortuosity distribution for the shortest path connecting each entrance voxel i to the exit face:
{t_i = min_j (t_ij) | j is an exit face MA voxel},
where i ranges over all MA voxels on the extrance face.
Fig. 10 shows the shortest paths in each direction obtained from a small
subvolume of the unwoven fiber data set for which these two tortuosity
distributions are computed. Fig. 11 shows the tortuosity distributions obtained.
| ||||||
| Fig. 10: Geomview plots of the medial axis paths from: (top) faces 3 -> 4 (x-direction); (middle) faces 5 -> 6 (y-direction); (bottom) faces 1 -> 2 (z-direction). For each direction, the left-hand plot shows the shortest paths between every inlet-outlet voxel pair while the right-hand plot shows only the shortest paths across the image from each inlet voxel. Note that the y-direction points toward the viewer. |
 |
| Fig. 11: Distributions of void space tortuosities in the (left) x- (middle) y- and (right) z-directions from the plots of Fig. 10. Top histograms correspond to the plots on the left of Fig. 10, bottom histograms to the plots on the right of Fig. 10. |
In an analogous manner, starting with the segmented files v_c_seg/fbr1.zzz.gz, one can also compute paths through the medial axis of the fiber nextwork and obtain fiber pathway tortuosities. (The changes to the files case2.5v.in through case4.6zv.in needed to do this should be straightforward.) In Fig. 12, we show the tortuosity distributions we obtain for the fiber network for this data set. Note the higher tortuosity of fiber paths in the z-direction. In contrast to the void space, the x-y layering of the fibers makes it difficult to find "straight" paths through the fiber mesh in the z-direction.
 |
| Fig. 12: Distributions of fiber network tortuosities in the (left) x- (middle) y- and (right) z-directions from the v_c_seg/fbr1.zzz.gz data files. |
H) References
( previous section) ( return to table of contents)
J. Bloomenthal. Polygonization of implicit surfaces. IEEE Comput. Graph. App. 5 (1988) 341.
T. C. Lee, R. L. Kashyap, and C. N. Chu. Building skeleton models via {3-D} medial surface/axis thinning algorithms. CVGIP: Graph. Models Image Process. 56 (1994) 462-478.
W.B. Lindquist. 3DMA General Users Manual.
W.E. Lorensen and H.E. Cline. Marching cubes: a high resolution 3-D surface construction. ACM Comput. Graph. 21 (1987) 163.
K.V. Mardia and T.J. Hainsworth. A spatial thresholding method for image segmentation. IEEE Trans. Pattern Anal. Machine Intell. 10 (1988) 919-927.
W. Oh and W. B. Lindquist. Image thresholding by indicator kriging. IEEE Trans. Pattern Anal. Mach. Intell. 21 (1999) 590-602. ftp://ftp.ams.sunysb.edu/pub/papers/1998/susb98_02.ps.gz
A. Sirjani and G. R. Cross. On representation of a shape's skeleton. Pattern Recogn. Lett. 12 (1991) 149-154.
H. Yang and W.B. Lindquist. Three-dimensional image analysis of fibrous material
H. Yang and W.B. Lindquist. A geometric analysis of 3D fiber networks from high resolution images.