Integrated approach to the data processing of four-dimensional datasets from phase-contrast x-ray tomography

: Phase contrast X-ray tomography (PCT) enables the study of systems consisting of elements with similar atomic numbers. Processing datasets acquired using PCT is nontrivial because of the low-pass characteristics of the commonly used single-image phase retrieval algorithm. In this study, we introduce an image processing methodology that simultaneously utilizes both phase and attenuation components of an image obtained at a single detector distance. This novel method, combined with regularized Perona-Malik ﬁlter and bias-corrected fuzzy C-means algorithm, allows for automated segmentation of data acquired through four-dimensional PCT. Using this integrated approach, the three-dimensional coarsening morphol-ogy of an Aluminum-29.9wt% Silicon alloy can be analyzed. morphological characterization of alloy during coarsening” (2014), Manuscript in preparation.


Introduction
Phase contrast x-ray tomography (PCT) enables the study of weakly absorbing samples, as well as systems consisting of elements with similar atomic numbers. This is because the real part of the refractive index δ dominates over the imaginary part β in PCT experiments [1,2].
In propagation-based PCT, a phase map is commonly obtained by applying phase-retrieval algorithms to the projections collected by a similar setup to Fig. 1 [3,4]. Then, filtered back projection [5] is applied to these phase maps in order to recover the refractive index decrement during reconstruction. This two-step approach [6] of phase-retrieval followed by backprojection will hereafter be referred to as PAG [3]. On the other hand, the one-step approach of using the filtered back projection algorithm to reconstruct the images directly from the traditional absorption-based projection images, collected at the same sample-to-detector distance (R 2 ) will be referred to as FBP. Robust segmentation of the PCT reconstructions is crucial for quantitative analysis of 3D structures, e.g., surface area of interfaces and interfacial curvature measurements [7].
The growing size of data collected during the experiments, typically on the order of terabytes, renders manual segmentation impractical. Furthermore, such an approach does not offer a high degree of reproducibility, accuracy, or consistency. However, automated segmentation of PCT reconstructions is nontrivial for the following reasons: • During the phase-retrieval step in the PAG approach, the tomograms undergo a smoothing operation. In other words, the single-image phase-retrieval algorithms that are conventionally used show inherently low-pass characteristics [8,9], which, in our study, leads to diffuse interfaces in the PAG reconstructions. Smoothing of the interface makes it difficult to determine accurately the interfacial morphology.
• On the other hand, the FBP images are characterized by dark-bright fringes at the interfaces, giving rise to the so-called halo effect [10]. In the near-field or short propagation regime, only the first-order Fresnel fringes are visible in the reconstructions [1,2]. These images are very challenging to segment.
Conventional image processing techniques, such as histogram thresholding and k-means clustering, fail to provide reliable results; for example, the "halo effect", described above, leads to spurious edges in the binary image [10]. Nevertheless, a few attempts have been made to semi-automatically segment PCT datasets. Ref. [11] used a weak watershed transform assembly to increase segmentation robustness. In this method, however, the centroids of the particles need to be marked by a user prior to segmentation. Ref. [12] presented a learning classifier to guide a constrained statistical shape model to fit the data; such a method is overly deterministic Fig. 1. PCT experiment, where R 2 is the detector distance, θ is the projection angle and λ is the wavelength of propagating wave. The frame (x, y, z) is the reference frame, while (r 1 , r 2 ) lie in the plane of the imaging detector [4]. and does not allow for particle shapes outside of the classifier to be segmented. It is for this reason that the rotating kernel transform [13] is also ineffective in segmenting multiple particles of variable shape and size. To circumvent the problems in segmenting PCT images, Ref. [14] proposed combining FBP and PAG images in Fourier space in their study of lung aveoli. Two-dimensional composite images were produced that were then easily thresholded [14]. We propose an integrated image processing methodology that utilizes both phase and absorption contrast, derived from applying PAG and FBP algorithms separately to the raw PCT data, along with a suite of data processing methods to allow the automated segmentation of large datasets acquired through 4D tomography. Unlike the work of Ref. [14], the images are combined in real space, and due to the complexity of the microstructure and the sensitivity of our measurements, a more robust image processing procedure is necessary. The hybrid images feature improved contrast-to-noise ratio and spatial resolution, thereby enabling the automated segmentation of such images by non-linear diffusion filtering and fuzzy logic. The binary images are then combined to reveal the 3D microstructures, in this case for an Al-29.9wt%Si alloy coarsening in time. To our knowledge, this is the first time that 4D PCT has been used to study the 3D interfacial morphologies of an alloy consisting of elements with similar atomic numbers.

