# American Mineralogist

- © 2000 American Mineralogist

## Abstract

Studies of crystal size distributions (CSD) can reveal much about how rocks solidify and under what conditions. Data from two-dimensional sections can be readily acquired at many different scales, from electron microscope images, thin sections, slabs, outcrops, and so on, but the conversion to true, three-dimensional values is complex. The widely used Wager method does not have a good theoretical basis and does not give accurate results. A modification of the Saltikov correction method is proposed here that is more accurate and can account for different crystal shapes and fabrics. Population densities determined by this method differ by factors of 0.02 to 100 from those determined by the Wager method. Published CSDs determined using other methods can be recalculated if the crystal shape and fabric parameters can be estimated. The method has been incorporated into a new program, CSDCorrections.

## Introduction

Textural studies in igneous petrology flourished in the nineteenth century, following the development of the petrographic microscope. At that time, textural observations were qualitative—average grain-sizes, grain relationship, and fabrics. Although qualitative data has its uses, it cannot constrain physical models of processes in the same way that quantitative data can. Textures are important because they result from the solidification of rocks, and hence if we want to understand solidification of igneous rocks then we must quantify and model texture. I use the word solidification rather than crystallization because solution of crystals may be important in many plutonic and even volcanic rocks (e.g., Higgins 1998; Marsh 1998).

There are many aspects of the texture of igneous rocks that can be quantified, but the most commonly quantified parameter is crystal size. The distribution of these sizes in three dimensions, the crystal size distribution (CSD), can give much more information on petrological problems than the mean, modal, or maximum size (e.g., Cashman 1990; Cashman and Marsh 1988; Marsh 1988). CSDs can reveal aspects of the thermal history of a magma and may give information on growth rates and magma repose times (e.g., Cashman 1993; Resmini and Marsh 1995).

There have been several studies of crystal sizes in igneous rocks during the last 40 years (one of the earliest was that of Jackson 1961), but the subject only started to mature in 1988 after the publication of the seminal papers of Marsh (1988) and Cashman and Marsh (1988). These authors established a theoretical basis for CSDs by the application of the industrial models of Randolf and Larson (1971). They also established the “CSD diagram” in which the natural logarithm of the population density is plotted against crystal size. Since then the subject has developed fast. However, some workers may have been reluctant to embark upon research in this field because there is no published general guide on how to determine CSDs. In this paper I try to fill that gap, based on my own experiences. In the first part I describe practical methods for acquiring data from two-dimensional sections, such as outcrops, slabs, and thin sections. The second part will concern the conversion of two-dimensional data to three-dimensional CSDs. I discuss aspects of the general problem and approaches to the solution. A computer program is presented that enables the conversion to be made simply.

## Methods of acquisition of crystal size distributions

Crystal size distributions can be measured directly when crystals can be extracted quantitatively whole from a volume of rock—a rather unusual situation for most igneous rocks, except perhaps carbonatites. Success also has been achieved using X-ray tomography of rocks for minerals at low concentrations, such as garnet or diamond (Denison and Carlson 1997; Rowe et al. 1997). In general, this method cannot separate touching crystals of the same mineral, although there has recently been some progress in this field (Proussevitch and Sahagian 2000).

Most CSD data are determined from two-dimensional sections through rocks—outcrops, slabs, or thin sections. The conversion from two-dimensional to three-dimensional parameters is not simple and will be discussed later. These two-dimensional sections are images that need to be processed to extract various parameters describing the intersections, such as length, width, area, perimeter, orientation, and centroid location. Processing can be either manual or automatic. In manual treatment, the different intersections are identified by eye, using color, birefringence, twinning, cleavage, and other properties. This technique is laborious, but can give high-quality data. Automatic image processing is much quicker, but is currently much less elaborate, using only colors for classification of particles. However, in some circumstances it can give very useful results. In these cases strict quality control must be enforced, for instance by manual analysis of a few selected samples.

