Quantitative analysis of X-ray microtomography images of geomaterials: Application to volcanic rocks

Analytical techniques commonly applied to study the rock textures are often restricted to two-dimensional (2D) measurements (see, for example, Ammouche et al., 2000). However, many research fields require the knowledge of the three-dimensional (3D) morphology and topology of the rock constituents. Several efforts have been made to extend 2D measurements to the third dimension, for example with stereological methods, but the information is mostly incomplete and strenuous to obtain (e.g., Stroeven et al., 2009). High-resolution X-ray computed tomography (μCT) supplies an effective method to directly acquire 3D information in a nondestructive way, providing a digital data set easily integrated with modern computer processing.

Beyond simple visualization, the main ambition of image processing is the quantitative morphological analysis of substantial data volumes. A number of parameters can be extracted from an X-ray tomographic data set if the spatial resolution is adequate to resolve the features of interest and the proper software is available, including geometry and organization of discrete rock components, such as crystals, vesicles, fractures, and zones with chemical alteration or with different composition.

In the geosciences (see Ketcham and Carlson, 2001; Mees et al., 2003; Carlson, 2006, for more applications), X-ray computed microtomography and 3D imaging have been used, for example, with metamorphic rocks to study porphyroblasts and make inferences about crystallization mechanisms (see, for example, Ketcham et al., 2005). With sedimentary rocks, soils, and sediments, interest has been directed toward the description of the 3D distribution of pores and fractures and the modeling of the fluid transport properties of rocks (see Van Geet et al., 2000). Additionally, as for paleontological samples, the nondestructive nature of μCT measurements has made it a valuable tool to image fine details of the internal structures of fossilized organisms (Tafforeau et al., 2006).

In the case of volcanic rocks, μCT can be used to image and to quantify the textural and morphological characteristics of the rock constituents, such as vesicles (former gas bubbles in solidified, erupted products), crystals, and glass fibers. For pyroclastic rocks, important parameters to characterize the vesicle portion are the pore size, size distribution, geometry and orientation, the pore-throat size and organization, the pore-surface roughness, and the topology of the overall pore and pore-throat network that are globally described by the tortuosity parameter (Wright et al., 2009).

Magma near the Earth's surface exists as a multiphase system, including gas bubbles and solid crystals in a liquid medium. The rheology of the magma and the processes that govern the transition between effusive and explosive eruptions can be fully understood if the gas permeability and flow through the bubble networks are quantified (Eichelberger et al., 1986; Klug and Cashman, 1996; Klug et al., 2002; Mueller et al., 2005; Burton et al., 2007; Namiki and Manga, 2008; Walsh and Saar, 2008; Bouvet de Maisonneuve et al., 2009). However, the dynamics of the magma flow in the conduit cannot be directly observed and thus requires indirect methods to be studied. The importance of pyroclasts in volcanology is due to the assumption that they are natural records of the magma state, in terms of texture and composition, during the last phases of conduit ascent. Most of the textural parameters are only slightly modified in the very last part of the pyroclasts’ life, due to successive bubble growth, coalescence, or collapse (Polacci et al., 2008, 2009a). Hence, the shape, size, orientation, and number density of vesicles and crystals possibly record parameters that are related to rheology, volatile saturation, flow and mechanical strain, and magma cooling rate (Rust et al., 2003).

In this paper, we present the results of the quantitative analysis of 3D images of volcanic rocks obtained by using X-ray μCT. The experiments have been carried out at the Elettra Synchrotron Light Laboratory in Trieste (Italy) and the reconstructed 3D images (or volumes) have been analyzed with a new software library, named Pore3D, custom-developed by the SYRMEP group at Elettra (Brun et al., 2010).

The purpose of this work is to present procedures available in Pore3D to extract quantitative information from μCT images of volcanic rocks. This information can be coupled to physical, rheological, and chemical properties of the parent magma and supply the bases to model and understand the eruptive dynamics.

Transmission X-ray μCT is a nondestructive technique based on the mapping of the linear attenuation coefficient of X-rays crossing the investigated sample (Baruchel et al., 2000a). When using conventional X-ray sources the beam is generally polychromatic and the axial 2D images (slices) reconstructed from μCT experiments in many cases display the so-called beam hardening artifacts (Baruchel et al., 2000a). These artifacts are due to the differential absorption of the X-ray spectrum by the sample and lead to a misleading recovery of the linear absorption coefficients, mainly appearing as brighter sample borders in the reconstructed slices. Because beam hardening is related to the density of the material and the geological samples are generally very dense, these artifacts can be strong.