Experimental methodology
Al-Si rods of composition 29.9 wt% Si were prepared by Ames Laboratory [15], see Ref. [16] for experimental details. The 4D propagation-based phase contrast tomography experiment was conducted ex situ at the TOmographic Microscopy and Coherent rAdiology ExperimenTs (TOMCAT) beamline of the Swiss Light Source (Paul Scherrer Institut, Switzerland) [17]. The sample-to-detector distance was set to 110.0 mm and optimized for the sample volume and a monochromatic X-ray energy of 28 keV. Such a setup satisfied the near-field condition for PCT [3,18].
The alloy consisted of large, interconnected Si laths in a eutectic matrix, see Fig. 2. The sample was placed in a custom-made isothermal furnace, and the Si laths were allowed to coarsen at 590 • C, just above the eutectic temperature of 577 • C. After 10 minutes, the sample was taken out of the furnace and tomographic projections were collected at room temperature when the sample was fully solid. The sample was then reheated to above the eutectic temperature, and kept at 590 • C for a subsequent 10 minutes, continuing the coarsening process. This cycle was repeated for six time iterations.
In between each iteration, 1001 projections were collected over 180 • . Phase-retrieval and subsequent reconstruction of the images were conducted on-site using Paganin's algorithm [3] and a modified Gridrec algorithm [19]. Additionally, FBP reconstructions of the data were produced at Northwestern University following the experiment, and used in the comparisons below. Each resulting dataset is 1525x1525x1598μm, with a voxel size of 0.74x0.74x0.74μm. Figure 2 shows the PAG reconstructions at various coarsening times.

Multimodal reconstruction technique
The phase map, φ (r 1 , r 2 ; θ ), of each projection is calculated during the phase-retrieval step of the PAG technique [3,18,20]. Taking the inverse radon transform [5] of φ (r 1 , r 2 ; θ ) gives the backprojected image, μ PAG (x, y, z), While it is easy to qualitatively differentiate between components in this δ -map, it is challenging to quantitatively characterize the system due to its low pass characteristics. On the other hand, when the projections are backprojected without the intermediate phase-retrieval step, the reconstructed image intensity, μ FBP (x, y, z), is such that [4] μ FBP (x, y, z) ∝ ∇ 2 δ (x, y, z) + μ atten (x, y, z) + μ mixed (x, y, z) where μ atten (x, y, z) is the linear attenuation coefficient given by and μ mixed (x, y, z) is a function of μ atten (x, y, z) and δ (x, y, z) [4]. The dominant term in Eq. (2) is the Laplacian of refractive index decrements, ∇ 2 δ (x, y, z), since the β (x, y, z) term in Eq. (3) does not provide appreciable contrast in the images. While the FBP image offers sharper interfaces due to ∇ 2 δ (x, y, z), the grayscale intensity levels of the components are very similar. A hybrid PCT reconstruction, μ + (x, y, z), is a linear combination of the PAG and FBP images, Here 0 ≤ c 1 ≤ 1, thereby combining the strong contrast present within the PAG image and the sharp interfaces found in the FBP image. In other words, the FBP image is a natural source of image sharpening. Thus, it is possible to extract two sets of data, PAG and FBP, from a single PCT experiment, and the linear combination of the two provides a hybrid reconstruction crucial for quantitative analysis. It is anticipated that this multimodal approach could be generally applicable to weakly absorbing samples (in which β ≈ 0) imaged in the near-field regime; these conditions would then give rise to the edge-enhancement of μ FBP and low-pass characteristics of μ PAG . However, the parameter c 1 , reflecting the contribution of the PAG image in the hybrid reconstruction, may be different and require sample-specific tuning for other datasets.