The simplest textural analysis involves only one section. If there is an appreciable fabric (foliation or lineation), then sections parallel and normal to the dominant structure are necessary or a quantitative estimate of the nature and quality of the fabric, and its orientation with respect to the section. The maximum information on the texture of a rock can only be obtained from many serial sections. However, the amount of work involved means that crystal numbers, and hence statistical reliability, may be sacrificed. If only the mean dimensions of the crystals are needed and the rock is massive (no fabric), then several direct stereological techniques can be used. Most use two parallel sections with a known spacing (Howard and Reed 1998; Royet 1991).

Lengths and widths of intersections can be defined in several different ways (Fig. 1A⇓). However, the exact definition used is not that important, provided it is used consistently. In the modeling presented here, the length is defined as being the greatest dimension of the intersection. The area of each intersection is uniquely defined, but unfortunately is not the best size parameter for stereological corrections (Higgins 1994).

Some large crystals can be measured directly from the outcrop, for example, megacrysts in granites (Higgins 2000) or oikocrysts in mafic rocks. Crystals can be measured with a ruler, or their outlines traced onto an overlay for measurement later, or a mosaic of photographs can be taken for later analysis. For glacially polished surfaces the intersections with the surface are measured. But where there has been extensive weathering the megacrysts may stand out from the surface. In this case the true 3D size of the crystals can be measured. However, the minimum crystal size that can be consistently measured must be clearly established.

For rocks containing large crystals, sawn slabs can be very useful. The surface can be polished, etched, or stained to enhance the contrast between the different minerals. Two of the best known stains are sodium cobaltinitrate following a hydrofluoric acid etch, which colors all potassium feldspar a deep orange, and amaranth, which gives a deep red to plagioclase. Other minerals also can be stained (Hutchison 1974). Crystals in slabs can be measured manually as for outcrops. Slabs can also be placed on a scanner and the images analyzed automatically or manually.

For automatic image analysis, the image first must be processed to reduce it to a binary (black and white) image of the mineral of interest. This is done by selecting the color values associated with the mineral using readily available general purpose, image-processing software (Adobe Photoshop, Corel Photo-Paint), shareware (LViewPro), or dedicated programs (Image-Pro Plus, SigmaScan Pro, IPLab, ImageTool). The image can be manipulated to remove “noise” (individual pixels or small groups of pixels) or to separate individual crystals. Several images can be produced for different minerals in the rock. The binary image can be analyzed by powerful freeware programs such as NIH Image (Macintosh OS) and Scion Image (NIH Image for Macintosh OS and Windows 95/98).

Thin sections can be measured manually in several ways. The simplest is to use the scale in the microscope eyepiece to measure crystal length and width. This method has the advantage that different magnifications and orientations of the stage can be used to identify the crystal. All crystals with intersections greater than the lower cut-off limit must be measured. However, it is difficult to keep track of which crystals have been measured, although a photograph can be useful.

A more complex method involves taking photographs of parts of a thin section, or the entire section. These can then be measured in three different ways: (1) the crystal length and width can be determined using a ruler; (2) the crystal outlines can be traced onto a transparent overlay, which is then scanned and analyzed automatically as mentioned before; and (3) the photographs can be put on a digitizing tablet and the outline of each crystal in the photograph traced using the digitizing puck. Several commercial programs (e.g., SigmaScan Pro) and shareware (Measure) can be used to control the digitizing tablet and reduce the raw positional data to morphometric parameters. In each case the actual outline of the crystal in the thin section should be verified with a polarizing microscope. Thin section images also can be analyzed by automatic methods, as for slabs. The image can be acquired directly from the microscope with a digital camera, or mosaics of photographs can be scanned. Some image analysis systems allow a combination of these methods: crystal outlines can be picked out manually from the processed image.

Images of very small crystals can be obtained using reflected light or with a scanning electron microscope, using the backscattered electron image or an electron microprobe with X-ray maps. These images can be analyzed in exactly the same way as for thin sections and slabs.

No matter what method is used to acquire intersection size data, the minimum size that can be recognized and measured must be established. If no data are listed for crystals smaller than a certain size, it is important to indicate whether this reflects an artifact of measurement or a real lack of intersections in that size range. It is also important to record the details of the data-acquisition method and steps taken to ensure adequate quality control.

## Results of measurements