Performing synchrotron X-ray μCT experiments at third-generation synchrotron sources, such as Elettra, highly improves the image quality, compared to conventional sources. Indeed, the X-ray beam used for synchrotron X-ray μCT has three main properties that significantly increase both the data quality and the imaging possibilities: the monochromaticity, the high intensity, and the high spatial coherence. A monochromatic high-intensity beam leads to reconstructed slices free of beam hardening effects and with a high signal-to-noise ratio. Moreover, thanks to the high transverse coherence length of the X-ray beam, it is possible to apply phase-contrast imaging techniques with a simple experimental setup. While in transmission X-ray imaging, the contrast is generated only by the absorption properties of the sample, in phase-contrast imaging it also arises from the interference between parts of the wave front that experienced different phase shifts. The free-space propagation between the sample and the detector allows the detection of X-ray interference patterns (Fresnel diffraction). They are associated with the sample features that induce phase modifications in the monochromatic, transmitted X-ray beam and that would be indistinguishable in absorption mode. The result is that in phase-contrast imaging, all the interfaces corresponding to abrupt phase changes in the sample show an enhanced visibility (Snigirev et al., 1995; Cloetens et al., 1996; Baruchel et al., 2000b). Compared with absorption mode images, the objects appear sharper and very small details can be detected (also smaller than the pixel size of the detector, see Cloetens et al., 1996; Mancini, 1998a). In this approach, in a μCT experiment, an important parameter is the sample-to-detector distance d. In the so-called edge detection mode (d << a²/λ, where a is the size of the investigated feature measured perpendicularly to the beam direction, and λ is the X-ray wavelength), the recorded images can be used to extract detailed morphological information from the sample (Cloetens et al., 1996; Cloetens et al., 1997; Mancini et al., 1998b). With modern cone-beam μCT systems, equipped with microfocus X-ray tubes, limited but clearly detectable phase-contrast effects can also be achieved, due to the small size of the focal spot (Wilkins et al., 1996).

The analysis of the volcanic samples was performed using the Pore3D software library. This software library allows a quantitative description of the morphology and topology of the sample components operating directly in the 3D domain, i.e., the 3D behavior is not inferred from the stacked 2D information. Pore3D is specifically designed and optimized for X-ray μCT images of porous media and multiphase systems, including geological materials, bones, biomedical compounds, or foamy structures (Brun et al., 2010; PORE3D, 2009).

Commercial software (e.g., MAVI-Fraunhofer ITWM; Avizo-VSG), publicly available libraries (e.g., DIPlib, Van Kempen et al., 2008; ITK, Ibáñez et al., 2003), and research codes (e.g., ImageJ, Abramoff et al., 2004; 3DMA, Lindquist and Venkatarangan, 1999; Blob3D, Ketcham, 2005b; Quant3D, Ketcham and Ryan, 2004) have been developed in recent years for 2D and 3D image analysis. Pore3D merges some of the features already offered by the above-mentioned software and adds new tools, e.g., for artifact reduction in the tomographic images (Kak and Slaney, 1987; Barrett and Keat, 2004). The main aim of this software project is to implement state-of-the-art methods for processing and quantitative analysis of 3D images. Pore3D includes the basic tools for performing prefiltering and segmentation of digital images and more refined procedures for analysis, such as skeleton, connected components, and anisotropy, as well as morphometric and fractal measures. These tools are organized in a modular architecture that through a choice of analysis modules leads the user to choose his/her own analysis strategy (Fig. 1).