Multimodal image analysis
The scalar c 1 in Eq. (4) is determined by optimizing μ + (x, y) with respect to two image quality metrics: contrast-to-noise ratio (CNR) and sharpness (SH). The notation (x, y) indicates a 2D slice of the 3D (x, y, z) volume. In other words, whereas Ref. [14] tuned the propagation distance for optimal image quality, we optimize the relative contributions of PAG and FBP images in μ + (x, y), at a single sample-to-detector distance. In this way, varying the scalar allows for robust segmentation of our PCT images. CNR is defined as where S and σ are the mean and standard deviation, respectively, of the pixel values in the foreground, f , and background, b, regions. For our images, the foreground refers to the Si laths and the background refers to the eutectic matrix. CNR was determined by manually tracing over the interfaces of the Si laths in a representative number of 2D images. Although technically SH lacks a precise definition, intuitively, sharpness is related to the fineness of the resolvable details. Ref. [21] developed an algorithm to determine the overall SH of an image; we use their global single parameter sharpness model, implemented as the ratio between the output energy of an ideal high pass filter and an ideal band pass filter [21]: whereξ = (ξ x , ξ y ) are the Cartesian frequency coordinates, and H and B are the high and low-band pass frequency ranges, respectively. Additionally, the resolution of the multimodal images is cross-checked using the method described by Ref. [22]. First, the power spectral density (PSD) of an arbitrary line profile in the image is computed. This PSD converges to a value defined as the noise baseline. Resolution is computed by taking twice the value of the PSD at the noise baseline, and matching it to a spatial frequency, k res [22]. Then, the spatial resolution, x res , is calculated as Practical implementation of this resolution criterion is met if the maximum spatial frequency of the image is less than one-half of the sampling frequency of the line profile (i.e., the test image is over-sampled) [22]. Once the optimal scalar c 1 is determined by the above mentioned metrics, the resultant multimodal image, μ + , is segmented according to the procedure described in subsequent sections.

Segmentation of hybrid images
Our work aims to characterize the evolution of primary Si laths in an Al-Si system. Advanced image processing techniques, such as the Expectation Maximization/Maximization of Posterior Marginals (EM/MPM) algorithm [23,24], which was recently applied to materials datasets, fails to segment PCT images using their current image models. As a result, this approach does not allow us to characterize the primary Si laths that require the removal of the smaller fluctuations in the matrix. These fluctuations are a result of the eutectic lamellae that appear with the same intensity level as the primary laths, and are an artifact of the quenching process. Moreover, it is necessary to enhance the edges of the primary Si laths in order to obtain a robust segmentation.
Blurring the eutectic constituent and enhancing the interfaces of the primary Si laths can be accomplished by using a nonlinear diffusion filter, such as a regularized Perona-Malik filter (RPM) [25,26], which is applied on all 2D (x, y) slices of the 3D (x, y, z) dataset. In this method, a filtered image u(x, y,t), where t is time, is obtained as the solution of the diffusion equation The notation D |∇u σ (x, y,t)| 2 is the nonlinear diffusion coefficient, κ is the gradient threshold parameter, K σ is Gaussian structuring element with standard deviation σ , and ⊗ denotes the convolution operation [25]. The input of the algorithm is the hybrid, multimodal image, see Eq. (11). The iterative convolution of the image u(x, y,t) with K σ , in Eq. (10), regularizes the Perona-Malik model such that RPM is robust against local noise at scales smaller than or equal to σ [25]. This means that gradients that result from lamellae are effectively removed, given that they are smaller than the Gaussian kernel. To minimize the local fluctuations in the eutectic, the images are also pre-processed with a combination of erosion and dilation operations, and median filtering. The gradient threshold parameter, κ, is commonly fixed at a user-defined value [26]. A fixed κ that is too small, however, may misinterpret large gradients due to noise as edges it should preserve, while a κ that is too large may delete edges and small structures during the diffusion process [27]. Thus, κ is updated such that it is proportional to the average noise in u(x, y,t) at any given time. Noise can be estimated from morphological filters, see Refs. [28,29]; however, these morphological operations, e.g., erosion and dilation, are computationally expensive and thus we estimate average noise as where ω is a constant. Equation (12) suggests that the noise in the image drops exponentially under many applications of the RPM filter; this is consistent with the decaying exponential form in Eq. (9). In this way, the gradient threshold parameter κ self-adapts to the image u after every iteration.