Most data acquisition methods produce the same result—list of crystal intersection properties, such as length, width, area, orientation, and the area of the surface measured. The first step is to calculate the frequency distribution of the area number density. CSD data commonly have an approximate logarithmic-normal distribution, hence it is best to use logarithmic size intervals for length and width measurements (Sahagian and Proussevitch 1998; Saltikov 1967). The number of size intervals (bins) per decade used will depend on the amount of data (number of crystals measured), the range in crystal sizes, and the stereological conversion method used (see later). In my experience, four to five bins per decade are usually suitable. More bins can introduce errors because there are fewer crystals in each bin and also because more cycles of correction are necessary during the stereological conversion (see later). The number of crystals in each bin is divided by the total area measured to give the area number density, *n*_{A} (*l*_{XY}), for each size interval, *l*_{X} to *l*_{Y}. [Throughout this paper, three-dimensional lengths are indicated by capital letters (e.g., *L*) and intersection lengths and widths by lower case letters (e.g., *l, w*)]

Some authors use the cumulative frequency of intersection lengths, *N*_{A} (*l*) instead of the number density. This is the sum of all intersections smaller than *l* (e.g., Cashman and Marsh 1988) or more rarely larger than *l* (Peterson 1996). The cumulative frequency distribution can be converted to *n*_{A} (*l*_{XY}) by differentiation over discrete intervals. Where data are sparse, for instance for large crystals, cumulative frequency methods have the advantage of eliminating “holes” (empty bins) in the data and smoothing it. In some cases, use of cumulative frequency methods eliminates a step in the calculations (see Peterson 1996). However, for the method proposed here, they introduce an extra, unnecessary step and hence will not be considered further.

## Stereology—conversion of 2D intersections to true (3D) dimensions

The raw data on crystal sizes are only a measurement of intersections of crystals with a plane (Figs. 1B⇑ and 2⇓). Clearly, three-dimensional size data are required for discussion of petrological processes. The problem of conversion of two-dimensional section data to three-dimensional CSDs is not simple for objects more geometrically complex than a sphere and there is no unique solution for real data. This subject is treated in a branch of mathematics called stereology. There are several reviews and books on the subject, mostly aimed at the biological sciences (e.g., Howard and Reed 1998; Royet 1991; Underwood 1970). Aspects of the problem relevant to geology have been discussed most recently by Sahagian and Proussevitch (1998), Cashman and Marsh (1988), and Peterson (1996).

Stereological solutions to this problem can be classed into direct or unbiased methods for which no assumptions about the shape and size distribution of the particles are necessary, and indirect or biased methods, in which some assumptions are needed (Fig. 3⇓). Indirect methods can be classified further into parametric methods, where a size-distribution law is assumed, and distribution-free methods that are applicable to all distributions.

In an ideal world the direct methods would be the best solution. However, although direct methods for the determination of mean size are very powerful and simple (Howard and Reed 1998), those for the determination of particle size distributions are extremely laborious because serial sections (many closely spaced sections) are necessary. Although possible in some cases, serial sectioning can be rarely done on large enough volumes of material to give statistically meaningful data. In addition, the smallest crystal that can be detected is equal to the section spacing. We are usually left, therefore, with indirect methods.

It should be mentioned that the techniques to be discussed below are based on the assumption that the section measured (or visible part of the surface) is thin compared with the dimensions of the object. For thin sections viewed in transmitted light, this assumption means that the crystals must be larger than about 0.03 mm, the thickness of a standard thin section. If crystals smaller that this limit are to be measured, then the crystal outlines are a projection and not a section, and hence different equations must be used that are not discussed here.

Many published CSD studies to date have used the following equation to convert section data to CSDs:

(1)

