Textures (sizes, shapes, distributions) of glass inclusions in crystals can provide important insight into magmatic processes. Three-dimensional X-ray tomography is an excellent nondestructive method for assessing igneous rock textures in general; however, standard absorption-contrast tomography is not useful for imaging glass inclusions because it typically does not provide enough contrast between the inclusion and crystal host to allow quantitative separation of these phases during image processing. Propagation phase-contrast tomography enhances contacts between different phases, allowing the use of edge-detection algorithms to separate glass inclusions from their host crystal. The sample-to-detector distance is increased to achieve edge enhancement at the cost of some loss in image resolution. We qualitatively assessed images acquired at a series of distances (30, 100, 150, 200 mm) to determine the best distance for imaging given this trade-off.

Typical image processing approaches used for absorption-contrast tomography are not useful for edge-detection problems, so we developed an IDL-based graphical user interface (GUI) to process propagation phase-contrast tomograms. A combination of grayscale filters, distance thresholding, erosion processes, and individual picking of outlier voxels is used to fit a convex hull to an individual inclusion, and the output provides information on its size and shape. The processed inclusions can be placed in the context of the original host, such that size, shape, and location of multiple inclusions can be assessed individually or in combination.

We apply this method to five glass inclusions in quartz crystals from early erupted Bishop Tuff. Results indicate that inclusions range in shape from nearly spherical to nearly ellipsoidal and from round to faceted. The relationship between degree of faceting, size, and location within the host is consistent with melt inclusion faceting under magmatic conditions over centennial to millennial time scales. We also show application of this method to other kinds of studies, such as glass inclusion size distributions.


Igneous rock textures have long been used to qualitatively assess the evolution of magmatic systems, and quantitative analysis of these textures has been a growing area of study in recent decades (Hersum and Marsh, 2007). Quantitative analysis has traditionally been limited to two-dimensional (2D) data sources that require application of stereological corrections (e.g., Sahagian and Proussevitch, 1998; Higgins, 2000; Jutzeler et al., 2012), or to crystal separates, which eradicate the context of an individual crystal and are prone to bias introduced in laboratory processing (e.g., crystal fragmentation). X-ray tomography is an increasingly important tool in both quantitatively and qualitatively assessing textures in three dimensions in rocks and other materials (e.g., Carlson et al., 2000; Mock and Jerram, 2005; Jerram and Higgins, 2007; Pamukcu et al., 2012; Gualda et al., 2012a) because it is a relatively nondestructive technique that allows for in situ analysis of features and avoidance of stereological corrections.

In parts 1–3 of this series, we outlined methods and applications of conventional and differential absorption X-ray tomography in the study of igneous rock textures (Gualda and Rivers, 2006; Pamukcu and Gualda 2010; Gualda et al., 2010). We have used these methods to qualitatively assess the spatial distributions of major and accessory minerals in pumice clasts (Gualda et al., 2010; Pamukcu et al., 2013) and to quantitatively determine crystal size distributions of major and accessory minerals to probe time scales of crystallization and the onset of eruptive decompression in giant magma bodies (Pamukcu and Gualda, 2010; Pamukcu et al., 2012, 2013; Gualda et al., 2012a). Here we focus on X-ray tomography as a tool to study textures of glass inclusions in crystals, the sizes, shapes, and distributions of which can have important petrologic implications, such as providing information on crystallization time scales (Skirius et al., 1990; Manley, 1996; Gualda et al., 2012a). Glass inclusion textures have been studied (e.g., Lowenstern, 1995; Manley, 1996; Bindeman, 2005; Gualda et al., 2012a) through observations using the optical microscope on crystals placed in refractive index oil; however, estimates of size, shape, and position of inclusions are limited by the nature of the method. We describe here a method that yields more accurate and precise measurements of inclusion size, shape, and position, and it also allows three dimensional (3D) visualization of inclusions.