Bias-corrected fuzzy c-means
The RPM filtered images are characterized by uneven background illumination. Intensity inhomogeneities or the bias field, arise from the extrinsic diffusional characteristics of RPM filtering as well as the intrinsic dark-light fringes suggested by Eqs. (2) and (4). This effect often results in under or over segmentation for fixed histogram threshold values. In order to estimate the bias field of an image, Ref. [30] introduced an algorithm based on fuzzy logic, known as bias-corrected fuzzy C-means (BCFCM). In particular, they modified the objective function, J m , of the standard fuzzy C-means algorithm as where c is the number of clusters, N is the number of pixels, x k is the k-th pixel of measured data, u ik is the degree of membership of x k in cluster i, p is the fuzziness coefficient, ν is the cluster prototypes or centroids, N k is the set of neighbors of x k , N R is the cardinality of N k , and α is a weighting parameter [30,31,32]. The notation * is the distance between x k and centroid ν i . More specifically, the regularizing effect of a pixel's local neighborhood is controlled by the parameter α [30]. Thus, the goal of BCFCM algorithm is to divide the data into two clusters, that of the primary silicon laths and that of the eutectic constituent, while taking into account the slow-varying bias field of the images. BCFCM was originally developed for magnetic resonance imaging, though it can be applied to PCT data by modeling the RPM filtered image as where x k and y k are the true and observed intensities of the k-th pixel, respectively, and β k the bias field at the k-th pixel. Inserting x k in Eq. (13) and minimizing J m with respect to constraints on membership, u, leads to an expression for β k [30]. The bias-corrected image x k has a bimodal histogram and therefore can be robustly segmented using conventional histogram-based segmentation methods, such as Otsu's method [33].

Digital inpainting of voids
BCFCM algorithm assumes a two-phase system. Presence of voids in the microstructure, which appear as dark regions in the image, result in misinterpretation of bias field of the image. Thus, prior to the bias-field correction, it is necessary to 1. determine if a given image has voids; 2. camouflage any voids.
Hartigan's dip test (HDT) is a statistical measure of the deviation of a distribution from unimodality [34]. It is a useful tool since statistically significant voids manifest as a separate peak in the image histogram. For each image, if HDT determines that the histogram is not unimodal to a confidence level of 95%, the voids are inpainted by solving the steady-state diffusion equation with constant diffusivity (i.e., Laplace's equation), only within the void, see Refs. [35,36] for details. Following inpainting of the voids, the images are processed using RPM filtering and BCFCM algorithm, respectively. Figure 3 summarizes the steps involved in processing of PCT data using the proposed methodology. All the algorithms discussed are applied in 2D; however, the methods can directly be extended to work in 3D. The total time for data processing is approximately 10 hours using a single node on Quest, the supercomputer cluster at Northwestern University, for a  296x296x159μm (34.6x10 6 voxels) stack of grayscale images. The node contains two Intel Nehalem Quad Core Xeon processors rated at 2.26 GHz. All codes are written in MATLAB R2012a [37]. It should be noted that the work of this paper is only to illustrate the proof-ofconcept of our highly integrated approach; as such, our execution time is only an upper bound on performance, and parallelization using a compiled language can drastically speed up convergence rates. where c 1 = 1. In general, the dark gray structures represent the primary Si laths and the light gray background along with the fine elongated darker features represent the eutectic. While there is sufficient contrast between Si laths and the eutectic in the PAG image in Fig. 4(e), the interfaces are quite diffuse, due to the low-pass characteristics of the single-image phaseretrieval algorithm. On the other hand, the FBP reconstruction in Fig. 4(a) features prominent dark-bright Fresnel fringes at the interfaces of the Si laths, which are due to ∇ 2 δ in Eq. (2).