where *n*_{V} (*L*_{XY}) is the number of crystals per unit volume in the length interval *L*_{X} to *L*_{Y}, and *n*_{A} (*l*_{XY}) is the number of intersections per unit area in the intersection length interval *l*_{X} to *l*_{Y}. This equation was first used by Wager (1961) to convert total numbers of crystals in sections to volumetric numbers, without any discussion of its origin or justification for its use. Wager (1961) may have chosen this equation because it is the simplest way to convert the dimensions of *n*_{A} (*l*_{XY}), 1/L^{2}, to that of *n*_{V} (*L*_{XY}), 1/L^{3}. Equation 1 is not discussed in general reviews of stereological methods (e.g., Howard and Reed 1998; Royet 1991; Underwood 1970) or those applied to geological problems (e.g., Sahagian and Proussevitch 1998). Peterson (1996) considered that Equation 1 had a weaker theoretical basis than the equations that he used (see later) and produced results that lay further from the true values. In my opinion, use of Wager’s equation should be discontinued. Data should be calculated using the methods of Peterson (1996), Sahagian and Prooussevitch (1998), or those presented in this paper. Published data determined using Wager’s equation is not wasted and can be recalculated if ancillary information (shape, fabric) on the crystals and section are available.

To return to the stereological problems of conversion, two aspects are relevant here (Royet 1991; Underwood 1970): the intersection plane rarely cuts exactly through the center of each particle; hence, even in a monodisperse population of particles (one true size), intersection sizes have a broad range about the modal value, from zero to the greatest 3D length (Figs. 4A and 4B⇓). This is known as the cut-section effect. For a polydisperse distribution (many different true sizes), smaller crystals are less likely to be intersected by a plane than larger crystals. This is known as the intersection-probability effect. Both problems compound to make the stereological conversion complex.

## Indirect methods—Parametric solutions

Peterson (1996) proposed a parametric solution to rock CSDs. He assumed that rock CSDs have a strict logarithmic variation in population density for a linear variation in length. This distribution is indicated by the theoretical studies of Marsh (1988) for simple igneous systems. However, it should be noted that most sedimentary rocks have a rather different distribution that is log-normal by volume frequency (Jerram 2000). Peterson first corrected the data for the intersection-probability effect (see Eq. 3 below). He then applied a further correction using three empirical parameters that depend on the shape and spatial orientation of the crystals. He found a good correlation between the CSDs determined by his methods and the true CSDs for synthetic data with logarithmic-linear CSDs. However, many published natural CSDs are not linear (e.g., Higgins 1996; Marsh 1998; Resmini and Marsh 1995; Waters and Boudreau 1996) and published linear CSDs calculated using Equation 1 may also be curved. Peterson was aware of this problem and applied the corrections derived from linear approximations to the actual curved CSDs. However, he did not verify whether the results for curved CSDs were accurate using synthetic or natural data. Hence, it is not clear if these parametric methods can be applied to many natural systems. It also seems unwise to use a technique to find CSDs that assumes a CSD shape beforehand. Another limitation of these methods is that they can only be applied to isotropic fabrics, unless other parameters are introduced.

## Distribution-free solutions—the intersection probability effect

The intersection probability effect is easily resolved for a monodisperse (one true size) collection of spheres (Royet 1991):

(2)

where *n*_{v} is the total number of spheres per unit volume, *n*_{A} is the total number of spheres per unit area, and *D* is the diameter of the spheres. This relationship also applies to polydisperse collections of spheres, and can be modified so that it applies to polydisperse collections of other shapes if a linear measure of the size of the particle is used instead of the diameter.

The size of spheres is uniquely determined by their diameter, but the size of non-spherical objects can be defined in several ways. The simplest is the *Length* of the longest dimension of the particle and is the measure used in this paper. The *Mean Tangent Diameter* is the diameter of a sphere of equal volume. This measure is most relevant if we want to determine the total amount of material in each particle. Another measure of size is the *Mean Projected Height* (MPH; Mean Caliper Diameter), which is defined as the mean height of the shadow of the particle for all possible orientations (Sahagian and Proussevitch 1998). In isotropic materials the MPH can be approximated as the ratio of the particle volume to the mean intercept area (Sahagian and Proussevitch 1998), which is close to the longest dimension of a parallelepiped (orthogonal solid) and hence the first definition of size. The MPH can also be measured empirically using geometrical models for both isotropic and non-isotropic fabrics.

The simplest formulation of the intersection probability effect for non-spherical objects is (e.g., Royet 1991)

(3)

