research papers
Insight into 3D microCT data: exploring segmentation algorithms through performance metrics
^{a}Computational Research Division, Lawrence Berkeley National Laboratory, Berkeley, CA 947208150, USA, ^{b}Berkeley Institute for Data Science, University of California Berkeley, Berkeley, CA 94720, USA, ^{c}Advanced Light Source, Lawrence Berkeley National Laboratory, Berkeley, CA 947208150, USA, ^{d}Materials Department, University of California Santa Barbara, Santa Barbara, CA 931065050, USA, and ^{e}Department of Mathematics, University of California Berkeley, Berkeley, CA 94720, USA
^{*}Correspondence email: tperciano@lbl.gov
Threedimensional (3D) microtomography (µCT) has proven to be an important imaging modality in industry and scientific domains. Understanding the properties of material structure and behavior has produced many scientific advances. An important component of the 3D µCT pipeline is image partitioning (or image segmentation), a step that is used to separate various phases or components in an image. Image partitioning schemes require specific rules for different scientific fields, but a common strategy consists of devising metrics to quantify performance and accuracy. The present article proposes a set of protocols to systematically analyze and compare the results of unsupervised classification methods used for segmentation of synchrotronbased data. The proposed dataflow for Materials Segmentation and Metrics (MSM) provides 3D microtomography image segmentation algorithms, such as statistical region merging (SRM), kmeans algorithm and parallel Markov random field (PMRF), while offering different metrics to evaluate segmentation quality, confidence and conformity with standards. Both experimental and synthetic data are assessed, illustrating quantitative results through the MSM dashboard, which can return sample information such as media porosity and permeability. The main contributions of this work are: (i) to deliver tools to improve material design and quality control; (ii) to provide datasets for benchmarking and reproducibility; (iii) to yield good practices in the absence of standards or groundtruth for ceramic composite analysis.
1. Introduction
Xray synchrotron facilities regularly produce terabytes of data, with imaging beamlines commonly storing data as twodimensional (2D) and threedimensional (3D) images (Bethel et al., 2015). The data volume generated daily and the variety of samples in terms of complexity and features is a challenge for effective data analysis of the experiments, particularly given frequent upgrades of the instruments' resolution and throughput. Specific characteristics of the image data such as the presence of heterogeneous structures in multiple scales and their complex architecture hinder the accuracy of current image processing algorithms. Nonetheless, microcomputed tomography (µCT) continues to be an essential imaging technique employed for the nondestructive 3D characterization of objects. This approach is widely used in academia and industry, including medical imaging, material science, electronics and geology.
In spite of the success of this imaging technique, some challenges still remain when analyzing these types of data, starting with the tomographic reconstruction and going through all image processing steps. For example, trying alternative acquisition schemes and experimental setups and/or quantitatively evaluating reconstruction methods are still a challenge. The work described by Ching & Gürsoy (2017) addresses this specific topic. The authors propose software that generates complex simulated phantoms and evaluates new or existing data acquisition schemes and image reconstruction algorithms for targeted applications. Following a similar strategy, the present paper addresses the problem of comparing and evaluating different image segmentation algorithms.
Given large data rates and sizes, machine learning for data acquisition (Yang et al., 2017) as well as for the automation of feature detection and extraction represents key steps in reducing data while gaining insight from µCT image structures. Many different categories of automated feature extraction exist (Hintermüller et al., 2010; Chen et al., 2012), but most of them are reliant on image segmentation techniques supported by unsupervised learning (Khanum et al., 2015). Broadly used to analyze experimental data, unsupervised segmentation algorithms (Chen et al., 2013) enable grouping picture elements, somewhat collecting tokens that `belong together'. Such algorithms can gather meaningful groups by searching for hidden structures from unlabeled data. Some of the advantages are dispensable training data and potential for data reduction, for example, removal of noncontributing image portions, such as background and artifacts.
Two important tasks during µCT data analysis are (1) selecting the best segmentation algorithm and (2) determining the best evaluation metrics. Those two choices highly impact the performance and quality of the results. Even harder to evaluate, the cases that lack groundtruth can lead to incomplete and/or ambiguous results, many times assisting only in qualitative terms. Therefore, we have employed a few strategies to evaluate image segmentation problems before we allude to the most appropriate algorithm for a specific application, for example, through direct comparison of algorithms based on a given groundtruth (Ushizima et al., 2011; Arbeláez et al., 2011; Perciano et al., 2016; Tassani et al., 2014; Sheppard et al., 2014) and combination of different algorithms (Polak et al., 2012) allied to confluence analysis as an indicator of result agreement.
In this work, we propose the Materials Segmentation and Metrics (MSM) dataflow that runs different algorithms separately, but checks for their agreement in terms of their segmentation results by using performance metrics. Suitability of metrics depends on the scientific goals of the experiment, that can be, for example, calculating porosity of a sample or counting the number of targeted objects. This article describes a multidisciplinary project involving experimental investigations of 3D µCT data of different materials performed at the Advanced Light Source (ALS) at the Lawrence Berkeley National Laboratory (LBNL), and development of algorithms for unsupervised image analysis, performed by investigators at the Center for Advanced Mathematics for Energy Research Applications (CAMERA), also at LBNL. We introduce a process for segmentation and analysis of 3D microtomography data, which offers strategies to evaluate algorithms and metrics applied to 3D microtomography experiments. Here, we use the proposed process to investigate geological samples and ceramic matrix composites (CMCs). These examples illustrate a strategy for executing different algorithms, extracting respective criteria, and parameter ranges to arrive at an appropriate answer.
Previous works on µCT data analysis described by Ushizima et al. (2011, 2012) address the problem of segmenting geological samples being studied to understand carbon sequestration, geologic storage of captured CO_{2} in underground rock formations. The authors developed tools for providing precise measurements of porosity and permeability, and for visualizing pore structures.
More recently, we have focused on yet another material sample: CMCs, which are composed of continuous silicon carbide (SiC) fibers and SiC matrices. With a high impact in industrial manufacturing, CMCs are an enabling element in the development of gas turbine engines that can operate at higher temperatures and therefore yield higher efficiency (Zok, 2016). Although there have been several recent studies related to computed tomography of CMCs, computational tools to automate and streamline image analysis are presently lacking (Bale et al., 2013).
In this article we: (a) introduce an analysis pipeline with different strategies for segmentation and evaluation; (b) demonstrate the use of this process in the detection of different structures in synthetic and experimental µCT datasets such as geological samples, glass beads and CMCs; and (c) make a critical assessment of the accuracy of the algorithms in determining various quantitative characteristics of the extracted structures using general and specific metrics. Fig. 1 presents the analytical dataflow of MSM. The figure highlights the transformations that the 3D image stacks undergo before MSM can deliver the materials characteristics list. The raw µCT image stack passes through a preprocessing step based on 3D nonlinear filtering. The preprocessed data serve as input for different segmentation algorithms [statistical region merging (SRM), kmeans, and parallel Markov random field (PMRF)]. The results from different segmentation algorithms undergo an analysis step based on both general and specific metrics. The evaluation takes into account reference data obtained from a semiautomated segmentation process.
The remainder of the article is organized as follows. The segmentation methods for sample partitioning and the metrics used are described in §2. In §3 we present different types of experiments using synthetic and experimental data to assess the accuracy of the results. In doing so, we show how the proposed process can be applied to evaluate different segmentation algorithms and to find the necessary metric in each case. Finally, §4 summarizes the key accomplishments of our investigation.
2. Materials and methods
This section describes the algorithms used in each step of the proposed process pipeline: image enhancement, image segmentation and measurements.
2.1. Image enhancement
Image enhancement techniques ensure that input data are well suited for a specific task, such as image segmentation or feature extraction. The challenge is to improve image quality while retaining essential information about the true structure. In the present study, we apply two main strategies: mathematical morphology (MM) (Pinoli & Debayle, 2012) and nonlinear edgepreserving filtering (Tomasi & Manduchi, 1998). Previously we developed `F3D', a graphicscard aware image processing plugin for Fiji (Schindelin et al., 2012) that employs MM and nonlinear filters and can handle data sets whose size exceeds the amount of RAM available in the computer system (Ushizima et al., 2014). Details on those techniques are out of the scope of this work.
F3D graylevel MM operators are onepass constanttime methods that can perform morphological transformations with structuring elements oriented in several directions. MM operators consist of two parts: (i) a reference shape or structuring element, which translates over the image, and (ii) a mechanism that defines the comparisons performed between the image and the structuring element (Van Droogenbroeck & Talbot, 1996). In this work, we use the closing operator, which is given by the combination of the dilation and the erosion operators (Gonzalez & Woods, 2006), to improve image contrast. The procedure involves the following steps:
(1) Define a structuring element based on the geometrical shape of the feature to be detected in the image.
(2) Apply the closing operator using the structuring element defined in step 1 [Fig. 2(b)].
(3) Subtract the input image from the result obtained in step 2 [Fig. 2(c)].
(4) Subtract the result obtained in step 3 from the input image [Fig. 2(d)].
Additionally, we apply a nonlinear image denoising filter with edgepreserving characteristics (Tomasi & Manduchi, 1998; Bethel, 2012). A weighted average of nearby pixels replaces each intensity value of the original image. This filter takes into account differences in intensity values in neighboring pixels to preserve edges while smoothing. Consequently, the influence between neighboring pixels depends on the similarity of their intensity values. It is defined as
where
I is the input image, and are spatial Gaussian kernels, p and q are pixel locations and S is the neighborhood.
The parameter , called range kernel, smooths differences within intensities, while the parameter , referred to as spatial kernel, smooths differences within coordinates. When the value of increases, the filter comes closer to a Gaussian convolution. The value of the spatial kernel is directly proportional to the size of the features to be smoothed. Fig. 3 shows an example of a bilateral filter applied to an image varying both parameters.
Upon application of the enhancement process, the images present improved contrast and reduced noise, making them more suitable for use as inputs to the image segmentation algorithms, described in the next section.
2.2. Unsupervised segmentation
In this section we present the three algorithms that are currently part of MSM: SRM, kmeans and PMRF. SRM and kmeans are among the most commonly used image segmentation strategies, yet due to the sensitivities of these methods the results can vary depending both on the dataset and parameters used. Additionally we evaluate PMRF, a novel graphbased algorithm, as another candidate to produce equivalent results, but addressing performance issues that traditional SRM and kmeans lack.
2.2.1. Statistical region merging (SRM)
Nock & Nielsen (2004) proposed an efficient region growing segmentation algorithm based on adaptive statistical threshold merging predicate on intensity levels. As in other region segmentation algorithms, it aims at associating a pixel to a region using a similarity criteria.
The iterative process starts with one region per pixel, followed by merging phases after the calculation of a statistical test that takes neighboring regions into account. This test considers an ascending order of intensity differences and checks if the mean intensities are sufficiently similar to be merged. The merging predicate, , regulates whether the observed regions R_{i} and R_{j} belong to the same statistical region. The merging predicate assumes that the pixels from a statistical region have the same expectation, and it is represented as
where the righthand side of the equation is the center value between R_{i} and R_{j}, and it is used as a merging threshold. The variable b is a function of g (the largest possible intensity value):
S_{l} is the set of regions with l pixels, δ is the probability error and takes values in . Q stands for the number of random variables, which somewhat translates the complexity of the scene and controls the coarseness of the segmentation. In other words, Q roughly estimates the number of regions in the image. The function represents the number of elements in the set of pixels. Even though the default µCT image output is 32bit, we propose the use of g = 255, therefore the methods take images in 8bit as input, given that we empirically verified that this is enough bit depth to represent more than 95% of the relevant intensity values from the original image.
2.2.2. kmeans
The kmeans algorithm (Macqueen, 1967) is a nonhierarchical unsupervised clustering method that classifies the input data points into k classes based on the inherent distance between point pairs. The algorithm assumes that the data features form a and tries to find natural clusters, iteratively minimizing the distance between the points and a set of centroids , . The minimization function is given by
where there are k clusters S_{j} , i = , and is the centroid of all points . The general kmeans algorithm for image segmentation consists of the following steps:
(1) Compute the histogram of pixel intensities.
(2) Initialize the centroids with k random intensities.
(3) Repeat the following steps until the cluster labels of the image converge:
(a) Cluster the pixels based on the distance of their intensities from the centroid intensities,
add i to S_{j}.
(b) Compute the new centroid of each cluster,
where k is the preferred number of clusters (usually set empirically and higher than the number of phases), i iterates over all intensities, j iterates over all centroids and are the centroid intensities.
2.2.3. Parallel Markov random field (PMRF)
Perciano et al. (2016) proposed a graphbased model called PMRF, which exploits the Markov random fields (MRF) framework to segment images. PMRF makes use of the linear and parallel (LAP) method (Mizrahi et al., 2014), a graph partitioning algorithm for MRF parameter estimation. In a MRF model, the optimization process uses a global energy function to find the best solution to a similarity problem, such as the best pixel space partition or the best matching. The energy function consists of a data term and a smoothness term. For image segmentation, we use the mean of the intensity values of a region as the data term. The smoothness term takes into account similarities between regions. The goal is to find the best labeling for the regions, so that the similarity between two regions with the same labels is optimal for all pixels (Mahapatra & Sun, 2012).
Given an image represented by = , where each y_{i} is a region, we seek a configuration of labels = , where and L is the set of all possible labels, L = . The MAP criterion (Li, 2013) states that one wants to find a labeling that satisfies
which can be rewritten in terms of the energies (Li, 2013) as
The prior probability is a Gibbs distribution, and the joint probability distribution is
where is a Gaussian distribution with parameters = and = is the parameter set.
The general PMRF segmentation framework works as follows. Initially, the input image passes through a feature extraction algorithm or an oversegmentation method to transform the voxeldomain image data into a less noisy representation. Next, the resulting regions compose the nodes of a graph representation of the image. The graph partitioning process is performed using the LAP algorithm, allowing simultaneous parallel parameter estimation and optimization for each subgraph. This parallelization strategy makes PMRF scalable, i.e. suitable for large experimental datasets. Finally, the iterative optimization process aggregates the graph nodes reaching an optimal segmentation through expectation maximization (EM) and maximum a posteriori (MAP) calculations.
2.3. Quantification metrics
In this section we describe metrics used to assess the accuracy of the unsupervised segmentation methods as applied to the µCT datasets.
2.3.1. General metrics
Our pipeline uses two sets of general evaluation measurements: binary segmentation metrics and material metrics.
The binary segmentation metrics consist of precision, recall and accuracy. Precision represents the proportion of voxels correctly classified as material, i.e. measures the performance of the algorithm with respect to false positives. It is given by: Precision = TP/(TP + FP), where TP stands for true positives and FP for false positives. Recall measures the performance of the segmentation algorithm with respect to false negatives (FN), and it is given by: Recall = TP/(TP + FN). Finally, accuracy gives the proportion of true segmentations among the total number of voxels. It is given by: Accuracy = (TP + TN)/(TP + TN + FP + FN), where TN stands for true negatives.
We calculate two additional metrics related to the material of the sample: the volume of the solid component of the sample, shown in cm^{3}, and the porosity of the material, given by = V_{v}/(V_{v} + V_{s}), where V_{v} is the volume of the void space (empty space) and V_{s} is the volume of the solid component.
All the metrics described in this section are calculated with respect to the reference data of each analyzed sample.
2.3.2. Characteristics of fiber beds
To evaluate the CMC samples containing fiber beds, in addition to the general metrics described in the previous section, we use more specific metrics to obtain additional scientific information about the samples. In a ) defines these cells. Fig. 4 shows an example of a Voronoi tesselation for a small region of fibers. The collection of cells can be interrogated in order to quantify various characteristics of the fiber bed. We compute three such characteristics for each through the stacks:
of a 3D µCT stack for these samples, each fiber is contained within a polygonal cell whose boundaries are defined by the perpendicular bisectors of lines joining the centroid of that fiber with the centroids of its nearest neighboring fibers. A Voronoi tessellation (Aurenhammer, 1991(i) The mean cell area, , given by
where is the area of one cell and is expressed in units of pixels.
(ii) The nonuniformity of cell areas (Shou et al., 2015), characterized by
where = is the entire set of cell areas.
(iii) The mean porosity, given by
where = is the set of unoccupied areas within the cells.
The α metric was proposed by Shou et al. (2015) for the analysis of permeability of the fiber reinforcement, and applied to simulated random fiber arrays. We propose to use the others to enrich the analysis of the samples. These metrics are calculated in 2D for each of the 3D stack.
Generally, the 2D cross sections used are the transverse slices of the 3D stack. However, the domain of analysis can be changed, i.e. the user can slice the data in different directions so that the metrics are calculated based on the selected domain. This flexibility allows taking into account that, due to the nature of a 3D tomographic reconstruction, characteristics and artifacts can vary depending on how the data are sliced.
3. Experiments
In this section we present a set of different in silico experiments using the proposed analysis pipeline. The experiments are carried out in an increasing order of complexity: (1) synthetic data corrupted with noise, (2) experimental data of geological sample, and (3) experimental CMC sample with fiber beds. We show how to use the proposed process to evaluate the segmentation algorithms and suitability of different metrics presented in §2.3 for each case.
3.1. Synthetic samples
3.1.1. Description
For the first set of experiments we use synthetic data corrupted with artifacts. In doing so, we have total control of the samples and the groundtruth to demonstrate the accuracy of the segmentation algorithms. We use benchmark images available at the Network Generation Comparison Forum (NGCF) (https://people.physics.anu.edu.au/~aps110/network_comparison), a forum that provides binary images with geometry that resembles porous media, which we contaminate with several levels of noise, similarly to Ushizima et al. (2012).
To produce images similar to the images obtained from tomographic experiments, we simulated a tomographic synchrotron experiment for each phantom image, and computed tomographic reconstructions. The ASTRA toolbox (Aarle et al., 2015) simulates tomographic projections for each phantom using 1025 angles over a 180° range. Afterwards, we added various sources of noise to the projections, simulating common artifact sources in practice, such as Gaussian white noise, miscalibrated detector pixels and shifting illumination, which typically result in ringlike artifacts in reconstructed images. The simulated projections were processed similarly to actual experimental data by applying a ringremoval algorithm (Münch et al., 2009) and computing a reconstructed image using the popular filtered backprojection algorithm (Kak & Slaney, 2001).
3.1.2. Analysis pipeline
For the synthetic µCT image stacks, we apply F3D and nonlinear image denoising to enhance image quality and then the three segmentation algorithms detect the phases of interest. Here, the algorithms separate two phases: structure and void space. To calculate all the measures described in §2.3.1, a voxeltovoxel comparison is made against the reference data (groundtruth) for each sample.
Empirically found values used for the parameters and by the edgepreserving filtering were 50 and 5. PMRF converges after five iterations on average, the SRM algorithm used Q = 8 and kmeans used k = 6.
3.1.3. Results
The first synthetic data simulate a glass bead column with µCT artifacts as shown in Fig. 5(a). Fig. 5(b) presents the reference data and Fig. 5(c) shows the 3D rendering of the corrupted data masked with the reference segmentation.
We apply the segmentation algorithms following the analysis pipeline described before. Figs. 5(d)–5(f) present the results for the 46th slice using PMRF, SRM and kmeans, respectively. The colors blue, red, yellow and black represent TP, FP, FN and TN.
Table 1 summarizes the quantitative results using the general metrics described in §2.3.1 for this dataset. Note that the quantitative results are compatible with the visual results shown in Fig. 5, i.e. the three algorithms perform similarly with small variations, and all of them approach the reference data. When these different algorithms and respective metrics somewhat coincide, our hypothesis is that the µCT segmentation worked properly.

The second experiment aims to evaluate the segmentation algorithms applied to a material with peculiar characteristics, including complex geometries and high curvature points. These synthetic data represent a geological formation similar to rocks. In this case, the material has a lower porosity, smaller and more connected structures. Additionally, the same approach for adding artifacts is applied to the synthetic data.
Figs. 6(d)–6(f) depict the segmentation results compared with the reference data, presented in Fig. 6(b). Fig. 6(a) shows the synthetic corrupted data and Fig. 6(c) presents the 3D rendering of the corrupted data masked with the reference segmentation. The colors blue, red, yellow and black represent TP, FP, FN and TN.
Note again the visual similarity of the results. Table 2 summarizes the general metrics values obtained for this experiment. Despite the fact that the simulated material in this case has a more challenging set of characteristics as described before, the results of the three algorithms also approach the reference data, with precision achieving values higher than 0.99, which is often enough to enable scientific interpretation of the sample.

3.2. Glass beads
3.2.1. Description
This real dataset presents the observation of a geological carbon sequestration experiment, which uses glass beads as a proxy for sand grains for in situ studies of infiltration and storage. Fig. 7 illustrates the experimental outcome using glassbeadpacked columns in biogenic mixture where calcite precipitation was induced by using the microbe S. pasteurii. The glass beads data set originally contains a 10 GB image stack and represents a 4.49 mm field of view, with a smaller core region corresponding to x × y dimensions of 1393 × 1398, as described previously by Ushizima et al. (2011).
3.2.2. Analysis pipeline
Following the procedure used for the synthetic datasets, firstly the F3D plugin combined with edgepreserving filtering enhances image quality and then the three segmentation algorithms detect the phases of interest. Here, besides separating two phases, material and void space, we apply the algorithms to find an additional phase representing the precipitation of calcium carbonate. The general metrics described in §2.3.1 are also used to evaluate the results on this dataset.
The values used for the parameters and by the edgepreserving filtering were 50 and 5. PMRF converges after five iterations on average, the SRM algorithm used Q = 8 and kmeans used k = 2 for the twophase segmentation and k = 6 for the threephase segmentation.
3.2.3. Results
We apply the pipeline described before to the experimental data of the glass bead column firstly aiming to separate material from void space. Fig. 7 presents a summary of the results for the first slice of the 3D stack. In the original data presented in Fig. 7(a), three phases can be identified by the difference of intensity values of the pixels: void space (darkest gray level), glass beads (intermediate gray level) and precipitate (brightest gray level). Fig. 7(b) presents the used reference for this slice for the twophases segmentation. Fig. 7(c) shows the 3D rendering of the original data masked with the reference data. Finally, Figs. 7(d)–7(f) present the results obtained by PMRF, SRM and kmeans, respectively. For the purposes of this work, the reference data used are the best results found in the literature for this dataset (Ushizima et al., 2011), which were validated by a material scientist expert. We observe a close agreement among the results compared with the reference data according to the values of TP, FP, FN and TN represented by the colors blue, red, yellow and black, respectively.
Table 3 presents the values for the metrics described in §2.3.1 calculated for this experiment. The results of the three algorithms are in close agreement with the reference data, and the accuracy metric achieved values higher than 0.94 indicating that the algorithms successfully segmented the solid component of the material from the void space.

Now we apply the same algorithms aiming to separate the three different phases: void space, glass beads and precipitate. Using the same previous slice, Fig. 8 shows the results for the three algorithms. In this case, the algorithms are able to separate the void space (black) and the glass beads (red), and obtain a phase which is a combination of precipitate and vessel (cyan).
Table 4 presents the volumes calculated for each phase. The reference data for the threephases segmentation are unavailable; however, we can notice the similarity among the results also indicating total agreement of the three algorithms, making the results valuable for visual analysis and characterization of this kind of material. As mentioned before, the phase representing the precipitation is combined with the vessel, which affects the calculation of the metrics for this phase. However, from an experimental point of view, it is possible to have the initial state of the experiment, i.e. the state of the glass bead column before infiltration. Consequently, it is feasible to subtract the initial state in order to obtain a phase which purely represents the precipitate.

The values obtained for the general metrics are calculated based on the full 3D binary results from the segmentation algorithms. The boundary shown in Fig. 7 was used for visualization purposes only, it is not a requirement for the metrics calculation. However, it is important to emphasize that, if a specific boundary is selected from the original image (a region of interest), then the calculated values for the metrics are going to vary accordingly.
3.3. Ceramic matrix composite
3.3.1. Description
Our last set of experiments focuses on the analysis of CMC samples containing fiber beds. The specimens analyzed in this work are unidirectional minicomposites. Each specimen contains approximately 5000–6000 aligned SiC fibers, each about 13 µm in diameter, surrounded by a cured preceramic polymer matrix. The matrix had been introduced through pressureassisted axial infiltration of a liquid polymer precursor into dry fiber beds encased in glass capillary tubes. The uniformity of fiber packing is of interest because it influences void formation during polymer impregnation and Octopus software (Inside Matters NV, 2016). Reconstruction required approximately 2 h when running on a standard desktop computer. The reconstructed results consist of sets of cross sections, recorded as 2D images. Brought together, the images can be used to recover the 3D structures.
(PIP) processing. µCT imaging of the infiltrated specimens was conducted at beamline 8.3.2 at ALS. The distance between the sample and the detector was 1.5–2.0 cm. The data were reconstructed using filtered back projection withFig. 9 presents a from one of the nine stacks studied. Three different regions are evident in the interior of the tubes: (i) fibers (lightest gray); (ii) cured precursor (medium gray) and (iii) voids (darkest gray). In this work, the precursor and the voids are together distinguished from the fibers.
3.3.2. Analysis pipeline
The segmentation and analysis pipeline were applied to nine 3D µCT image stacks, each approximately 2000 × 2000 × 2160 voxels. For each stack, F3D was first used to enhance image quality. Then, the image segmentation algorithms recover binary masks representing the separation of two phases: fiber beds and void space (combination of cured precursor and voids). From the binary masks, we calculate the general metrics described in §2.3.1.
Additionally, we use the binary masks to construct the Voronoi tesselation, which, combined with the segmentation results, enables calculation of the three specific metrics in equations (11), (12) and (13). The parameter Q for the SRM algorithm is set to 8, kmeans was applied with k = 4 and PMRF executes an average of five iterations. Automated threshold estimators are applied to obtain the binary results from SRM (Otsu) and kmeans (Isodata). The values used for the parameters and by the edgepreserving filtering were 50 and 5.
The results from the three segmentation algorithms are compared with corresponding results of the reference segmentation. The latter segmentation was obtained using MATLAB scripts with numerous parameters optimized for these specimens. The parameter values obtained from the reference segmentation are deemed to be the best estimates of their true values.
3.3.3. Results
Fig. 10 presents the segmentation result for a region of interest from a 2D plane of a stack. In this set of results, differences among the segmentation results can be observed, indicating that the algorithms are not in agreement and further analysis is necessary.
Following the same procedure used for the previous experiments, firstly we analyze the results using general metrics. Fig. 11 presents the graphics for the porosity and volume obtained for the nine stacks. The graphics show the curves in increasing order for both metrics. There is a weak indication that the PMRF algorithm approaches the reference results slightly more.
Table 5 presents the average segmentation performance of the algorithms applied to the nine stacks, which is in agreement with the other metrics results, where the PMRF algorithm achieves a slightly higher average accuracy. In general, it is clear that the three algorithms perform well considering these metrics. It is possible to have an overall understanding of the samples. However, it is difficult to draw any precise conclusions from these results regarding the fiber beds.

Taking a closer look at the results in Fig. 10, the fiber beds obtained by the kmeans algorithm suffer from missing fiber regions due to inaccurate pixel segmentation [note the spread yellow regions in Fig. 10(e)]. In contrast, being graphbased, both the SRM and the PMRF algorithms allow use of higherlevel information, correctly merging similar homogeneous regions.
However, the SRM algorithm is negatively impacted by the final threshold, leading to occasional missing or misclassified fiber regions [note the yellow regions in Fig. 10(d)]. So, in fact, there are more details about the results that the general metrics are not able to capture. Consequently, we need more specific criteria for a more precise evaluation.
The three fiber bed characteristics defined in §2.3.2 are now used. Fig. 12 shows the variation in the average value of each characteristic across each sectioning plane with position in the stack. The result based on the PRMF method is consistently in close agreement with the one from the reference segmentation. In contrast, the result from kmeans differs considerably from the reference values, overestimating both the nonuniformity of cell area, α, and the porosity, ρ, while underestimating the cell area, .
Table 6 summarizes the errors in the three computed characteristics from the three segmentation methods for all nine data sets. The PMRF method consistently yields the best results. For example, the error in α is a mere 0.0096. (For reference, α = 1 corresponds to uniform fiber packing; values obtained in the present composite specimens are typically about α = 1.1.) The errors in and ρ from this method are 0.012 and 0.0073, respectively. (For reference, the average values of and ρ are about 430 and 0.28.) Higher errors (in some cases by an order of magnitude) are obtained for the other methods.

Note that the results obtained using the specific metrics are different and more clear from the ones obtained using the general metrics. The reason behind this is that the specific metrics capture details of the fibers individually for each 2D plane throughout the 3D stacks, giving a more precise analysis of this structure. In doing so, the PMRF algorithm is in fact considerably more precise in this case than the other two algorithms, meaning that this algorithm appears to be more suitable to the analysis of this kind of sample. The main reason behind the better results obtained by the PMRF algorithm is its ability to estimate and separate more precisely the regions of interest based on a contextual model, not only on the intensity values of the pixels.
The set of experiments carried out was only possible and feasible given the proposed pipeline. The methodology described here uses different strategies for segmentation and evaluation of µCT samples, and also enables the critical assessment of the accuracy of the algorithms in obtaining scientific information about the samples. The results demonstrated the impact of choosing the suitable segmentation algorithms and evaluation metrics when analyzing 3D µCT data, and how the proposed pipeline facilitates this exploration and evaluation process.
4. Concluding remarks
Segmentation algorithms enable quantitative characterization of different materials in µCT data through measurements that extract physical properties from samples. However, finding the most appropriate segmentation algorithm for the data and finding the best evaluation metrics are a constant challenge, while those two choices heavily impact the performance of the results in obtaining valuable information from the data.
We have presented an automated pipeline called MSM for the analysis of 3D µCT images that offers different strategies to evaluate combinations of segmentation algorithm and metric for different scientific problems. We performed a broad and detailed set of experiments covering samples with different types of materials: synthetic data of glass beads and a rocklike sample, a real glass beads sample, and real CMC samples containing fiber beds. The impact of choosing the right segmentation algorithm and metric depending on the scientific goals is demonstrated throughout the experiments. Although MSM comes with a predefined set of segmentation algorithms and metrics, the pipeline can be easily extended with additional methods for comparison.
We show that the three unsupervised segmentation algorithms detect different phases of the materials and give comparable results for all the experiments. However, specimens with complex geometries and fine details may diverge considerably with regards to the segmentation results. The analysis of fiber beds requires tailored segmentation algorithms, and a more specific set of metrics so that the necessary scientific questions can be answered. In this case, the experiments indicate that the PMRF algorithm is more suitable and more precise when segmenting the fiber beds.
Future directions of this work include studies of scalability of the segmentation methods in order to improve the efficiency of the analytical process. This includes the use of an underdevelopment MPIbased version of the PMRF algorithm. In doing so, we intend to perform a detailed analysis of the fiber beds in CMC samples, modeling the fibers, and measuring fiber deformations and the matrix deformation evolution. Additionally, we will develop metrics that suggest whether or not the input data need enhancement. Currently, the image enhancement step in the proposed pipeline is optional and dependent on the user's choice.
Acknowledgements
This work was supported by: the Office of Naval Research, monitored by Dr David A. Shifler (see funding details below); the grant `Towards Exascale: High Performance Visualization and Analytics' and the Early Career Research project, both under Advanced Scientific Computing Research (ASCR), under program manager Dr Lucy Nowell (see funding details below); the Center for Applied Mathematics for Energy Related Applications (CAMERA), under management of ASCR and Office of Basic Energy Sciences (BES) of the US Department of Energy; an ALS Doctoral Fellowship in Residence (supporting NML) with funding provided by the US Department of Energy; and a National Science Foundation Graduate Research Fellowship (NML) (see funding details below). Any opinion, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. The experimental work was performed at beamline 8.3.2 at the Advanced Light Source (ALS), a Division of Lawrence Berkeley National Laboratory. The authors gratefully acknowledge Alastair MacDowell and Harold Barnard (ALS) for assistance with planning and executing the experiments. ALS is supported by the Director, Office of Science, Office of Basic Energy Sciences, of the US Department of Energy under Contract No. DEAC0205CH11231.
Funding information
Funding for this research was provided by: Office of Naval Research (award No. N000141310860); Advanced Scientific Computing Research (award No. DEAC0205CH11231); National Science Foundation Graduate Research Fellowship (award No. 1144085).
References
Aarle, W. van, Palenstijn, W. J., De Beenhouwer, J., Altantzis, T., Bals, S., Batenburg, K. J. & Sijbers, J. (2015). Ultramicroscopy, 157, 35–47. Web of Science PubMed Google Scholar
Arbeláez, P., Maire, M., Fowlkes, C. & Malik, J. (2011). IEEE Trans. Pattern Anal. Mach. Intell. 33, 898–916. PubMed Google Scholar
Aurenhammer, F. (1991). ACM Comput. Surv. 23, 345–405. CrossRef Google Scholar
Bale, H. A., Haboub, A., MacDowell, A. A., Nasiatka, J. R., Parkinson, D. Y., Cox, B. N., Marshall, D. B. & Ritchie, R. O. (2013). Nat. Mater. 12, 40–46. CrossRef CAS PubMed Google Scholar
Bethel, E. W. (2012). Exploration of Optimization Options for Increasing Performance of a GPU Implementation of a ThreeDimensional Bilateral Filter, Technical Report LBNL5406E. Lawrence Berkeley National Laboratory, Berkeley, CA, USA. Google Scholar
Bethel, W. et al. (2015). DOE ASCR Workshop, pp. 2–30. DOE. Google Scholar
Chen, R.C., Dreossi, D., Mancini, L., Menk, R., Rigon, L., Xiao, T.Q. & Longo, R. (2012). J. Synchrotron Rad. 19, 836–845. Web of Science CrossRef IUCr Journals Google Scholar
Chen, W., Ostrouchov, G., Pugmire, D., Prabhat & Wehner, M. (2013). Technometrics, 55, 513–523. Google Scholar
Ching, D. J. & Gürsoy, D. (2017). J. Synchrotron Rad. 24, 537–544. CrossRef IUCr Journals Google Scholar
Gonzalez, R. C. & Woods, R. E. (2006). Digital Image Processing, 3rd ed. Upper Saddle River: PrenticeHall. Google Scholar
Hintermüller, C., Marone, F., Isenegger, A. & Stampanoni, M. (2010). J. Synchrotron Rad. 17, 550–559. Web of Science CrossRef IUCr Journals Google Scholar
Inside Matters NV (2016). Octopus imaging, https://octopusimaging.eu/. Google Scholar
Kak, A. C. & Slaney, M. (2001). Principles of Computerized Tomographic Imaging. Philadelphia: Society for Industrial and Applied Mathematics. Google Scholar
Khanum, M., Mahboob, T., Imtiaz, W., Abdul Ghafoor, H. & Sehar, R. (2015). Int. J. Comput. Appl. 119, 34–39. Google Scholar
Li, S. Z. (2013). Markov Random Field Modeling in Image Analysis, 3rd ed. Springer Publishing Company. Google Scholar
Macqueen, J. (1967). 5th Berkeley Symposium on Mathematical Statistics and Probability, pp. 281–297. Google Scholar
Mahapatra, D. & Sun, Y. (2012). IEEE Trans. Image Process. 21, 170–183. CrossRef PubMed Google Scholar
Mizrahi, Y. D., Denil, M. & de Freitas, N. (2014). International Conference on Machine Learning (ICML), 21–26 June 2014, Beijing, China. Google Scholar
Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567–8591. Web of Science PubMed Google Scholar
Nock, R. & Nielsen, F. (2004). IEEE Trans. Pattern Anal. Mach. Intell. 26, 1452–1458. CrossRef PubMed Google Scholar
Perciano, T., Ushizima, D. M., Bethel, E. W., Mizrahi, Y. D., Parkinson, D. & Sethian, J. A. (2016). IEEE International Conference on Image Processing (ICIP), 25–28 September 2016, Phoenix, Arizona, USA, pp. 1259–1263. Google Scholar
Pinoli, J. C. & Debayle, J. (2012). IEEE J. Sel. Top. Signal. Process. 6, 820–829. CrossRef Google Scholar
Polak, S. J., Candido, S., Levengood, S. K. L. & Wagoner Johnson, A. J. (2012). Comput. Med. Imaging Graph. 36, 54–65. CrossRef PubMed Google Scholar
Schindelin, J., ArgandaCarreras, I., Frise, E., Kaynig, V., Longair, M., Pietzsch, T., Preibisch, S., Rueden, C., Saalfeld, S., Schmid, B., Tinevez, J.Y., White, D. J., Hartenstein, V., Eliceiri, K., Tomancak, P. & Cardona, A. (2012). Nat. Methods, 9, 676–682. CrossRef CAS PubMed Google Scholar
Sheppard, A., Latham, S., Middleton, J., Kingston, A., Myers, G., Varslot, T., Fogden, A., Sawkins, T., Cruikshank, R., Saadatfar, M., Francois, N., Arns, C. & Senden, T. (2014). Nucl. Instrum. Methods Phys. Res. B, 324, 49–56. CrossRef CAS Google Scholar
Shou, D., Ye, L. & Fan, J. (2015). J. Compos. Mater. 49, 1753–1763. CrossRef CAS Google Scholar
Tassani, S., Korfiatis, V. & Matsopoulos, G. K. (2014). J. Microsc. 256, 75–81. CrossRef CAS PubMed Google Scholar
Tomasi, C. & Manduchi, R. (1998). Proceedings of the Sixth IEEE International Conference on Computer Vision, Bombay, India, pp. 839–846. Google Scholar
Ushizima, D., Morozov, D., Weber, G. H., Bianchi, A. G. C., Sethian, J. A. & Bethel, E. W. (2012). IEEE Trans. Vis. Comput. Graph. 18, 2041–2050. CrossRef CAS PubMed Google Scholar
Ushizima, D., Parkinson, D., Nico, P., AjoFranklin, J., MacDowell, A., Kocar, B., Bethel, W. & Sethian, J. (2011). Proc. SPIE, 8135, 813502. Google Scholar
Ushizima, D. M., Perciano, T., Krishnan, H., Loring, B., Bale, H., Parkinson, D. & Sethian, J. (2014). IEEE International Conference on Big Data, 27–30 October 2014, Washington, DC, USA, pp. 683–691. Google Scholar
Van Droogenbroeck, M. & Talbot, H. (1996). Pattern Recognit. Lett. 17, 1451–1460. CrossRef Google Scholar
Yang, X., De Carlo, F., Phatak, C. & Gürsoy, D. (2017). J. Synchrotron Rad. 24, 469–475. CrossRef IUCr Journals Google Scholar
Zok, F. W. (2016). Am. Ceram. Soc. Bull. 95, 22–28. CAS Google Scholar
This article is published by the International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.