Multimodal reconstruction
To assess the quality of the multimodal images, we measured CNR and SH for each image. Since we are interested only in relative values of CNR and SH for comparison purposes, normalized values, denoted byĈNR andŜ H, respectively, will be used whenever possible. Specifically,ĈNR in Fig. 4(f) is normalized with respect to the minimum and maximum values of CNR, at c 1 = 0 and c 1 = 1, respectively. It is also important to note thatŜ H is sensitive to the filter cut-off frequencies in Eq. (6). Thus, we investigate three high-pass frequency cutoffs, denoted k 1 > k 2 > k 3 in Fig. 4(f), where k 1 , k 2 , and k 3 are 0.32, 0.28, and 0.24 pixels −1 , respectively; allŜ H values are normalized with respect to the minimum and maximum values of SH at cutoff k 1 . In all cases, the band pass frequency range is defined as [0.02, 0.2] pixels −1 . As expected,Ŝ H increases with decreasing high-pass cut-off frequency, for a given value of c 1 , since the output energy of the high-pass filter increases. Regardless of the cut-off frequency selected, the results indicate that FBP images are sharper than PAG images, while PAG images offer greater contrast-to-noise ratios. The optimal reconstruction is a trade-off between these two image quality metrics, see Fig. 4(f). For all practical purposes, we use c 1 = 0.5 in our hybrid images. The resulting hybrid image preserves the contrast between Si laths and the eutectic, while also providing sharper interfaces for quantitative analysis.
Furthermore, spatial resolutions calculated using the power spectral density (PSD) approach are in good qualitative agreement with the relative sharpness estimates. Figure 5 shows the calculation of this resolution criterion, where |C(k)| 2 is the spectral power of the detected signal, μ s is the noise baseline, and k res is the maximum spatial frequency from Eq. (7). The PSDs shown in Fig. 5(a) and Fig. 5(b) have been calculated for line profiles in FBP and PAG images, respectively. It can be observed that the FBP image has a spatial resolution, x res , that is roughly 40% greater than that of the corresponding PAG image; as such, the FBP image provides a  natural source of image sharpening, but at a cost of CNR.
The enhanced resolution of the FBP image can be understood by exploiting the Wiener-Khintchine theorem [38], which states that the autocorrelation function is the Fourier transform of the power spectrum (i.e., autocorrelation and PSD are a Fourier transform pair). According to Fig. 5(a), then, the PSD of the FBP image shows high autocorrelation for long wavelengths. For this reason, k res for μ FBP is greater than k res for μ PAG . This correlation in pixel intensity values over large distances is expected since only high-frequency, interfacial fringes manifest in the FBP reconstructions.