where *n*_{V} (*L*_{XY}) is as defined in Equation 1, *n*_{A} (*L*_{XY}) is the number of crystals per unit area in the length interval *L*_{X} to *L*_{Y}, and *H̄* _{XY} is the mean projected height for a length interval *L*_{X} to *L*_{Y}.

## Direct solutions—cut-section effect

The intersection probability effect uses the assumption that each intersection passes exactly through the center of the object along the longest axis. This is very unlikely and hence we have to contend also with the cut-section effect.

The cut-section effect can be solved analytically only for spheres (e.g., Royet 1991). The probability of a random intersection of a sphere of unit diameter, *S*_{XY}, having an intersection diameter between *d*_{X} and *d*_{Y} is

(4)

This function rises to a maximum near the diameter of the sphere, hence the most likely intersection diameter is close to the maximum (Fig. 4A⇑).

The distribution of intersection dimensions for shapes such as parallelepipeds must be derived numerically and can be very complex with several modes (Fig. 4B⇑; Higgins 1994; Sahagian and Proussevitch 1998). There are two aspects to this problem: what are the most likely intersection dimensions (length and width); and what is the form of the intersection function for smaller and larger intersections (tails)?

As noted above, for spheres and near-equant forms, the mean intersection length is close to the true 3D size of the object. However, for monodisperse populations of randomly oriented anisotropic figures, such as parallelepipeds, the most-likely intersection length is close to the intermediate dimension (Fig. 4B⇑; Higgins 1994). That is, for a particle 1 × 2 × 10 mm, the most-likely intersection length is 2 mm. The most-likely intersection width is close to the short dimension (Fig. 4B⇑). The same is true for non-rotational ellipsoids that have the same dimensions (Sahagian and Proussevitch 1998). The same argument also can be applied to polydisperse populations (many 3D sizes). These modal values for intersection length and width can be corrected to the true length of the crystal, or any other size parameter, if the shape of the solid is known. How to determine crystal shapes will be discussed later. However, there is another problem: intersection lengths and widths for a monodisperse population tail out to smaller and larger values respectively around the most-likely intersection length. If the simple conversions of Equations 2 or 3 are used, then the effects of tailing become very important. This is because small intersections of large objects (corners) will be converted using their apparent length and not their true length. This problem is illustrated for prisms and tablets in Figure 5⇓.

Saltikov (1967) proposed a method of unfolding a population of intersection lengths into the true length using a function of the intersection lengths. This method is well illustrated in the paper of Sahagian and Proussevitch (1998). For a polydisperse population, the 3D distribution of lengths can be found by applying the function for a monodisperse population of the same shape to the 2D length distribution. The values of *n*_{V} (*L*_{XY}) for a series of bins from 1, 2, 3, 4 … can be calculated sequentially from the following equations (slightly modified from those of Sahagian and Proussevitch):

(5)

where *n*_{V} (*L*_{1}) and *n*_{A} (*l*_{1}) refer to the first (largest) size interval, *P*_{12} is the probability that a crystal with a true size in interval 1 will have an intersection that falls in the interval 2, and H̄_{1} is the mean projected height for interval 1. [Note simplified notation for *n*_{V} (*L*_{1}) and *n*_{A} (*l*_{1})]. The size intervals are logarithmic—both Saltikov (1967) and Sahagian and Proussevitch (1998 ) proposed ten size bins per decade, that is, each bin is 10^{−0.1} smaller than its neighbor. The probabilities can be calculated from numerical models of different crystal shapes (Higgins 1994; Sahagian and Proussevitch 1998). For isotropic fabrics, the shape is constructed mathematically and sectioned using randomly oriented planes placed at random distances from the center of the crystal model. The mean projected heights also can be calculated using the same models.