The difference in X-ray attenuation between a glass inclusion and its crystal host is often small, such that they can be difficult to impossible to distinguish in tomograms collected using standard absorption-contrast methods. Similarly, differential absorption-contrast tomography cannot be used to isolate inclusions because of low concentrations of the high atomic number elements that can be analyzed by this method (for details see Gualda et al., 2010). Given these limitations, we have explored a different method to distinguish and quantify glass inclusion textures, i.e., propagation phase-contrast X-ray tomography. This method enhances edge contrast of features, which makes separation of features with similar or same linear attenuation coefficient much more feasible than with standard tomography. We show here that this technique can be successfully used to quantitatively characterize the textures of rhyolitic glass inclusions in quartz. We describe the protocols used for imaging and processing, including a program we developed for image processing, and we show a few examples of inclusions imaged in quartz from the Bishop Tuff (California, USA). Our treatment is currently restricted to convex objects, so we focus on fully enclosed inclusions.


X-ray tomography is a 3D imaging technique that yields 3D maps (or tomograms) of specific sample properties. Visualization can be achieved by converting the value in each volume element (voxel) into grayscale intensity, which also allows for use of image processing techniques. Most commonly, X-ray tomography is based on absorption contrast, where the measured property is the linear attenuation coefficient, a function of the X-ray energy, phase density, and elemental composition (see Gualda and Rivers, 2006). During image processing different phases can be quantitatively distinguished using grayscale filters (see Ketcham, 2005); however, distinction of different phases is hampered or made impossible when the difference in attenuation between phases of interest is small. This is often a problem in studies of volcanic rocks, where the linear attenuation coefficient of glass or a mineral phase is similar to or the same as that of another mineral phase (e.g., quartz, feldspar). Similarly, glass inclusions in a crystal can often be qualitatively distinguished from the surrounding host, but the attenuation contrast between the inclusion and host may not be sufficient to allow quantification. It is sometimes possible to overcome these issues by using differential absorption contrast, which yields 3D maps of element concentrations, but unfortunately imaging energy requirements are such that the method is mostly limited to heavy elements such as Zr and rare earth elements (see Gualda et al., 2010).

To overcome these limitations, we use propagation phase-contrast X-ray tomography to enhance edge contrast between crystals and their glass inclusions. In propagation phase-contrast tomography, the contrast is generated by refraction of X-rays at the interface between materials with different refractive indices. X-rays incident on the interface are bent toward the material with the higher refractive index, which will be air in the case of a solid-air interface. This is analogous to the Becke line phenomenon in visible light microscopy. The refraction angle is quite small, <1 mrad, so propagation distances on the order of 0.1 m are required for the Fresnel fringes to be displaced by more than a pixel on the detector (Snigirev et al., 1995; Cloetens et al., 1997; Tafforeau et al., 2006).


We performed all tomographic imaging on the bending magnet beamline of the GeoSoilEnvironCARS (GSECARS) sector of the Advanced Photon Source at Argonne National Laboratory (Chicago, Illinois). Details of the experimental setup of the GSECARS beamline can be found in Rivers et al. (1999), Sutton et al. (2002), and Gualda and Rivers (2006). Individual crystals containing glass inclusions were mounted on metal pins and placed on a rotating stage (Fig. 1). Beam energy was 20 keV for all analyses. During imaging the sample was rotated a total of 180°, and a radiograph was recorded by a charge-coupled device (CCD) camera at each step, for a total of 720 (2 × 2 binning in CCD camera) or 900 (no binning) radiographs collected for each tomogram. The resulting image resolution is ∼2.8 (binned) and ∼1.4 (unbinned) μm/voxel, respectively. The tomo_display program (Rivers and Gualda, 2009) was used to reconstruct these radiographs to obtain a 3D map of linear attenuation coefficient (Rivers et al., 1999).