Proposed segmentation method
For non-linear diffusion filtering, we choose κ 0 = 0.1, ω = 0.05, σ = 0.5, K σ to be defined as 10x10 pixels, and the number of iterations of RPM to 1,000. The parameter σ is set in reference to the standard deviation of grayscale intensity values in the eutectic. It is important to note that the quality of the filtered images is most sensitive to the gradient threshold parameter, see Fig. 6 for comparison of static and dynamic κ. Static refers to the fact that κ is fixed for all iterations, while dynamic indicates that κ is of the form given by Eq. (12). While static κ may lead to the smoothing of semantically important edges, and in the worst case, to the deletion of small structures during the diffusion process, dynamic κ preserves such edges, as previously discussed.
After approximately 1,000 iterations of the RPM filter, the image converges to a steady-state. Figures 7(a)-7(d) shows the effect of RPM filtering on the pre-processed image. Successive iterations of the RPM filter allow for the removal of intra-phase noise while still preserving interface positions. Interface width decreases with iterations of RPM filter because these stronger fluctuations are above the gradient threshold κ(t) and thus diffusion is inhibited. With RPM filtering, it is possible for diffusion and edge detection to interact in one process.
In the BCFCM algorithm, values of α = 1 × 10 −5 , p = 1.4, N R = 8, and 5,000 iterations were used to represent the image using two clusters, corresponding to the Si laths and the eutectic. The neighborhood parameter α is estimated as the inverse of the signal-to-noise ratio in the RPM-filtered image. Theoretically, the fuzziness coefficient p ∈ [1, ∞) [31], although the ideal value of p is problem-specific, and was empirically selected for our dataset. However, increasing p beyond 2 causes the clusters to overlap significantly such that cluster boundaries become ill-defined. Using these parameters, the bias-corrected image is presented in Fig. 7(e). Figure 8 summarizes the 2D segmentation steps. Figure 8(b) shows the inpainting of the void shown in the hybrid image ( Fig. 8(a)). Following inpainting of the voids, the images are processed using RPM filtering and BCFCM algorithm, see Figs. 8(c)-8(f). The evaluation of the segmentation results, e.g., Fig. 8(f), is difficult since the ground truth is unknown. However, the automated segmentation approach can be compared to manual segmentation in terms of the adjusted Rand index (ARI) [39,40]. ARI is a measure of similarity between the two methods, and varies between -1 and 1, where 1 indicates a perfect match. It is defined as where a i = ∑ j n i j , b j = ∑ i n i j , and n = ∑ i j n i j (16) and where n i j is the number of voxels belonging to class i in the manual segmentation, and to class j in the automated approach [39,40]. For our images, there are two classes corresponding to the Si lathes and the eutectic matrix. The notation * * represents the binomial coefficient. Manual segmentation was conducted by tracing out the interfaces of the hybrid images. To overcome experimenter's bias, we collected eight hand-generated segmentations of three hybrid PCT images from eight different people. Then, using Eqs. (15)-(16), we found that ARI = 0.91 ± 0.02. This ARI value is not intended to make a statement on the absolute accuracy of the automated segmentation, but to provide a relative comparison to the manual case. As such, the high ARI value indicates that using the proposed automated method, we can achieve results similar to that of manual segmentation with significantly less effort and much greater reproducibility. This segmentation technique takes into account the inherent structures and properties of the material being analyzed. For instance, the lamellar structure of the eutectic phase manifests as small (i.e., 2.5 ± 0.5 μm in width) intensity fluctuations in the hybrid images, while the Si laths are, on average, considerably larger. In order to isolate the Si laths from this finer lamellar structure, it is necessary to smooth the smaller fluctuations while enhancing the larger ones, through RPM filtering. This toolbox of algorithms for segmenting PCT images also has the potential to be tailored for other materials systems.

3D reconstruction
Once binarized, the 2D images are combined to reveal the 3D microstructure at different time steps. Figure 9 reflects the extraordinary morphological and topological complexity of the primary Si laths that evolve during coarsening. Characterization of such structures for comparison to coarsening theory requires a fully three-dimensional analysis (Fig. 9). The increase in size scale of the structure is consistent with a coarsening process, and an increasingly isotropic structure. For instance, Ref. [16] demonstrated, qualitatively, that the highly anisotropic Si laths do not evolve through a series of equilibrium Wulff shapes, and that there are interfaces with low mobility in the structure. For quantitative microstructural characterization of the coarsening morphologies in Fig. 9, see [41].

Conclusion
The proposed technique allows for the automated processing of PCT images. In agreement with Ref. [14], we find that there are advantages of near-field imaging, since it permits both phase-contrast and attenuation-based reconstructions of the same microstructure. The linear combination of these reconstructions, at a single sample-to-detector distance, offers the possibility of tuning the contrast-to-noise ratio and spatial resolution, thereby utilizing the advantages of both imaging modalities. The multimodal images can then be robustly segmented by using algorithms from computer vision, biomedical imaging, and art restoration communities; in particular, we use RPM filtering followed by BCFCM method to achieve binarized datasets. Finally, these datasets reveal the 4D microstructural evolution of an Aluminum-29.9wt% Silicon alloy during coarsening.
utilized the Quest high performance computing facility at Northwestern University, which is jointly supported by the Office of the Provost, the Office for Research, and Northwestern University Information Technology.