The Saltikov method works well for spheres and near equant objects because the modal 2D length lies close to the maximum 3D length (Fig. 4A⇑; e.g., Armienti et al. 1994; Sahagian and Proussevitch 1998). However, this method is less successful for complex objects because the largest 2D intersections are not very common (Fig. 4B⇑), that is, the modal 2D length is much less than the maximum 2D length. Therefore, the less abundant, and hence less accurately determined, larger intersections are used to correct for the most abundant intersections (Fig. 6A⇓), introducing large errors in the corrected 3D length distributions. The problem is actually worse than might be expected from the works of Saltikov (1967) and Sahagian and Proussevitch (1998 ) as both authors placed the maximum intersection length at the upper limit of the size bin. Because true particle sizes can be distributed throughout the size bin, at the very least, the maximum intersection size should be placed at the center of the bin. In the CSDCorrections program, which uses the methods proposed in this paper, a more complex algorithm is used and the maximum intersection sizes are actually evenly distributed across the size interval. However, this does not make a large difference to the calculated probabilities.

Sahagian and Proussevitch (1998) reduced this problem by ignoring the class of largest intersections and shifting the probabilities to lower classes. Here, I propose using the most likely intersection length or width (modal value) to correct for the tailing to other intersections (Fig. 6B⇑). However, it is not possible to correct tailing in both directions using this technique—tailing to intersections either larger or smaller than the modal value must be ignored. It is clear from Figure 5⇑ that tailing to smaller intersections is a much more important problem than tailing to larger intersections for most CSDs. Hence, the first size interval is centered on the mode of the intersection length or width.

Some workers have considered the Saltikov technique to be impractical because of the accumulation of terms, and hence errors, in the smaller bins. This problem is not always very serious as many of the terms are not significant. It can also be reduced by using fewer, wider bins—in this paper 4 or 5 bins per decade are used. Figure 5⇑ shows how this is applied to a polydisperse collection of 1:5:5 tablets. Intersections greater than the first interval cannot be corrected precisely with this method. However, if they are added to the first interval, then this error is minimized.

Clearly, the fewer corrections that are needed the greater the accuracy of the final data. In many cases use of intersection widths instead of lengths reduces the amount of corrections necessary. Figure 7⇓ suggests a possible strategy for this choice.

Once the mode of the intersection length or width is used instead of the maximum values, then the intersection length scale no longer corresponds to the true length scale: that is, *L*_{X} ≠ *l*_{X} and *L*_{Y} ≠ *l*_{Y}. For an orthogonal block with a short dimension S, an intermediate dimension I, and a long dimension L, the modal value of the intersection lengths on randomly oriented planes is I (Higgins 1994) hence:

(6a)

Similarly for intersection widths:

(6b)

It has been found that for tablets, where L/I = 1, a much better fit is obtained to test data (see later) if the mean intersection length is used rather than L in Equations 6. This gives a maximum difference of 20% in the length scale for 1:10:10 tablets and decreases to zero for cubes. Equations 5 do not need to be changed as *H̄*_{1} is determined in terms of the true crystal length.

Other refinements have been also made to the Saltikov method. Tuffen (1998) has pointed out that the inverse of the Mean Projected Height must be averaged over the true length interval and proposed that

(7)

where *X* = upper limit of interval and *Y* = lower limit of interval. Equation 7 can be integrated to give

(8)

The original equations of Saltikov (1967) used logarithmic size bins. Sahagian and Proussevitch (1998) followed this idea, but introduced more complex equations for linear size bins. However, if the probabilities used in Equations 5 are calculated for each cycle of corrections then the calculations remain simple for any bin size. This method of calculation has been implemented in CSDCorrections, so that previously published data can be corrected. Nevertheless, the use of logarithmic size bins is strongly recommended.

The Saltikov method can only be used in its simplest form for isotropic fabrics. However, the probabilities of intersection, and the factors for converting intersection widths and lengths to true crystal sizes can be readily modeled for anisotropic fabrics if the fabric parameters are known or can be estimated. The intersection probabilities also depend on the shape of the crystals. Hence crystal shape and fabric will now be discussed.

## Crystal shapes and fabrics

The shape of crystals has its own interest (e.g., Sunagawa 1987), but is also necessary for the stereological corrections. Here overall crystal shapes are simplified and defined in terms of the aspect ratio and degree of rounding.

