1.4937540

Identifying layers in random multiphase structures
Kevin Mader, and Marco Stampanoni
Citation: AIP Conference Proceedings 1696, 020046 (2016);
View online: https://doi.org/10.1063/1.4937540
View Table of Contents: http://aip.scitation.org/toc/apc/1696/1
Published by the American Institute of Physics
Identifying Layers in Random Multiphase Structures
Kevin Mader∗ and Marco Stampanoni†
∗
†
4Quant Ltd., Switzerland & Institute for Biomedical Engineering at University and ETH Zurich, Switzerland
Institute for Biomedical Engineering at University and ETH Zurich, Switzerland & Swiss Light Source at Paul
Scherrer Institut, Villigen, Switzerland
Abstract.
X-Ray microscopic methods, benefiting from the large penetration depth of X-rays in many materials, enable 3D investigation of a wide variety of samples. This allows for a wide variety of physical, chemical, and biological structures to be seen and
explored, in some cases even in real time. Such measurements have lead to insights into paleontology, vulcanology, genetics,
and material science. The ability to see and visualize complex systems can provide otherwise unobtainable information on
structure, interactions, mechanical behavior, and evolution. The field has, however, led to a massive amount of new, heterogenous, difficult to process data. We present a general, model-free approach for characterizing multiphase 3D systems and show
how the method can be applied to experimental X-ray microscopy data to better understand and quantify layer structure in
two typical systems: investigation of layered fibers and clay samples.
Keywords: x-ray imaging; x-ray microscopy; k-means; morphology; 3D analysis; 3D structure
PACS: 68.37.Yz, 82.80.Ej
INTRODUCTION
The primary challenge of for the analysis of 3D images is reduction. A standard measurement at the TOMCAT
Beamline at the Swiss Light Source results in a 2560x2560x2160 array of absorption values for each time point
measured (48 GB of data, 3 times that for DPC phase measurements). Ultimately for quantitative analysis, comparison,
and statistical analysis, these data need to be collapsed into a few metrics which describe different aspects of the
system. This presents a great challenge in properly extracting the useful information from the system. Because each
sample type is different, there are a huge number of sample-specific, model-based procedures for analyzing the data.
This presents a challenge in not only implementation, but also application, comparsion, and understanding of a wide
variety of non-standard metrics with many caveats, levels of parameter tuning, and resolution sensitivity. Some of
the methods commonly used are: for lungs and vascular networks a branching tree structure [1, 2, 3], for trabecular
bone a plate and rod structure [4], for round pores and cells an ellipsoid shape analysis. Furthermore there are many
structures which do not yet have a suitable model or are generally difficult to enclose in such a framework such as
mixtures like geological stone, ice cream, bone scaffolds, and even cell agglomerations. In particular structures with
significant heterogeneity and randomness in volume, shape, and spatial distribution are very difficult to characterize
using these standard approaches.
Interest in understanding the structure of heterogeneous materials has existed long before 3D imaging methods and
computational power have been able to measure them well. These materials are regularly simulated and modeled
using analytical and computational approaches. In particular a family of functions known as N-Point correlation
functions have been the basis for much of this work [5, 6, 7, 8]. The analysis results in a significant reduction
amount of information by collapsing the expansive, dense 3D data to one or more simple correlation functions. This
reduction makes patterns, anisotropies, and distributions in the samples much more evident and easier to characterize
and ultimately quantify. Additionally this transformation allows for the deep understanding of the N-point correlation
functions as applied to a number of different real and synthetic systems to be applied to 3D microscopy measurements.
From these functions it is much easier to visualize and extract single scalar quantities since they are single curves
rather than full 3D datasets. Furthermore, these quantities can further be used to compare thousands of samples, or
evolution in a time-elapsed and separate sample experiments and can be verified and validated through robustness and
reproducibility studies. Finally since the full 3D information is still available these findings can be mapped back to the
3D structure where they provide a guide for what type of changes are significant.
In this manuscript, we first demonstrate how the algorithm can be efficiently implemented and applied to a wide
variety of multi-phase arbitrarily structured 3D data. We then apply the technique to several synthetic to show its
XRM 2014
AIP Conf. Proc. 1696, 020046-1–020046-13; doi: 10.1063/1.4937540
© 2016 AIP Publishing LLC 978-0-7354-1343-6/$30.00
020046-1
relation to existing metrics. Finally we apply the technique to real systems and show how understanding of the physical,
chemical, and biological processes driving the structural evolution can be better understood through this tool.
MODEL
The basic theoretical model has been summarized from [5]. We follow a shorter version of the derivation as demonstrated in [6] with the primary definitions and derivations available in the supplemental materials section. The model
developed can then be used to characterize a number of different aspects of a given material. The radial distribution
function is the most simple of these functions. It involves taking the 2 point correlation function for the 3D case.
S2A (R) = S2A (x, x + R)
This function can be rewritten in terms of a conditional probability, which makes later analysis easier to interpret and
implement. Effectively the radial distribution function is the probability that point (x + R) is part of V̄A given point (x)
belongs to V̄A
S2A (R) = P(I A (x + R)|I A (x))
Since we originally defined a mask phase which is the region containing the sample, we can make a correction to this
function to account for the fact that we are only interested in the distribution of phase A inside the mask / sample
volume. Thus adding the additional condition to the probability that the point (x + R) is part of V̄mask . Particularly
important around sample borders, the condition plays an important role as otherwise these points would contribute
non-isotropically to the overall distribution.
P(I A (x + R)|I A (x)I mask (x + R))
This correction introduces the potential singularity at positions where S2mask (R) = 0, we therefore assume the B material
is evenly enough distributed (which was the case in all experimental data) to have very few of these cases and when it
is we simply leave the point as undefined.
MULTIPHASIC SYSTEMS
The set of correlation functions are very well suited for analyzing multiphase systems because many of the standard
analysis approaches are only designed for 1-2 phases and lack the ability to quantify more integrated complicated
systems.
Composite Fiber Analysis
In many different types of manufacturing processes quality control is a critical issue. In many different fields
thickness, separation of components, and variability within and between components determines the success of a
product. We take an example of composite fibers (fig. ) to illustrate how phase-contrast tomography and the 2 point
correlation function can provide quantitative information on the thickness and separation of layers. The system consists
of carbon fibers (rods) covered in two layers. The first is a low-density SiC compound and the second a higher
density SiC compound. Finally there are void regions filled with nothing. Phase contrast methods provide strong
contrast between the carbon, low, and high density SiC layers. We are most interested in the coatings of the wires and
specifically how thick each layer is and if there are gaps between the layers. To solve this we employ the surface radial
distribution function using the points on the surface of the carbon rods as the I A (x). We furthermore define the low
density as I L , the high density as I H and the void as IV and characterize
PA→L (R) = I L (x + R)I A (x)
PA→H (R) = I H (x + R)I A (x)
PA→V (R) = IV (x + R)I A (x)
020046-2
C-Core
Higher Density SiC layer
Low Density SiC layer
100µm
1: On the left is a typical phase-contrast scan (performed by S.C. Irvine from Paul Scherrer Institute, Villigen,
Switzerland, sample provided by S. Stock from Institute for BioNanoTechnology in Medicine, Northwestern, USA).
Due to the strong contrast in the sample, basic segmentation using thresholds (right image) are able to separate the
image into the 3 distinct phases (carbon and low and high density SiC).
We then use these equations and their final values to assess the structure of the fibers (fig. ). We determine the
thickness by looking at the range of r values which each phase occupies and we assess variability by examining the
range of the values shown. We find that the first layer (low density) has a high degree of variability but surprisingly
the second (high density) appears to have very low variability since the values are additive (the second layer is on top
of the first layer). If the second layer contributed to the variability, the range of the values for the second layer surface
radial distance would have to exceed the range of the first layer and they appear to be very nearly the same. Finally
the void fraction shows a small peak near 0 indicating that near the surface of the rod or fiber there is in around 10%
of the cases (by surface area) a gap between the core and the first layer.
Distribution of clay around radiomarkers
The long-term storage and management of radioactive waste is a very important issue and understanding how these
dense materials integrate themselves into the ground and spread on a microscopic scale is important for understanding
the large trends and estimating the time-scales of leakage and potential contamination for the waste stored. While
X-ray based tomographic techniques can very well distinguish between the different phases in a standard sample
environment consisting of dense materials, clay, and calcite, the task of determining the relationships between these
phases is very difficult and eludes visual inspection (fig. ). The first two issues we examine are the propensity for the
dense materials to embed themselves in clay or calcite. The second is the anisotropy of the samples and what sort of
role compaction forces and gravity might play on the distributions.
Both of these questions can be very well addressed using the surface radial distribution function using the points on
the surface of the dense materials as the I A (x). We furthermore define the clay as IY , the calcite as I T .
020046-3
Phase Probability
/ Composition
Probability of Material / Composition
R
Low density SiC
High
density SiC
Inner
Outer ring
Void
/
Void Gap
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0
100
25µm
200
300
400
500
Distance from rod Material Surface (microns)
50µm
Core Surface Distance
Up to 10% surface
detachment between
core and inner layer
600
75µm
R
(µm)
Maximum inner layer thickness
Minimum inner layer thickness
2: An ROI single rod is extracted on the left for visualization purposes and the surface radial distribution function is
used to characterize the low and high density SiC layers in reference to the surface of the rod. The graph on the right
shows the 1D representation of this function and the values when looking at each phase (low - blue, high - red, and
void -green). The low shows a rather variable thickness from the surface of between 30μm and 60μm, while the high
shows a consistent thickness (the distribution is not broadened) of around 30μ (the distance between the peak of the
red curve and the peak of the blue curve). Finally the void distribution shows a value approaching 0.1 near the surface
indicating that around 10% by area of the surface is not attached to the first layer.
PA→Y (R) = IY (x + R)I A (x)
PA→T (R) = I T (x + R)I A (x)
(1)
(2)
Since the material consists of only these three phases (meaning their sum is total volume) and the volume fraction
of dense materials is low compared to the other two materials, we can interpret the probability values in a much more
concrete manner as volume fractions (fig. ). The result of the analysis shows that the dense materials preferentially
embed in calcite and the distribution of calcite and clay around the materials is anisotropic indicating another factor
at play to arrange the phases and consistently bias orientations. Combining this knowledge with the in-silico models,
we see that both the absolute difference and the anisotropy are indicative of a layered distribution. Specifically the
orientation of the layers is perpendicular to the direction of anisotropy. Given the rather small visibility or peak to
trough magnitude the variation in the sample is high and similar to the models tested with higher variation. Visual
inspection hints at this layered behavior as well, but would be very difficult to characterize, quantify and compare.
020046-4
Opalinus Clay
Calcite
Dense Minerals
Pyrit FeS2, zinc blende ZnS and siderite FeCO3
3: A ground sample was measured and segmented by A. Jakob of the Laboratory for Waste Management and D.
Grolimund from the microXAS Bealine at the Paul Scherrer Institute. Here the segmented image (using thresholds)
is shown of the calcite, clay, and dense materials samples. The different phases are very complex and discerning their
interconnectivity and spatial relationships from an image alone is quite difficult.
CONCLUSION
In this text we have shown how the radial distribution function can be used to quantify layered structures in multiphasic
random samples. We examine two very different samples, the first with a clear radial structure and well defined phases.
The second with no apparent structure. The method is equally applicable to both images and reduces the dimensionality
of the 3D imaging data substantially.
We then perform identical analyses on a synthetically generated set of samples with known layered structure which
we use as a reference for interpreting differences and variations in anisotropy.
The application of the Radial Distribution Function has long been used in theoretical material science applications
and for selecting representative volumes (a volume is representative when the RDF matches the RDF of the entire
volume within a given tolerance). With this text, we expand the utility of the tool for easily analyzing samples in a
model-free manner. Furthermore it is substantially less sensitive to noise or parameters than other metrics, particularly
surface or meshing-based analyses of structures. While we have demonstrated the utility on layered structures, there is
a large variety of other samples and materials which can be analyzed following the in-silico coupled analysis presented
in this manuscript.
020046-5
Volume Fraction
Volume Fraction
Asymmetry in distribution
- full 3D analysis
required
Distance from Dense Interface
Radial Angle
4: The graphs show the characterization of the complex 3 phase material and in particular the relationships between
the dense materials and the clay and calcite phases. The left graph shows the volume fraction of the two materials as
a function of the distance away from the surface of the dense minerals. The second analysis, a determination of the
degree of anisotropy in the system, shows the volume fraction as a function of radial angle. The system shows distinct
anisotropic behavior because the calcite abundance over clay is only present at certain radial angles and at -60 degrees
(left arrow) the two phases are equally present. The systems returns to the original difference at the other angles (right
arrow).
Acknowledgements
We would like to acknowledge S.C. Irvine (PSI), J. Andreas (PSI), and D. Grolimund (PSI) for the measurements
and analysis and the Swiss Light Source for providing beamtime and support for these experiments.
REFERENCES
1.
2.
3.
4.
5.
6.
7.
8.
9.
S. Heinzer, G. Kuhn, T. Krucker, E. Meyer, A. Ulmann-Schuler, M. Stampanoni, M. Gassmann, H. H. Marti, R. Müller,
and J. Vogel, NeuroImage 39, 1549–58 (2008), ISSN 1053-8119, URL http://www.ncbi.nlm.nih.gov/pubmed/
18077185.
S. Heinzer, T. Krucker, M. Stampanoni, R. Abela, E. P. Meyer, A. Schuler, P. Schneider, and R. Müller, NeuroImage 32,
626–36 (2006), ISSN 1053-8119, URL http://www.ncbi.nlm.nih.gov/pubmed/16697665.
S. I. Mund, M. Stampanoni, and J. C. Schittny, Developmental dynamics : an official publication of the American Association
of Anatomists 237, 2108–16 (2008), ISSN 1058-8388, URL http://www.ncbi.nlm.nih.gov/pubmed/18651668.
M. Stauber, and R. Müller, Bone 38, 475–84 (2006), ISSN 8756-3282, URL http://www.ncbi.nlm.nih.gov/
pubmed/16338187.
S. Torquato, Random Heterogeneous Materials, Springer, 2001, URL http://www.amazon.com/
Random-Heterogeneous-Materials-Salvatore-Torquato/dp/0387951679.
Y. Jiao, F. Stillinger, and S. Torquato, Physical Review E 76 (2007), ISSN 1539-3755, URL http://pre.aps.org/
abstract/PRE/v76/i3/e031110.
C. Yeong, and S. Torquato, Physical Review E 58, 224–233 (1998), ISSN 1063-651X, URL http://adsabs.harvard.
edu/abs/1998PhRvE..58..224Y.
S. C. Blair, P. A. Berge, and J. G. Berryman, Journal of Geophysical Research 101, 20359–20375 (1996), ISSN 0148-0227,
URL http://www.agu.org/pubs/crossref/1996/96JB00879.shtml.
K. Mader, High-throughput, synchrotron based tomographic microscopy tools for the quantitative characterization of complex
structures (2013), URL http://e-collection.library.ethz.ch/view/eth:6915.
020046-6
10. K. S. Mader, P. Schneider, R. Müller, and M. Stampanoni, Bone 57, 142–154 (2013), ISSN 1873-2763, URL http:
//dx.doi.org/10.1016/j.bone.2013.06.026http://www.ncbi.nlm.nih.gov/pubmed/23871748.
11. K. Mader, R. Mokso, and C. Raufaste, Colloids and Surfaces A: . . . 415, 230–238 (2012), ISSN 09277757,
URL http://dx.doi.org/10.1016/j.colsurfa.2012.09.007http://www.sciencedirect.com/
science/article/pii/S0927775712006188.
SUPPLEMENTARY MATERIALS
Model
The more complete derivation of the model summarized below is shown in the Thesis of K. Mader [9]. For the
first case the system examined a two-phase system (rock and air, bone and marrow, liquid and bubbles, etc.). We shall
call the firstphase A and the second B. We define the mask (V̄mask ) as being the region which includes both phases
(V̄mask = V̄A V̄B ) and therefore every point in V̄mask must be of either phase A or phase B.
Indicator Functions
We define the regions belonging to phase A and B as V̄A
indicator functions as
1,
A
I (x) =
0,
1,
I B (x) =
0,
1,
I mask (x) =
0,
and V̄B respectively. From these regions, we define the
x ∈ V̄A
x ∈ V̄A
x ∈ V̄B
x ∈ V̄B
x ∈ V̄mask
x ∈ V̄mask
For the rest of the section we will focus on phase A, knowing that the functions created are equally applicable to
phase B. ### N-Point Correlation With the indicator function defined, it is now possible to create N-point correlation
functions
N
SNA (x1 , x2 , · · · , xN ) = ∏ I A (xi )
i=1
If the distribution of phases is statistically homogeneous, we can simplify the correlation function by examining
relative displacements rather than distinct points (Ri = xi − x1 )
N
SNA (R2 , R3 , · · · , RN ) = ∏ I A (x + Ri )
i=1
For distributions, which are locally homogeneous, but globally heterogeneous, it is possible to make the approximation
for each of these local regions and calculate the function for that region resulting in many correlation functions. This
approach is particularly well suited for characterizing the variations in structure within a sample since these correlation
functions can then be directly compared.
Simulated Data
Interpretation of radial distribution functions in measured samples with high levels of variability is challenging. We
approach this problem by synthetically creating (in-silico) layered samples with known thicknesses and orientations,
and then comparing the RDF to the known values. To the synthetic samples, various degrees of random fluctuations
can be added to simulate a progressively more and more disordered sample. Several example samples are shown below
and their corresponding RDF functions. For simplicity rather than having generic phases A and B, we used the clay
and calcite names from the examples so the results are less abstract.
020046-7
Input Images
•
A well-ordered vertically layered material with 3 phases. The red point represents the first phase and the green
and purple the other two.
020046-8
•
A less ordered vertically layered material with 3 phases. The green and purple phases are much more strongly
intermixed in this example making the structure more difficult to initially identify.
020046-9
•
A radially ordered layer system. The layers from radial shells instead of vertical slabs, but once the coordinate
system has been changed the analysis is identical
020046-10
Calculated Radial Distribution Functions
10%
20
30%
50%
Z − XY Plane Angle
10
0
−10
−20
70%
20
val
1.00
0.75
0.50
0.25
0.00
90%
10
0
−10
−20
−20 −10
•
0
10
20−20 −10
0
10
XY Angle
20
The graph shows the probility of the clay phase as a function of radial angle for increasing levels of random
fluctuatings added to the image (as shown above). The image remains fairly robust until 70% of the pixels are
changed at which point the image itself is no longer recognizeable as a layered structure.
5.85
13.2
15.8
17.7
19.2
45
val
1.00
0.75
0.50
0.25
0.00
90
•
30
Z − XY Plane Angle
0
150
100
50
0
150
100
50
0
150
100
50
0
150
100
50
0
−100 0 100
−100 0 100
−100 0 100
XY Angle
−100 0 100
−100 0 100
The graph shows the probability of the clay phase against the radial angles (XY and in Z-XY) in different panels.
The vertical panels show the orietnation of the structure itself while the horizontal panels show the radial distance.
020046-11
20
−10
0
10
Y Displacement
Calcite
10
0
−10
−20
20
10
Clay
0
Probability
0.00
0.25
0.50
0.75
1.00
Phase
Calcite
Clay
−10
−20
−20 −10
•
0
10
20
−20 −10
0
10
X Displacement
20
−20 −10
0
10
20
The radial slab has a radial distribution function as expected with the probabilities of each phase mirroring the
shape of the sample
Interpretation
The radial distribution function of the progessively noisier images show the robustness of the analysis against noise
and high degrees of variability which is common in experimental data. The next chart shows how the information on
both the orientation and thickness of the layers can be infered from the 3D information in RDF by examining the
anisotropy of the distribution at various radial distances. The degree of fluctuation indicates the angle between 0 (no
flucatuation at all) and 90 (a dot pattern). The angles in present themselves as progressively higher amplitude sinus
curves. The thickness can be infered from the relative differences between the distribution at different radial positions.
RDF Implementation
The implementation of the series of correlation functions was done as a plugin as part of the Spark Image Layer
framework ([10, 11]). The plugin takes a boolean or label image as the first input and has an option second input which
can be a grayscale image. The plugin then requires the size of the discretized box (B = [−r, −r, −r] → [r, r, r]) and
creates two temporary arrays. The first is the cumulative array (Ĉ), which stores the total count or score for that position.
The second is the incidence array (N̂), which stores how often that value was visited without being incremented. The
Ĉ
final value is then the cumulative array divided element-wise by the incidence array ( N̂i j ). The entire code base can
ij
then be easily parallelized without lock issues by storing the immutable input array(s) in a shared memory space and
creating and storing the mutable Ĉ and N̂ arrays separately for each thread or task. At the end the cumulative arrays
can be summed and the incidence arrays summed before their division.
The arrays are filled in one of two ways. The stand method involves scanning across the entire I A (x) == 1 space and
scanning through (x + R)∀R ∈ B. As it is extremely computationally expensive to iterate over the entire data set and
we have already assumed spatial statistical homogeneity, we use a Monte Carlo technique to get a meaningful sample
of the distributions. We thus randomly sample the function at between 50,000-100,000 randomly selected positions
and build up the Ĉ and N̂ this way. The iteration is run until the all voxels in the function change between iterations
less than tol = 1e − 6 relative to the mean value (a parameter which can be tuned depending on the need for accuracy).
This means values at or below 1e − 6 of the mean are not reliable. The exhaustive approach requires weeks on a 2000
020046-12
x 2000 x 2000 data set with a 100 x 100 x 100 window, while the Monte Carlo runs approximately 20-80,000x faster
depending on the needed sensitivity.
020046-13