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, beneﬁting 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 ﬁeld has, however, led to a massive amount of new, heterogenous, difﬁcult 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 ﬁbers 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-speciﬁc, 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 difﬁcult to enclose in such a framework such as mixtures like geological stone, ice cream, bone scaffolds, and even cell agglomerations. In particular structures with signiﬁcant heterogeneity and randomness in volume, shape, and spatial distribution are very difﬁcult 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 signiﬁcant 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 veriﬁed and validated through robustness and reproducibility studies. Finally since the full 3D information is still available these ﬁndings can be mapped back to the 3D structure where they provide a guide for what type of changes are signiﬁcant. In this manuscript, we ﬁrst demonstrate how the algorithm can be efﬁciently 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 deﬁnitions 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 deﬁned 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 undeﬁned. 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 ﬁelds thickness, separation of components, and variability within and between components determines the success of a product. We take an example of composite ﬁbers (ﬁg. ) 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 ﬁbers (rods) covered in two layers. The ﬁrst is a low-density SiC compound and the second a higher density SiC compound. Finally there are void regions ﬁlled 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 speciﬁcally 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 deﬁne 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 ﬁnal values to assess the structure of the ﬁbers (ﬁg. ). 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 ﬁnd that the ﬁrst 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 ﬁrst 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 ﬁrst 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 ﬁber there is in around 10% of the cases (by surface area) a gap between the core and the ﬁrst 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 difﬁcult and eludes visual inspection (ﬁg. ). The ﬁrst 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 deﬁne 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 ﬁrst 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 (ﬁg. ). 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. Speciﬁcally 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 difﬁcult 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 difﬁcult. 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 ﬁrst with a clear radial structure and well deﬁned 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 ofﬁcial 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 ﬁrst case the system examined a two-phase system (rock and air, bone and marrow, liquid and bubbles, etc.). We shall call the ﬁrstphase A and the second B. We deﬁne 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 deﬁne 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 deﬁne 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 deﬁned, 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 ﬂuctuations 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 ﬁrst 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 difﬁcult 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 ﬂuctuatings 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 ﬂuctuation indicates the angle between 0 (no ﬂucatuation 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 ﬁrst 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 ﬁrst 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 Ĉ ﬁnal 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 ﬁlled 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

1/--страниц