The aspect ratio has three parameters, the short, intermediate, and long dimensions (S:I:L). If crystals can be separated from the matrix, then their shape can be measured directly. It is not necessary to have a quantitative separation, as would be needed to determine the CSD. Typically it is not possible to separate crystals, and their shape must be estimated from the parameters of the intersection width/length ratio (*w/l*) distributions. Higgins (1994) showed that for a massive rock (no preferred orientation) and crystals modeled as parallelepipeds, the mode of the *w/l* is equal to S/I. The ratio I/L is more difficult to determine, but can be estimated from the skewness of the *w/l* distributions as follows:

(9)

and

These equations are not accurate for near equant shapes. If the rock is foliated or lineated, then the mode of the *w/l* of sections parallel or normal to the fabric can give much more precise values for both S/I and I/L (Table 1⇓).

The degree of rounding can vary from euhedral (angular) to ellipsoid. Concave and branching crystal forms cannot be corrected by the methods proposed here and must be measured by direct stereological methods. A special feature of the intersection distribution graphs for parallelepipeds is that tailing to intersections smaller than the mode are almost independent of length or width. The same distribution for spheres ramps down to zero abundance at zero length or width (Fig. 4A⇑). Hence the effect of rounding can be approximated by modifying the part of the distribution to the left of the mode. The exact form of the function is not clear, but has been approximated here in the following way: a roundness factor is defined to be between one (ellipsoids) and zero (angular crystals). The function follows that of the parallelepiped down to the point where the roundness factor equals the ratio of I/S. From then on the function decreases linearly to zero at zero size intersections.

Fabric in rocks can be defined by a preferred orientation of crystal shapes or lattice orientations. Only shape-preferred orientations will affect the stereological conversion of CSDs, although lattice-preferred orientations have interest in their own right. Fabric can be lineated or foliated, and the quality can vary from perfect alignment of all crystals to no alignment (equivalent to a massive rock). The effects of preferred orientations on the probabilities of intersection dimensions can be calculated from a numerical model as before, except that the orientation of the intersection planes is constrained according to the type and quality of the fabric and the orientation of the section.

## Population density

All methods so far described have produced a set of values of the volumetric number density *n*_{V} (*L*_{XY}). However, the value of this parameter depends on the width of the interval *L*_{Y} - *L*_{X}. Hence it is usually divided by the width of the interval to give the population density .

(10)

Note that the units of population density are 1/Length^{4}, as the volumetric number density of crystals has the units of 1/Length^{3} and it is divided by a length, the bin width.

If the CSD is continuous, that is, it has no empty bins, then the value of population density will not change for different bin widths. If the gaps appear to be due to scarcity of data then the bin widths may be widened or more crystals counted. However, there is no intrinsic reason why all CSDs should be continuous and so we must be able to accommodate such gaps. (However, it should be noted that intersection histograms are usually continuous.) Although commonly drawn as line graphs or as a series of unconnected points, CSDs are actually histograms (Fig. 8⇓). Therefore, the line type CSDs should not be drawn across empty bins. Similarly, unconnected points give no idea if there are intervening empty bins.

The most convenient units of population density for most studies are 1/mm^{4}, although some authors (e.g., Cashman and Marsh 1988) use mm for linear dimensions (Size) and 1/cm^{4} for population density. CSDs are generally plotted on a graph of Ln (population density) versus Size following Marsh (1988).

## Error analysis

Inaccuracy in the determination of population densities arises principally from three sources. The easiest to understand and quantify is the counting error. This is taken to be the square root of the number of intersections within an interval. It is usually only significant for larger size intervals with fewer than 20 intersections.

The second source of error is in the value of the probability parameters *P*_{XY} used in Equations 5. Although these parameters are precisely defined for fixed convex shapes, crystals in most natural systems have more irregular and variable shapes. Another source of error is that tailing to intersections larger than the modal interval is included in the modal interval (Fig. 6B⇑). Hence, it is difficult to estimate the contributions from this source to the total error. However, it is easy to calculate the contribution of the counting errors of other intersection intervals to the total correction of an interval. This source of error is most important for small size intervals, where corrections are most significant.

The third source of error lies in the conversion of intermediate crystal dimensions (for intersection length measurements) or short crystal dimensions (for intersection width measurements) to crystal lengths (the size criteria used in this paper) using Equations 6. Errors in the determination of the crystal shape will produce systematic errors in both the population density and the size.

## Verification of the method