The implementations of the algorithms in Pore3D were cross-checked with other software. For instance, numerical results computed by the anisotropy analysis toolbox of Pore3D were compared with those produced by publicly available software such as Quant3D and the code provided by Simmons and Hipp (1997). Similarly, standard histomorphometric analysis results were compared with those by the commercial software VG Studio Max, the open source software Microview (http://microview.sourceforge.net/), and the already cited code provided by Simmons and Hipp (1997). In the same way, the skeleton analysis was compared with the skeletonization tools available in the 3DMA package. The code for the computation of Minkowski functionals is provided by Ohser and Mücklich (2000).

The capability to interact with Matlab and IDL is another characteristic of the Pore3D software library allowing users to present and further process numerical results with these well-known software packages. The Pore3D engine is available for 32- and 64-bit computer architecture, exploiting hardware with more than 4 GB of RAM. This aspect is crucial since high resolution X-ray μCT produces huge data sets. For some procedures, code execution is faster on multicore and multiprocessor machines thanks to the adoption of OpenMP programming extensions. Future versions of the software will offer support for computer clusters. At present, Pore3D is only available to the Elettra users. In the near future, it will be offered as Software-as-a-Service to external users. More details about the software as well as access policy and the project evolution can be found at the official web page (PORE3D, 2009).

Generally, the first step of the image analysis process consists of the extraction of a volume of interest (VOI) having dimensions suitable for the available computing resources but preserving the sample representativeness. If the VOI includes multiscale heterogeneities of the investigated sample, and it is able to represent values of statistical significance averaging the measurement of highly variable properties, it is usually denoted as the representative elementary volume (REV). In practical cases, finding a REV can be difficult, especially with strongly nonhomogeneous samples or with insufficient spatial resolution (Degruyter et al., 2010).

After the extraction of the VOI, the volume is usually preprocessed with one or more grayscale-to-grayscale filters. These filters are used to ease further processing of the image by reducing noise and potential small artifacts, to enhance the edges of the objects, and to distinguish the different phases, corresponding to different chemical compositions (and/or density, in the case of X-ray CT). Next, the image is segmented, i.e., partitioned into sets of voxels (the 3D equivalent of pixels) corresponding to particular components of interest. Several approaches exist: simple methods which compare the grayscale values with a global threshold, more refined thresholding techniques with locally adaptive thresholds, and region-growing methods (Pratt, 2007). After the segmentation, a “cleaning” procedure may be performed in order to remove small defects and undesired or unphysical objects. Cleaning procedures are typically based on morphological (erosion and dilation) as well as topological (based on the object connectivity, as the object labeling) operators (Serra, 1982).

For the segmented volume, several structural, morphological, and textural parameters can be computed. For instance, if the investigated pore space is assumed to be composed by disconnected “blobs,” the connected components analysis allows separating, labeling, and counting each element. Then volume, shape, and characteristic orientation of each element can be measured. On the other hand, in the presence of an interconnected porous space, an analysis based on the concept of the skeleton is typically more suitable (Lindquist and Lee, 1996). It consists of tracking and studying what can be considered as the medial axis of the porous network. By applying one of the skeletonization algorithms, followed by a “pruning” processing if unnecessary parts need to be removed, several features of the complex skeleton can be extracted and measured, thus retrieving information about the topology of the porous space.

The samples analyzed in this study are the explosive products of different eruptions, with compositions spanning the basalt-to-trachyte range (see Table 1 for detailed information). The samples have been collected at volcanoes and/or volcanic areas characterized by eruptive styles from mild Strombolian to Plinian. We selected an assortment of specimens with different compositions and textural features to test the Pore3D performance and to possibly relate such features to the eruptive processes that generated the rocks.

Samples STR1 and STR2 represent, respectively, a scoria and a pumice from Stromboli 198 volcano (Aeolian Archipelago, Italy). This is a mildly explosive, steady-state active volcano, displaying a high-K basaltic composition and a normal, daily activity of intermittent explosions every 10–20 min with continuous noneruptive, quiescent degassing (Rosi et al., 2000). The basaltic scoria STR1 is the product of a normal-style explosion of May 2006, when highly porphyritic (mostly phenocrysts of feldspar, pyroxene, and olivine), poorly to moderately vesiculated high-K basalts were erupted (Polacci et al., 2009a).

Sample STR2 was collected after the 15 March 2007 paroxysm, when an eruption column formed as a consequence of a renewal in the volcanic intensity, and pumice products mingled with black scoriae were emitted (Landi et al., 2009). These products are crystal-poor, highly vesiculated pumices from a high-K, calc-alkaline basaltic magma.

Sample AGN is a pumice clast from the Agnano-Monte Spina eruption at Campi Flegrei (Neapolitan area, Italy): this event represents the highest-magnitude eruption of the past 5000 yr (De Vita et al., 1999). The sample has a trachytic composition and belongs to a pumice-fallout deposit (B1 unit of layer B). Macroscopically, it presents a clear anisotropy of the bigger vesicles, which are stretched coherently in the direction of the longest axis of the sample (tube pumice). The pyroclast is porphyritic with phenocrysts of feldspar and clinopyroxene as dominant phases.

Finally, sample AMB is a scoria from Ambrym volcano, a large basaltic volcano of the Vanuatu Islands. The sample, of basaltic composition, belongs to the past activity of the Benbow crater and was collected during a field campaign in October 2008. Macroscopically, it is characterized by a homogeneous distribution of isolated spherical vesicles of various sizes and very low crystal content, mostly with sizes at the submillimeter scale.

Synchrotron X-ray μCT experiments were performed at the SYRMEP beamline of Elettra (SYRMEP, 2010). The source is a bending magnet and the beamline provides, at a distance of ∼23 m from the source, a monochromatic, parallel, laminar-section X-ray beam with a maximum area of ∼160 mm × 6 mm. The monochromator is based on a double-Si (111) crystal system working in Bragg configuration and the energy E of the X-ray beam can be tuned between 8.3 and 35 keV with an energy resolving power on the order of 10⁻³.

The parallel-beam microtomography setup is described by Polacci et al. (2009a). For each sample, between 900 and 1800 radiographic images (or projections), determined by the sample lateral size and texture, were acquired with equiangular steps over a range of 180°, depending on the beam characteristics. The projections were recorded with a water-cooled, Photonic Science XDI-VHR charge coupled device (CCD) camera, consisting of a full frame CCD imager coupled to a gadolinium oxysulphide scintillator by a fiber-optic taper. The camera has a 12-bit dynamic range and a 4008 × 2672 pixel chip. The effective pixel size of the detector is (4.5 × 4.5 μm²). In our experiments, the pixel size has been fixed to 9.0 μm (2 × 2 binning), thus yielding a field of view of 18 mm (horizontal) × 12 mm (vertical). A monochromatic X-ray beam energy between 25 and 30 keV, depending on the sample absorption, provided suitable contrast; the duration of each microtomographic scan was on the order of 1 h. For all the investigated specimens, the sample-to-detector distance d was set to 200 mm (edge detection mode). The custom software Syrmep_tomo_project (see for details Montanari, 2003), based on a filtered back-projection algorithm (Kak and Slaney, 1987), was used to reconstruct the 2D slices from the sample radiographs. The slices were stacked to obtain volumes where the isotropic voxel has an edge size of 9.0 μm.

In addition to the samples analyzed by the synchrotron source, one sample (AGN) was imaged at the TOMOLAB station at Elettra (TOMOLAB, 2010; Mancini et al., 2007; Polacci et al., 2009a). The TOMOLAB is a μCT system with a polychromatic X-ray cone-beam, complementary to the SYRMEP beamline both for the energy range and the vertical size of the beam. The instrument is equipped with a Hamamatsu sealed, microfocus X-ray tube operating in a voltage range from 40 to 130 kV, with a maximum current of 300 μA and a minimum focal spot of 5 μm. The detector is a water-cooled Photonic Science XDI-VHR CCD camera, consisting of a full frame CCD imager coupled to a gadolinium oxysulphide scintillator by a fiber-optic taper. The camera has 12-bit dynamic range and a 4008 × 2672 pixel chip, providing a good combination between a small effective pixel size (12.5 × 12.5 μm²) and a large field of view (maximum 50 mm horizontal × 33 mm vertical). Exploiting the magnification effect offered by the cone-beam geometry (Feldkamp et al., 1984; Kak and Slaney, 1987), the source-to-sample (d₁) and source-to-detector (d₂) distances can be varied in order to achieve a spatial resolution close to the focal spot size of the source, on samples from a few millimeters to a few centimeters in size. For sample AGN the d₁ and d₂ distances have been set to 80 mm and 300 mm, respectively. This, combined with the 2 × 2 binning applied to the CCD camera pixels, results in an equivalent pixel size of 6.7 × 6.7 μm². From the tomographic projections acquired by the CCD camera during a 360° rotation of the sample, a set of 2D slices was reconstructed with the commercial Cobra-Exxim Computing software package. The ring artifacts present in the slices were reduced with an algorithm custom-implemented in the Pore3D software (Brun et al., 2009).

Preliminary volume visualization with volume rendering procedures is generally performed to qualitatively characterize the sample microstructure. In this work, the volume renderings were obtained with the commercial software VGStudio MAX 2.0. A more accurate quantitative study of the volcanic samples and their morphological parameters was then carried out with the software library Pore3D.

In Figure 2, 2D axial slices reconstructed from the central part of the investigated samples are shown, with the analyzed subvolume positions. In Figure 3 the volume renderings corresponding to the analyzed subvolumes of the different samples are presented.

The tomographic experiments on the selected pumices and scoriae have given us the opportunity to reconstruct and study the 3D internal structure of very different samples that originated at volcanoes with different eruptive behavior and hazard potential. In the following, to illustrate the code potential we present several of the Pore3D utilities, their importance, and how they have been applied to volcanic samples.

The AMB scoria (Fig. 2A) was chosen to particularly focus on porosity-related parameters, by processing a subvolume of 268 × 268 × 268 voxels, equivalent to a cube of 2.4 mm side length (Fig. 3A). The sample was selected to study the morphology and the topological characteristics of the voids (vesicles), such as their connectivity, because it contains a homogeneous distribution of spherical to subspherical vesicles of different sizes separated by glassy matrix walls.

The 3D image was first preprocessed with a grayscale-to-grayscale filter to reduce the noise and enhance its edges. We selected a bilateral filter (Tomasi and Manduchi, 1998), which is more refined and demonstrated better results than the more “classical” mean, median, and Gaussian filters. The resulting image was then segmented to obtain the binary distribution for the objects of interest. We adopted a manual segmentation process among several available choices (for example, Otsu filter, Otsu, 1979; and Niblack filter, Niblack, 1986) to optimally identify the vesicles. Quantitative analysis was carried out on this preprocessed image.

In the Basic Analysis module, the porosity parameter (i.e., the fraction of voids in the total sample volume) and the topological characteristics, as described by the Minkowski functionals (Steele, 2007; Ohser and Schladitz, 2009), were computed: specific surface area (surface area of all voids divided by the total VOI volume), integral of mean curvature (index of the dominance of convex or concave shapes, see later and Russ and DeHoff [1986]), and Euler characteristic (index of connectivity of the object network, see later and Odgaard and Gundersen [1993]; Ohser and Mücklich [2000]).

In the Texture Analysis module, the fractal dimension of the sample was calculated (Mandelbrot, 1982) to assess the self-similarity of order at different scales in a random arrangement of objects. The fractal computation was performed on the void fraction, which is supposed to be a fractal with the same dimension as the void-matrix interface, similar to sedimentary rocks (Katz and Thompson, 1985). The fractal dimension lies within the range of 2 to 3: smooth surfaces have a value close to 2, while higher values represent an increasing roughness and/or complexity of the objects (Risovič et al., 2008).

Then, with the component labeling, each voxel was assigned to a specific Look-Up Table and the image prepared for the analysis. Pore3D calculates the number of connected elements, their volume, sphericity (ratio of the surface area of an equivalent sphere to the surface area of the object), and diameter of the maximal inscribed sphere.

The Morphometric Analysis module performs a series of measurements derived from bone histomorphometry (Accardo et al., 2005) and here first applied to rock samples. In this sample the glassy matrix between pores is a few voxels thick and shows features that are well resolved in the tomographic image. A characterization of the morphology of the matrix network is hence straightforward and relevant. The analysis is based on the parallel plate model (Parfitt et al., 1983) and it introduces morphometric parameters that have become a standard “de facto” for quantitative analysis of matrix-supported materials. Although the applicability of the parallel plate model is not easy to assess, the calculated parameters provide a good insight into the sample microstructure, and they have been widely adopted in medical research, where a wide bibliography exists (see, for example, Odgaard, 1997) as well as in other sciences (for instance, for food science, see Falcone et al., 2005). An additional advantage of this approach is the availability of the software for the calculation, as most 3D visualization programs include this package. Structural parameters that can fully describe the network complexity, such as structure thickness, structure separation, and fragmentation index (also called trabecular pattern factor, a measurement of the convexity and/or concavity for the pattern connectivity, Hahn et al. [1992]), have been calculated for this specimen.

The sample AMB shows an overall structural simplicity given by the physical disconnection of most vesicles. However, voids statistics often need to consider fully separated objects, for example, to compute vesicle number density and vesicle size distributions (Proussevitch et al., 2007; Polacci et al., 2006; Polacci et al., 2008). For this reason, a series of algorithms have been implemented in Pore3D to separate voids and to reconstruct the vesicles as they were before the eruptive coalescence process connected them together. The watershed transformation is used to close the vesicles by simulating a “water flooding” of the image starting at the lowest gray values, which are computed on the distance-transformed image (Meyer, 1994; Borgefors, 1996).

Moreover, when determining the object morphological characteristics, objects intersecting the boundaries should not be included. Therefore, the code also performs a “modified” connected components labeling starting from the image margins, where intersecting objects are labeled with the background value, i.e., they are removed from the image. The connected components of sample AMB have been labeling-processed to verify the presence of void-objects of different dimensions (Fig. 4A). Later, a “modified” labeling process following a watershed transformation was able to show the abundance of few-voxels-sized objects not intersecting the volume borders (Fig. 4B).

The STR2 pumice sample (Fig. 2B) was chosen to investigate the connected void fraction topology by processing two subvolumes of 200 × 200 × 200 voxels each, equivalent to a cube of 1.8 mm side length. This sample is significant because it contains a complicated network of vesicles, which constitutes a challenging task for the topological characterization of the interconnected spaces. Moreover, the sample STR2 shows a macroscopic gradient of porosity throughout the pyroclast. The interior of the investigated sample (STR2a) has wide, expanded pores, strikingly in contrast with the part close to the outer surface (STR2b), where the porosity is lower and vesicles are of smaller size (Figs. 3B and 3C). This observation lets us approach the afore-mentioned problem of the REV extraction in strongly heterogeneous samples. Indeed, the same analysis performed on the inner and outer sectors of this sample produces unsurprisingly different results. For more simple cases, Pore3D implements its own algorithm to find the optimal “3D sampling window,” but visual inspection by the user remains valuable.

The connected void fraction has been characterized by considering its skeleton, which is a thinned 1D representation of 3D vesicles, with their pore spaces channels and throats (Lindquist and Lee, 1996), and can describe the connectivity of the voids in terms of pathway tortuosity. Since there is no precise definition of the skeleton of a 3D digital image (Cornea et al., 2007), different skeletonization algorithms (Fig. 1), such as the one proposed by Palágyi and Kuba (1999), the so-called “LKC” (Lee et al., 1994), the DOHT (Pudney, 1998), and a modified version of the Cornea et al. (2005) algorithm, have been implemented in Pore3D to let the user choose the most suitable one.

In this study, we applied the LKC algorithm because it is less sensitive to the “noise” at the pore boundaries that often generates redundant parts that may disturb the topology of the skeleton (Bai et al., 2007). Moreover, by visual inspection, the LKC skeletonization seems appropriate for the description of the texture of most of our samples. The numerical results are also in agreement with the Euler characteristic and branches/nodes ratio values. In some cases, the skeleton was extracted from preprocessed images, such as binary smoothed versions of the input data, and it was improved by postprocessing pruning. Pruning algorithms are used to make the skeleton simpler for the analysis by removing the unnecessary elements that fall under a threshold length. The skeleton obtained for sample STR2 is shown in Figure 5.

The Anisotropy Analysis module can investigate another parameter often used in material sciences to characterize the sample's structure. It corresponds to the preferential orientation of the objects of interest, in our case the vesicles, for which we can assume that the original spherical shape has been modified before magma quenching. Sample AGN (Fig. 2C) shows a preferential elongation of bigger pores parallel to the pyroclast major axis and it has been selected to verify how Pore3D processes anisotropic distributions of vesicles. A VOI of 300 × 300 × 570 voxels (2.0 × 2.0 × 3.8 mm³) is large enough to preserve this peculiarity (Fig. 3D), which is shown by the crystals as well, organized coherently with the vesicles’ anisotropy (Fig. 6). In the Anisotropy Analysis module, the code uses the mean intercept length (MIL) method (Whitehouse, 1974; Ketcham and Ryan, 2004), which is based on the average length of directional chords intercepting the objects of interest and defines an elongation index as the measure of the preferred fabric orientation. We have obtained the elongation index for the whole void space and then, by a simple algorithm that returns the volume containing the biggest connected component in the volume, we have iteratively performed the same morphometric analysis on each of a hundred biggest connected voids. For the void fraction, Pore3D also computes the isotropy index, which is a measure of the similarity of the fabric to a uniform distribution (Benn, 1994).

Finally, we have studied a multiphase system, sample STR1 (Fig. 2D), which presents at least two different classes of crystals embedded in the glassy matrix together with the voids. Pore3D implements several automatic multiple thresholding methods to segment gray levels belonging to different phases (as the “extended” version of the Otsu method or the K-means clustering multiphase segmentation). For sample STR1 we have adopted a multistep simple thresholding approach and we have segmented the crystalline fraction from the background, in a VOI of 314 × 314 × 204 voxels, equivalent to a volume of 2.8 × 2.8 × 1.8 mm³ (Fig. 3E). While the void fraction can be easily disregarded, the glassy matrix is partially segmented together with the objects of interest. To overcome such problems, Pore3D offers morphological operators such as erosion and dilation, through which voxels of the object boundaries can be removed or added (Serra, 1982). With these tools we have first removed spurious borders of the object and then restored the original edges, so as not to loose any significant volume. Moreover, the use of phase-contrast μCT has been particularly effective for this sample. In fact, due to the contrast enhancement at their boundaries, we were able to visualize objects in the sample, featuring very small variations in density and chemical composition. In the processed and enhanced image, different crystalline phases were easily recognized (most probably feldspars and pyroxenes) and separated from vesicles and glassy matrix to quantify their volumetric parameters and morphometric properties (Fig. 7). The pyroxene volume crystallinity is ∼28% and the feldspar crystallinity is around 12% of the VOI.

In the study of volcanic materials, the textural information obtained with X-ray μCT represents a pioneering and attractive area of research. In this section we discuss in more detail the measured parameters, with reference to their possible volcanological implications. Significant information can be extracted and related to the volcanology of the specimens, making the textural data a complementary tool in the study of the eruptive dynamics (see, for example, Andronico et al., 2008); however, this application is well beyond the objective of this paper.

The main steps of the quantitative analysis that we introduced in the previous section have been applied to the whole set of volcanic products. In Table 2 we list the analysis results for each sample in order to compare them and make their significance clear in the volcanological context. The first parameter in Table 2 is the porosity, described as the void-to-total volume fraction. It is high in the frothy pumice, especially in its inner part, being lower in the scoria samples and in the tube pumice. The porosity variability through the sample types supports the idea that this parameter is more related to the eruptive dynamics than to magma composition, which is almost identical for the samples from Stromboli volcano. In the case of the highly vesiculated Stromboli pumice, we ascribe the gradient in the sample porosity to further volatile exsolution after magma fragmentation. In the case of the tube pumice sample AGN, the porosity is close to that of the scoria samples, possibly a consequence of vesicle elongation that may have affected the ability of vesicles to nucleate, grow, or coalesce (Wright et al., 2009; Bouvet de Maisonneuve et al., 2009).

Of the other Minkowski functionals, calculated on the void objects, the integral of the mean curvature spans a wide range of values, from positive to negative, indicating a dominance of convex or concave surface shapes, respectively (Russ and DeHoff, 1986; Velichko et al., 2008). For example, the value for sample STR2 is very low and this points to very rough surfaces, with a large number of coalesced pores and concave surface elements. The other samples have positive values of the integral of the mean curvature, i.e., a high fraction of convex surfaces, as occurs for samples dominated by isolated spherical voids. The Euler characteristic has the attributes of a topological order parameter describing the spatial connectivity: configurations with a positive value typically consist of isolated objects dispersed in the host component, while negative values indicate connected aggregates (Mecke, 1996; Odgaard and Gundersen, 1993; Velichko et al., 2008). In our data, for example, the lowest values are found for sample STR2, which shows a high connectivity of pores, and, as confirmed by the connected component analysis (see later), the lowest number of vesicles.

With the Connected Component Analysis, complete information about the shape and size of the measured objects can be obtained (not listed in Table 2). Parameters like volume, diameter of the maximal inscribed sphere, diameter of the equivalent sphere, and sphericity can be calculated for each of them. The knowledge of the shape and surface area of the components is especially desirable in particle systems, where they are important characteristics for the understanding of surface processes such as adsorption and dissolution.

After isolation of individual pores, several measurements can be performed on each vesicle in order to obtain statistics about the overall vesicle size distribution. For sample AMB the component analysis has been performed again after the vesicle separation and the suppression of the vesicles connected to the image margins. The vesicle volume distribution (Proussevitch et al., 2007) can be fit by a power law and the exponent governing the size distribution can be associated with specific vesicle nucleation, diffusion, and decompression processes (Polacci et al., 2009a; Gaonac'h et al., 2005; Blower et al., 2003; Blower et al., 2001; Gaonac'h et al., 1996a; Gaonac'h et al., 1996b). Moreover, knowledge of the vesicle volume distribution can be critical for the REV selection, when the bulk vesicularity needs to be physically measured. The knowledge of the rate at which the frequency of large vesicles decreases with increasing size contributes to the definition of the minimum representative sample size.

The abundance of more or less isolated vesicles can be also evaluated on the basis of the morphometric analysis of the glassy walls. For pumices with highly interconnected voids, as sample STR2, the structure thickness is low and septa are densely and randomly arranged, as indicated by the structure linear density, which is high for sample STR2 compared to the other samples. In addition, the trabecular pattern factor is used to quantify the connectivity of the matrix septa and hence to summarize the material microarchitecture. This factor calculates the relation between convex and concave surfaces of the matrix frame and describes a well-connected matrix as a prevalence of concave surfaces, with a low factor value. On the other hand, a prevalence of convex surfaces points to a matrix frame with few interconnections among surfaces and gives a higher value of the trabecular pattern factor. Among our samples, this value is low for the scoriae STR1 and AMB and for the tube pumice AGN, where isolated vesicles and concave surfaces of the matrix prevail, being higher in the porous pumice STR2, characterized by a scarcely connected frame and by higher structure separation (Table 2). In fact, the structure separation parameter measures the separation of the solid septa and gives an indication of the average linear pore dimensions.

The skeleton analysis allows us to make further comparisons between our samples and to investigate the complexity of the void fraction. This information can be used as an additional indication of the permeability of the material, as it reflects, together with the pore shape, the efficiency of percolation paths for fluids (Wright et al., 2009). Once skeleton joints (nodes in Table 2) and paths (branches) are calculated, the voids which have their center at the nodes are processed and merged in a single pore, if they overlap each other. In Table 2, we observe that the sample STR2b shows the largest number of pores. However, the most outstanding feature of sample STR2 is the large number of branches, which, together with the number of pores, give remarkably large coordination numbers (i.e., number of branches connected to each pore-body). This observation indicates a high skeleton complexity and a high connectivity, which are consistent with a Euler characteristic much less than 0 and a low number of connected components.

Finally, we focus on the elongation and isotropy indices as parameters often used to better model the pore structure and its relationships with the permeability and the topological effects of the deformation (Nakamura et al., 2008; Wright et al., 2009). The elongation index measures the preferred orientation of a fabric in the plane containing the maximal axes of the deformation ellipsoid (Harrigan and Mann, 1984), and varies between 0 (no preferred orientation on such plane, but with an overall shape ranging from isotropic to pancake-like shapes) and 1 (perfect preferred single orientation). The isotropy index measures the similarity of a fabric to a uniform distribution and varies between 0 (all observations confined to a single plane or axis) and 1 (perfect isotropy).

The anisotropic pattern recognized in the tube pumice sample AGN is particularly intriguing. In fact, the overall elongation of the macroscopic porosity is confirmed by the elongation index for the voids, which is much higher than in the other samples (0.33 compared to values <0.1). The isotropy index points to a similarity of the “deformation” directions in the plane normal to the elongation axis. Vesicle deformation is attributed to stresses during magma ascent in the conduit, which deforms the voids from their spherical shape (Rust et al., 2003). Recent experimental studies have demonstrated that deformation increases the rate of vesicle coalescence in magma, so producing a high void connectivity (Okumura et al., 2006; Nakamura et al., 2008; Okumura et al., 2008). In the tube pumice sample AGN, the relatively high Euler characteristic and the total number of connected components indicate an abundance of isolated pores, which seems to contradict this observation. As the basic parameters are statistical measures, we think that their overall values are affected by the quantity of smaller vesicles (connected components number), which are mainly isolated, despite the volumetric dominance of elongated, bigger vesicles. In fact, if the anisotropy analysis is carried out only on the 100 larger vesicles, the elongation index is higher than for the whole set of components, with an average value of 0.37 against the general one of 0.33. The low sample porosity could be related to the deformation and coalescence of the largest vesicles, which create the pathways for gas flow (Polacci et al., 2006; Polacci et al., 2009b). This would reduce the porosity and create collapsed vesicles textures (Mongrain et al., 2008), only partially reconstructed by a second-phase nucleation of smaller vesicles.

In the study of geological materials, the morphological and textural information obtained from a nondestructive imaging technique represents a promising field of research. X-ray μCT, combined with 3D quantitative analysis, provides an extensive microstructural analysis tool and supports the investigation of parameters of great interest in volcanology (see Degruyter et al., 2010, for extended bibliographic references). In particular, the analysis of vesicle size, shape, distribution, orientation, and degree of interconnectivity, quantifies aspects that are directly related to the magma nature and dynamics (Polacci et al., 2009a).

In this paper we present several procedures to extract quantitative parameters from μCT images of volcanic rocks. These analyses have been performed by using the Pore3D software library, custom-developed at Elettra in the last few years and especially devoted to the study of porous media imaged with X-ray μCT techniques.

Pore3D can perform filtering, segmentation, skeletonization, and quantitative analysis of 3D data with a flexible approach, by choosing and implementing the appropriate set of state of-the-art functions and procedures. As there is not a unique methodology to extract quantitative information from 3D data sets, and porous media with vesicle-like pores (volcanic rocks, polymeric foams, or food products) are very different from fibrous materials (felts or natural wood), trabecular structures (trabecular bone) or randomly distributed materials (reservoir sandstones and rocks), Pore3D has been designed to take into account this heterogeneity. For this reason it offers different alternative approaches and it can be easily extended to the study of other porous materials.

The application of Pore3D to volcanic clasts allowed us to identify the distinctive features of each type of pyroclast we analyzed. Highly vesiculated (frothy) pumices show vesicles arranged in an isotropic pattern of coalescing aggregates, especially in the interior of the sample. The scoriae show lower porosity and an abundance of isolated vesicles. In the more crystalline scoria sample, most of the vesicle separation is due to the presence of crystals of different phases and morphologies. Finally, the tube pumice sample confirms at the microscale the anisotropic distribution of its constituents, observed also at the scale of the hand sample. Larger, elongated vesicles are seen together with a generation of small, mainly spherical and isolated vesicles.

The recent development of Pore3D has been aimed at averaging the textural parameters and extracting a single value to describe the general spatial distribution of some specific features. As further improvement of the software library, we envisage the visualization of a part of the local parameters in their spatial distribution, on “layer” maps, as done by the GIS (Geographic Information System) software (Barraud, 2006). A better knowledge of the geomaterials can be achieved if the local property variations throughout the sample are known. Moreover, the modeling of parameters like tortuosity and permeability represents an interesting field of research in the study of porous systems as volcanic rocks. The implementation in the Pore3D code of algorithms calculating these parameters is in progress.

The understanding of the rock physical properties can be greatly improved by the recognition of object populations (vesicles, grains, crystals, glass heterogeneities, fractures, etc.) with various morphological characteristics and by mapping them in order to match different textural information. The knowledge of textural parameters such as the porosity and vesicle anisotropy, the matrix wall roughness, and the skeleton architecture, with their spatial distribution, can greatly improve the study and quantification of flow characteristics.

However, a comprehensive understanding of critical volcanic processes such as vesiculation, development of magma permeability, fragmentation, and ascent toward the Earth's surface can be reached only by combining the textural information with other geochemical and geophysical results.

We thank the project ANR-CNRS VOLGASPEC-ANRCTT2006 “Tracking Volcanic Eruptions from the Dynamics of Magma Degassing processes”; project PRIN MIUR “Thermodynamics and kinetics of magmatic degassing: observations, experimentation and modeling”; and project INGV-DPC “Paroxysm” under which part of this research, including sample collection, was performed. We also thank the Regione Friuli Venezia Giulia project “Development of an X-ray portable system for non-destructive analysis of archeological and artistic materials” and Dr. Stefano Favretto for his assistance during the tomographic measurements.