To achieve enhanced edge contrast by propagation phase-contrast, the sample-to-detector distance, kept to the minimum possible in absorption-contrast tomography (∼30 mm in the GSECARS setup), was increased to as much as ∼200 mm. Increasing this distance has the adverse effect of reducing image resolution; therefore, we tested four sample-to-detector distances (30 mm, 100 mm, 150 mm, 200 mm) in order to assess the trade-off in edge contrast and image resolution. Tomograms show an increase in edge contrast with increasing distance (Fig. 2), as clearly shown by tomograms filtered using edge-detection algorithms (Fig. 3). The 200 mm distance was the farthest distance achievable without requiring significant physical changes to the setup. Qualitative assessment of edge-detected images at the various sample-to-detector distances shows that the largest distance (200 mm) resulted in the greatest edge enhancement with comparatively little loss in resolution, and the results were sufficient for our image processing procedures to be successful. Thus, we imaged all remaining crystals with a sample-to-detector distance of 200 mm.

We further compared the results of imaging performed using 2 × 2 binning of pixels in the CCD camera versus no binning; 2 × 2 binning results in a factor of 4–5 reduction in acquisition time (from ∼1 h to <15 min per crystal), at the cost of a factor of 2 reduction in resolution (from 1.4 to 2.8 µm per voxel). Comparison of binned and unbinned images shows some qualitative difference in edge detection (Fig. 4). A decrease in image resolution was noticeable in the binned images, but this decrease did not hamper our image processing. We considered this decrease in image quality from binning to be small compared with the significant decrease in imaging time that resulted.

Image Processing

In our previous three papers on this topic, we used the Blob3D program (Ketcham, 2005) for much of the image processing; however, Blob3D has more limited use in edge-detection problems, and its output does not include a number of relevant quantities for our work with glass inclusions. Consequently, we have developed an IDL-based GUI program, vol_inclusion_hull, for processing these tomograms that will be made available free of charge (see following for details).

We first apply an edge-detection filter to the whole tomogram; a Sobel filter yields the best results among several options we tried (the program vol_sobel can be used to apply a Sobel filter to a tomogram in 3D). Single glass inclusions are clipped from both the original and the edge-detected tomograms using the vol_trim program (Figs. 5A, 5B; Rivers and Gualda, 2009). Simple grayscale thresholding is used within vol_inclusion_hull to select the voxels describing the inclusion edge (Fig. 5C). In order to construct the surface of the inclusion, we determine a convex hull around the selected voxels, which yields the smallest convex 3D envelope that contains all the selected points (Figs. 5C, 5D); because a convex hull is used, only convex inclusions can be properly processed with our current approach. We adjust the fit of our convex hull interactively in vol_inclusion_hull using one or more of the following methods: (1) adjusting the grayscale threshold values to select brighter or darker pixels; (2) applying a distance filter such that the convex hull is created using only voxels that are within a designated distance from the centroid of the object; (3) applying an erosion filter to the selected voxels to reduce noise within and thickness of the boundary region; (4) manually selecting individual voxels for removal.

Outliers are the main source of difficulty in processing these volumes (Fig. 6). We use three methods to monitor for them and to assess the quality of the fit of the convex hull.

  1. We display 2D slices through the original volume overlain with a polygon representing the 2D trace of the convex hull for that slice (Figs. 6A, 6D). This allows us to see the location of individual outliers and remove them.

  2. We create a histogram describing the distance of selected voxels from the object centroid (Figs. 6B, 6E). This allows us to determine the distance at which the majority of the selected voxels are located; voxels at distances substantially farther than typical are considered potential outliers.

  3. We create a 3D image that displays all of the points selected by the grayscale filter and the polyhedron created based on the current values for the distance and erosion parameters (Figs. 6C, 6F).