The method of data reduction proposed here must be verified to show that it does yield results that are not only precise, but also accurate. The key algorithm of the program CSDCorrections has been inverted to produce the crystal intersections shown in Figure 2⇑. Such intersections could then be measured to check that the original population density is recovered. However, an independent check on the methods is much more valuable. The modeling of Peterson (1996) was based on a very different approach, and should provide such a check.

Peterson presented data for populations of spheres and orthogonal parallelepipeds with linear CSDs and several different slopes. Lists of intersection lengths were obtained at http://www.nrcan.gc.ca/~tpeterso/csd_e.html. Figure 9⇓ shows the results for two families of solids with very different slopes and intercepts. In both cases there is a close agreement between the original population density and that calculated with the program CSDCorrections, with maximum errors of ± 0.5 ln units in the population density.

## Application of the method to rocks

Published CSDs calculated according to other conversion equations can be readily recalculated if the crystal shape, degree of roundness, fabric, and orientation can be estimated.

The CSDs of two lavas and two plutonic rocks are shown in Figure 10⇓ and the original data in Table 2⇓. Cashman and Marsh (1988) determined the CSD of plagioclase from the Makaopuhi lava lake, Hawaii. These data have been converted to mm and recalculated using crystal dimensions of 1:3:4 and massive fabric (Fig. 10A⇓). For larger crystals there is a difference of a factor of 50 in the population density [Ln(4)]. There is much less difference in the overall slope, 0.75, but the new CSD has a more-pronounced curvature. The CSD of microcline megacrysts from the Cathedral Peak Granite, Sierra Nevada (Fig. 10B⇓; Higgins 2000) have population densities one fiftieth [Ln(−4)] that calculated using Wager’s equation, and the turn-down for small crystals is more evident. The CSDs of plagioclase crystals from the Kameni islands, Greece (Fig. 10C⇓; Higgins 1996) have somewhat different slopes using the two equations, but actual values converge for greater lengths. Plagioclase from the Kiglapait Intrusion, Labrador has a pronounced hump-shape typical of textures produced by textural coarsening (Higgins, in preparation). In this case both methods give exactly the same result (Fig. 10D⇓).

## Concluding remarks

The stereology of conversion of the two-dimensional data obtained from the study of thin sections etc., is complex, but can be approximated for convex shapes. The approach here is based on a parallelepiped model with an aspect ratio the same as that of the crystals. It is modified to account for rounding of the crystals. The fabric of the rock must also be determined as well as the quality of the fabric and the orientation of the section with respect to the fabric. CSDs calculated by other popular methods can have population densities that are very different from those produced by the methods suggested here. However, the slope of many CSDs does not change much. Published CSDs calculated using other methods can be recalculated if the crystal shape and fabric can be estimated. Finally, possible advances in stereology make it essential that in all studies, the raw 2D data of CSD measurements together with the areas measured, crystal shapes, and fabrics are made available, so that the data may be recalculated in the future.

The calculations described above are incorporated into a Windows 95/98 program called CSDCorrections. Data are entered either as lists of intersection widths or lengths, or as intersection interval lengths and numbers of intersection lengths, or widths per unit area. The latter method is to accommodate recalculation of published data. In addition, the crystal shape (block or sphere), overall dimensions, degree of roundness, fabric, quality of fabric, orientation of section, and area measured must also be entered. The program CSDCorrections, other programs, and information on CSDs, are available from the American Mineralogist web site http://www.minsocam.org/MSA/AmMin or http://wwwdsa.uqac.uquebec.ca/~mhiggins/CSD.html.

## Acknowledgments

I thank the Université du Québec à Chicoutimi for awarding me a sabbatical and Tim Druitt and the Département de Géologie, Université Blaise Pascal, France for welcoming me. Discussions with Hugh Tuffen helped the development of the program. I also thank Sarah Barnes, Bruce Marsh, Dougal Jerram, and two anonymous reviewers for their comments. This research was supported by the Natural Science and Engineering Research Council of Canada.

## Footnotes

Paper handled by Robert F. Dymek

- Manuscript Received October 7, 1999.
- Manuscript Accepted March 25, 2000.