These operations can all be performed with vol_inclusion_hull, and parameters can be easily changed with results displayed in real time. This allows for quick processing of individual inclusions. We find that most inclusions can be processed in ∼5–10 min. Output of the program includes (1) the parameters used to create the hull, including a list of manually excluded voxels that can be imported into the program and used for fast reprocessing if needed; (2) basic metrics characterizing the inclusion, including volume, surface area, sphere-normalized surface area-to-volume ratio (SNSVR; see Ketcham, 2005), best-fit ellipsoid properties (minimum, intermediate, and major axes), properties of the best-fit ellipsoid scaled to have the same volume as the inclusion, and ellipsoid-normalized surface area-to-volume ratio (ENSVR); and (3) a file containing IDL variables used and created during processing that can be read into IDL and used to create new visualizations or perform additional calculations.

The method described here is not only useful for characterizing inclusions in crystals but any objects that can have their edges emphasized by edge-detection filters. We use our procedure to characterize and create visualizations of the host crystal (see following).

Multiple inclusions can be combined into a single visualization using the program vol_inclusion_gather, which uses the output of vol_inclusion_hull and known locations (recorded as output by vol_trim) of each subvolume (i.e., volumes containing individual clipped inclusions) within the original volume.

The location of an inclusion relative to the nearest edge of the crystal host can be determined using the program vol_inclusion_distance. This program first repositions the inclusion within the crystal host as in vol_inclusion_gather; it then determines the distance between the inclusion centroid to all the triangles that make up the convex hull of the crystal, and the smallest calculated distance is output as the distance from the crystal edge.

The programs vol_inclusion_hull, vol_inclusion_gather, and vol_inclusion_distance, as well as other tools assembled under the name “vol_tools” (Rivers and Gualda, 2009) can be obtained from us or from http://my.vanderbilt.edu/ggualda/vol_tools/. All of these programs can be run free of charge using the IDL Virtual Machine (see www.exelisvis.com/ProductsServices/IDL/DevelopmentEnvironment.aspx).


Faceting Time Scales

Our motivation for developing this method was to assess glass inclusion size, shape, and location for determination of the longevity of giant magma bodies (see Gualda et al., 2012a). A full discussion of the implications of our results is beyond the scope of this manuscript and will be provided in subsequent contributions. The intent of this section is to show the kinds of data that can be derived using the method described here.

We focus on fully enclosed glass inclusions in quartz crystals from the early erupted Bishop Tuff (California), using some of the same samples studied in Gualda et al. (2012a). We show inclusions from two quartz crystals (F815_xtl827, F815_xtl831) selected from a single pumice clast from fall unit F8 (following the classification of Wilson and Hildreth, 1997) of the Bishop Tuff. Pumice clasts were first crushed lightly and crystals containing clear, nondevitrified glass inclusions were picked and mounted on metal pins with epoxy (Fig. 1). Effort was made to ensure that no part of the crystal was embedded in the epoxy. Imaging and processing followed the protocols described herein.

Studied inclusions show a wide range of shapes, from nearly spherical to more ellipsoid, and also vary substantially in the degree of faceting, as revealed by the deviations of SNSVR or ENSVR from 1 (Fig. 7). Considering only the four inclusions from sample F815_xtl827 (Fig. 8), we find that the largest inclusion (inclusion 1) is the most faceted and located nearest to the center of the crystal (0.422 mm from crystal edge). The smallest inclusions (inclusions 3 and 4) are the least faceted and nearer to the edge of the crystal (0.200 mm and 0.299 mm from crystal edge, respectively), while the intermediate-sized inclusion (inclusion 2) is faceted and positioned between the largest and smallest inclusions (0.304 mm from crystal edge). These observations are consistent with the idea that melt inclusions gradually facet after entrapment during crystal growth under magmatic conditions (Skirius et al., 1990; Manley, 1996; Gualda et al., 2012a).

These results may initially seem counterintuitive, as smaller inclusions are expected to be more faceted than large inclusions due to the smaller volume that needs to be diffused. However, position is also an important variable in the problem: given inclusions of the same volume, those closer to the center of the crystal are expected to be more faceted than those nearer to the edge. Considering these two factors, size and position, in combination, our finding of large inclusions that are more faceted than small inclusions is reasonable and does not conflict with the approach of using faceting to obtain estimates of residence time.

Using the kinetic expression provided in Gualda et al. (2012a), we calculate the time required for full faceting, i.e., the time required for the inclusion to obtain the full negative crystal shape (e.g., a hexagonal bipyramid, in the case of quartz), of these five inclusions. This faceting time is a function of temperature (estimated between 720 and 780 °C for the Bishop Tuff; Hildreth, 1979; Bindeman and Valley, 2002; we use 750 °C; Gualda et al., 2012b; Gualda and Ghiorso, 2013), diffusion of silica in melt (1.28 × 10–14 m2/s; Baker, 1991), a series of constants, and the change in volume. It is important that the faceting process is modeled assuming volume invariance, i.e., there is no net gain or loss of volume of the inclusion. In this sense, faceting is attained by diffusion of material from curved regions, where dissolution takes place, toward corner regions, where reprecipitation occurs (Gualda et al., 2012a). It is this volume, i.e., the volume that has moved in the faceting process, that we call the “change in volume”; to calculate the change in volume, we compare the volume of the inclusion determined with tomography to that of a hexagonal bipyramid of the same volume. The errors from this calculation were thoroughly examined (Gualda et al. (2012a), and the error was estimated as ∼125%.

All five inclusions give full faceting times in the range of 102–103 a (Figs. 7 and 8), consistent with those obtained in Gualda et al. (2012a). Considering only the inclusions in sample F815_xtl827, we find that the largest, most faceted, most centrally located inclusion would require the longest time for full faceting (inclusion 1: 674–2359 a), while the smallest, less faceted, more edgeward inclusions would require the shortest times for full faceting (inclusion 3: 164–573 a, inclusion 4: 222–777 a). The inclusion that is intermediate in size, position, and faceting gives times for full faceting intermediate between the largest and smallest inclusions (inclusion 2: 422–1476 a). Given that none of the inclusions considered is fully faceted, these results are consistent with centennial to millennial crystallization times for quartz from the Bishop Tuff, in agreement with the findings in Gualda et al. (2012a) and Pamukcu et al. (2012).

We are currently working on comparing time scales obtained from glass inclusion faceting with those from diffusion profiles along boundaries containing inclusions. Toward this end we are imaging crystals with tomography, subsequently sectioning and polishing them to reveal inclusions of interest, and imaging them in 2D using a cathodoluminescence detector attached to a scanning electron microscope to reveal zoning. Preliminary results suggest that faceting times compare very well with diffusional relaxation times; presentation of the results will be the subject of an upcoming contribution.

Glass Inclusion Size Distributions

Bindeman (2005) suggested that glass inclusion size distributions and the radius ratio of a glass inclusion to the crystal host (Ri/Rh) can provide information on inclusion capture and crystal fragmentation. Such size distributions and quantities can be readily determined using the methods we describe here. To demonstrate this, we determine the Ri/Rh for the four inclusions from crystal F815_xtl827 using tomography data (Fig. 8).

For this application, we calculate the equivalent radius (i.e., the radius of a sphere with the same volume) of inclusions and the host crystal; the host crystal is processed in the same manner as the inclusions. Embayments in the crystal result in an overestimate of the volume, but it is clear from Figure 8 that the crystal host is also not a whole crystal. We expect both factors to partially offset each other. Further, even if the crystal volume we determine is larger or smaller than the true volume by as much as a factor of two, the equivalent radius we obtain is only overestimated or underestimated by 20%. As a result, we consider this to be an adequate estimate of the crystal size.

The Ri/Rh of the four inclusions in crystal F815_xtl827 ranges from 0.05 to 0.08. These results are consistent with those of Bindeman (2005), who found that 90% of glass inclusions in quartz crystals from late erupted Bishop Tuff and the Lower Bandelier Tuff have small Ri/Rh (0.02–0.08), suggesting that inclusions with larger Ri/Rh are less likely to survive in the face of even small overpressures. The data we present here are too few to draw similar conclusions on inclusion capture and crystal fragmentation, but the method is clearly useful and appropriate for such studies.

Further Work

At this juncture, our processing routine is only appropriate for use on single, fully enclosed inclusions due to complications in the following.

  1. Creating convex hulls of multiple objects in a single image: currently, when multiple inclusions are processed together (i.e., not clipped and processed individually) they are combined into a single convex hull. For the purposes of this work, this is not a prohibitive problem as we do not require data on every inclusion in a given crystal, and we can process single inclusions in tens of minutes.

  2. Appropriately describing concave objects: the convex hull operation is only appropriate for inclusions that do not have concave embayments. In the case of an embayed or concave object, the convex hull operation overestimates the volume of the object because the embayment is effectively ignored.

Work to devise methods to process multiple inclusions at once and objects with concave shapes is ongoing.

In our work thus far, we have focused on the characterization of glass inclusions in quartz crystals. Future work will include imaging and quantification of textures of glass and/or mineral inclusions in other minerals (e.g., feldspar, pyroxene, olivine, amphibole, and other accessory minerals). Beyond the characterization of the texture of inclusions, the ability to locate inclusions within crystals in 3D may prove helpful in efforts to expose specific inclusions through sectioning of the host crystals (e.g., Fourier transform infrared spectroscopy sample preparation).


X-ray tomography is an excellent method for quantitatively and qualitatively studying textures of igneous rocks, and here we demonstrate its application in assessing glass inclusion textures. Processing images obtained using standard absorption-contrast tomography is difficult to impossible when features have the same or similar linear attenuation coefficient, as is often the case with glass inclusions in a quartz host. We demonstrate that propagation phase-contrast tomography eases these processing problems by enhancing edges of these features such that edge-detection algorithms can be successfully applied to quantitatively describe glass inclusion shapes. Edge enhancement is achieved by increasing the sample-to-detector distance; while some reduction in image resolution is expected, our results suggest that it is minor compared to the increased ability to enhance edges.

The typical approach for image processing that we have used in previous work, using the Blob3D program (Ketcham, 2005), is not suitable for edge-detection problems. As a result, we have developed an IDL-based GUI program for processing propagation phase-contrast tomograms of glass inclusions that uses a combination of grayscale filters, distance thresholding, erosion filtering, and individual selection of voxels to fit a convex hull to the inclusion of interest. This processing method is applicable to single, concave inclusions, yielding quantitative estimates of size, shape, and position of inclusions within the host crystal. Additional work is needed to process multiple inclusions at once and to characterize nonconcave inclusions. Nonetheless, there are diverse applications that the method presented here can be used for, such as quantifying textures of glass and mineral inclusions in other silicate minerals, and determining positions of inclusions in crystals for subsequent sectioning and analysis.

We demonstrate the application and success of this method in assessing glass inclusion faceting of five inclusions in quartz crystals from early erupted Bishop Tuff. Results show that glass inclusions have a wide range of shapes, from nearly spherical to nearly ellipsoidal, from round to faceted, and that their textures suggest quartz residence times on centennial to millennial time scales, consistent with results in Gualda et al. (2012a). We also show that this method can be used for studies of glass inclusion size distributions, akin to work done by Bindeman (2005).

Funding for this work was provided by a Vanderbilt University Discovery Grant and National Science Foundation grants EAR-1151337 and EAR-0948528 to Gualda. Portions of this work were performed at GeoSoilEnviroCARS (sector 13), Advanced Photon Source (APS), Argonne National Laboratory. GeoSoilEnviroCARS is supported by the National Science Foundation–Earth Sciences (grant EAR-1128799) and Department of Energy–Geosciences (DE-FG02-94ER14466). Use of the Advanced Photon Source was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences contract DE-AC02-06CH11357.