Mineralogica l Magazine, October 2000, Vol. 64(5), pp. 945–959

Three-dimensional study on the interconnection and shapeof crystals in a graphic granite by X-ray CT and image analysis
S. IKEDA1 , *,{, T. NAKANO2 AND Y. NAKASHIMA3
1

Geological Institute, Graduate School of Science, University of Tokyo, Hongo, Bunkyo-ku , Tokyo 113-0033 ,Japan
2Geological Information Centre, Geological Survey of Japan, Higashi, 1-1-3, Tsukuba, Ibaraki 305-8567 , Japan3
Geophysics Department, Geological Survey of Japan, Higashi, 1-1-3, Tsukuba, Ibaraki 305-8567 , Japan

ABSTR ACTThe technique of investigating 3-dimensional interconnections and the shapes of crystals in a rock by
X-ray computerized tomography (CT) and image analysis was developed using a graphic granite
specimen as an example. Fifty 2-dimensional tomographic images (slices) of the graphic granite wereobtained ‘non-destructively’ using a medical X-ray CT scanner. Since a CT value of the specimen was
decreased with increasing cross-sectional sample area by the effect of beam-hardening, the CT value
was corrected using the area of each slice. Binary images of the slices were made comparing one ofthem with a thin-section of the slice. Using the binary images, connection analysis of quartz rods in thegraphic granite specimen was performed on the basis of percolation theory (cluster labelling). This
analysis showed that at least 89.9% of the quartz rods were connected in three dimensions.
Furthermore, the 3-dimensional shape of the quartz rods was analysed using the 2-point correlation
function calculated from the binary images. The average shape of the quartz rods was obtained byfitting an ellipsoid to the high-value region of the 2-point correlation function. The elongation axis of
the ellipsoid agreed well with the crystallographic c-axes of the quartz rods.
K EY WORDS : X-ray CT, 3-dimensional image analysis, graphic granite, cluster labelling, 2-point correlation

function.

Introduction
IN the Earth Sciences, textures in rocks are
examined in two dimensions using thin-sections
or polished sections. Sometimes, 2-dimensional
investigation is enough to quantify 3-dimensional
characteristics of texture. For example, if anappropriate model of crystal shape is assumed, we
can determine a 3-dimensional crystal size
distribution from a 2-dimensional distribution bystereology (e.g. Cashman and Marsh, 1988 ).
However, many features of the 3-dimensional
structure cannot be investigated by 2-dimensional* E-mail: ikeda@chem. s.u-tokyo.a c.jp{
Present address: Department of Complexity Scienceand Engineering, Graduate School of Frontier Sciences,
University of Tokyo# 2000 The Mineralogical Society

observation. A 3-dimensional network structure,such as a pore network, is often interpreted as
separate pores by 2-dimensional observation.
Also, interconnection of crystals in three dimensions cannot be recognized by 2-dimensional
analysis. The interconnections of crystals
contain signi cant information relating to theorigins of the rocks, e.g. the process of magma
differentiation (Philpotts et al., 1998). Therefore,
it is important to develop the technique to analyse
a 3-dimensional network structure.
There are some techniques to investigate the 3dimensional internal structure of rocks. For
example, Bryon et al. (1995), Cooper andHunter (1995) and Marschallinger (1998) demonstrated 3-dimensional reconstructions of theinternal structures in rock samples using successive digitized images acquired by serial thinsectioning or lapping. Philpotts et al. (1998)

S. IKEDA ETAL.

discovered the 3-dimensional network structure(chains) of plagioclase phenocrysts in a basalt
specimen by a partial melting experiment.
However, the technique of serial polishing needs
skill and time. Also, the techniques of partial
melting or selective dissolution with acids or
alkalis are not appropriate for every sample.
Nowadays, we can observe the internal
structures of rocks by non-destructive methods,
such as X-ray computerized tomography (CT).
X-ray CT was invented and developed in the
medical  eld in the 1970s and has been applied togeological samples in recent years. For example,
by the X-ray CT method, we can analyse directly
the 3-dimensional crystal-size distribution, crystal
spatial disposition (Denison et al., 1997; Denison
and Carlson, 1997), and structure of cracks or
pores (Raynaud et al., 1989; Nakano et al., 1992;
Nakashima et al., 1997). However, the techniques
to investigate the 3-dimensional connection andshape of crystals have not become well established. Furthermore, X-ray CT images generally
have an error resulting from beam-hardening (seebelow). This error affects precise analysis. In thisstudy we developed a technique to correct theerror and to investigate the 3-dimensional
connection and shape of crystals in a rock by

combining the X-ray CT method and image
analysis. The connection analysis was based on
the technique of pore-network analysis proposed
by Nakano et al. (1992). Graphic granite wasselected as an example for these analyses because
quartz rods in graphic granite are known toconnect to each other as illustrated in thefollowing section.
Specimen and interconnection of quartz rods
We used the specimen of graphic granite shownin Fig. 1 (Sample GSJ-R6052, from theGeological Survey of Japan). This specimen
was obtained by H. Matsubara (Matsubara, 1956)
from the Shiozawa district in Ishikawa town,
Fukushima Prefecture, Japan. Around Ishikawa,
there are many granitic pegmatites within agranitic body intruded into the Abukum ametamorphic belt. The specimen was collectedfrom a graphic granite layer in one of thepegmatite bodies (Matsubara, 1956). The K-Ar
age of the surrounding granitic body is
90 100 m.y. (Kawano and Ueda, 1967).
The specimen consists of quartz rods and hostalkali-feldspar (microcline microperthite). In thespecimen, most of the quartz rods show

FIG. 1. Specimen of graphic granite (GSJ-R6052) with coordinat e system for the X-ray CT (left), and three slices(right). Fifty slices parallel to the XY-plane were obtained along the Z-axis. The dark parts are quartz rods both in thephotograp h of the specimen and in the CT images. The areas surrounde d by black lines in the CT images were used
for the three-dimensiona l analysis and correspond to the areas of images shown in Fig. 4.

946

SHAPE OF CRYSTALS IN GRAPHIC GRANITE

FIG. 2. Photographs of the thin-section of the X-ray CT image no. 16. in crossed polars. (a) Simultaneous extinctio nof alkali feldspar host. The bright parts are quartz grains. (b) Simultaneous extinctio n of most of the quartz rods. The
quartz rod surrounde d by the broken line has three domains with different crystallographi c orientation s. The dark
domain in the rod shows the simultaneous extinctio n with most of the quartz rods. The dark circular region near thecorner of the box (broken line) was lost in the process of making the thin-section.

simultaneous extinction (Fig. 2) as reported inmany previous studies (e.g. Wahlstrom, 1939;
Simpson, 1962). Furthermore, we con rmed,
using a precession camera, that the elongation
direction of the quartz rods correspond approximately with the crystallographic c-axis of quartz.
In order to  nd the cause of the simultaneous
extinction, Simpson (1962) performed successive
surface grinding of a graphic granite specimen
and showed interconnections of the quartz rods.
He concluded that the simultaneous extinction
occurs because all quartz rods are parts of thesame crystal. Therefore, we concluded thatgraphic granite was a good choice of material toinvestigate 3-dimensional crystal interconnection.
X-ray CT
X-ray CT is a non-destructive technique used todetermine tomographic maps of the X-ray linear
attenuation coef cient in an object. We describe
the 2-dimensional X-ray CT images as ‘slices’.
The brightness in the slices (e.g. in Fig. 1) shows
the linear attenuation coef cient, which depends
on the density and chemical composition of theobject and the energy of the incident X-ray beam
(Koch and MacGillavry, 1962; McCullough,1975; Denison et al., 1997).
The linear attenuation coef cient obtained by
X-ray CT is generally called the ‘CT value’. It is adimensionless number normalized by the linear
947

attenuation coef cient of water as de ned by thefollowing equationm mw
CT value1000
1mw
where m and mw are the linear attenuationcoef cients of the object and water, respectively.
We used a medical X-ray CT scanner (Hitachi
CT-W2000) to obtain the slices of the specimen.
The scanner belongs to the third-generation typeadopting a continuous rotation method. The
continuous (polychromatic) X-ray beam from atarget of Mo-W alloy was used. The voltage andcurrent of the X-ray tube were 120 kV and
175 mA respectively. The X-ray exposure timefor each slice was 4.0 s. We used the 3-dimensional coordinate system for acquisition of the
X-ray CT images shown in Fig. 1 (XYZ-axes).
The slices lay in the XY-plane. The size of thevoxel (volume element in three dimensions;corresponding to ‘pixel’ in two dimensions) on
the XY-plane was 0.31360.313 mm2 , which wasconsidered as the in-plane resolution in the slices.
The thickness of each slice was 1.0 mm. The
specimen was automatically moved along the
Z-axis by 0.5 mm after acquisition of each slice.
Thus each slice overlaps by 0.5 mm in thickness
with the adjacent slices. Fifty 16-bit images(slices) along the Z-axis were obtained in all. It
took 28 min to acquire 50 slices.

S. IKEDA ETAL.

With the graphic granite specimen, the CT
value of alkali feldspar (~2080) was greater thanthat of quartz (~1980), though the density of alkali
feldspar is lower than that of quartz. The higher
CT value of alkali-feldspar depends on highattenuation by K (Raynaud et al., 1989). As
described above, the linear attenuation coef cientis a function not only of density but also ofchemical composition and the energy distributionof the X-ray beam. We con rmed the attenuationcoef cients of alkali feldspar and quartz using theformulae and attenuation coef cients of
McCullough (1975) and Hubbell et al. (1980) at
an X-ray energy of <120 keV.

is homogeneous. This beam-hardening causes
two types of errors in the CT value. One error
arises from the heterogeneity of the CT value
within each slice. The CT value at the centre of
the slice is lower than that at the periphery evenin the homogeneous material. The heterogeneity
within each slice was not seriously affected byomitting the periphery from the analysed area
(Fig. 1 right). The other error is the difference in
the CT value between different slices. A slicewith large cross-section area (e.g. no. 50 inFig. 1) has a lower CT value than a slice withsmall area (e.g. no.1 in Fig. 1). In this study, thedifference in the CT values amongst the slices
was very large, and we had to correct for it. As
shown in detail later, we found a simple
relationship between the error (difference ofCT value) and the cross-sectional area of theslice. The increase of the CT value wasproportional to the decrease of the area of theslice. We corrected the CT value of each slice inproportion to the decrease of the area.

Image analysis techniques
Correction of beam-hardening
The CT value obtained using a continuous(polychromatic) X-ray beam is in uenced by
the beam-hardening (Denison et al., 1997). The
low-energy component of X-ray photons is moreabsorbent than the high-energy component.
Thus, most of the low-energy photons areabsorbed at the periphery of a specimen justafter entry of the X-ray beam into the specimen.
In the X-ray CT technique, image reconstruction
is performed using the sum total of the photons
belonging to all energy levels as intensity data;thus the linear attenuation coef cient at the
periphery, where many photons are absorbed, iscalculated as a high value even if the specimen

Thresholding
In order to perform 3-dimensional analysis, wemade binary images of the  fty slices using an
ordinary thresholding (segmentation) technique ofa grey level image (e.g. Dilks and Graham, 1985),
and identi ed the quartz rods as follows. Firstly,
the binary image of slice no. 16 was made bycomparing it carefully with the thin-section
corresponding to the slice (Fig. 2). Since the CT
value at the edge of the specimen becomes muchlarger owing to the effect of the beam-hardening(and because of the Gibbs’ phenomenon relating
to the algorithm of image reconstruction) as
described above, we cut out an area of
28628 mm2 (90690 voxels) from the originalslice and used it for the 3-dimensional analysis.
Binary images of the other 49 slices were made in
the same way. The cut-out areas on several slices
are shown in Fig. 1 (right).
Connection analysis (Cluster labelling)

FIG. 3. De nition of connectio n of voxels occupied byquartz crystal in the three-dimensiona l image.

948

We developed software to visualize and analyse
objects from 3-dimensional X-ray CT images(Nakano and Fujii, 1989a, 1991; Nakano et al.,
1992, 1997). The binary images were used toanalyse the connection of quartz rods. In
percolation theory, the analysis is known as the‘cluster labelling method’ (Stauffer, 1985). In thisanalysis, connection of two adjacent voxels in theSHAPE OF CRYSTALS IN GRAPHIC GRANITE

FIG. 4. Grey level images and binary images before and after the correctio n of CT value. Only 11 of 50 slices areshown (every 5 slices). Each image consists of 90690 = 8100 voxels (28628 mm2 ). (a) 16-bit grey level images
before the correction. (b) Binary image before the correction. (c) 16-bit grey level images after the correctio n. (d)Binary image after the correctio n. In the binary images, the black region indicates the quartz rods which belong to
the largest cluster while the grey indicates the quartz rods which do not belong to the largest cluster in Fig. 6f.

quartz rods is evaluated as shown in Fig. 3. Whentwo voxels shared a face, we considered them as
connecting voxels. When two voxels contacted
only at a vertex or an edge, we considered them as
separate voxels. This criterion is commonly used in
the connectivity analysis (Stauffer, 1985; Nakano
and Fujii, 1989b). We call connected voxels a‘cluster’. The extent of interconnection of quartz
rods was assessed quantitatively using the term‘connectivity’ which is a volume fraction relating
the largest cluster of quartz rods to the volume ofall quartz rods (Nakano and Fujii, 1989b):
Connectivity

Volume of the largest cluster
Total volume of all quartz rods

spatial correlation features. One is the longdistance correlation such as the spatial distribution of minerals in rocks (Morishita and Obata,
1995; Morishita, 1998). In this case, we canobtain information on the degree of randomnessfor individual minerals, the mean distance
between different minerals, etc. The otherfeature of the function is the short-distancecorrelation about the origin, and it represents theaverage shape of an object. We used the latterfeature of the function. In the graphic granite
specimen shown in Fig. 1, the quartz rods seem toextend (grow) in the same direction; thus it was
expected that the 2-point correlation function of
the quartz rods would show meaningful resultswith respect to their average shape. Furthermore,
it is known that the 2-point correlation function is
not strongly in uenced by the noise of the image
data; thus it was expected that we could perform
precise shape analysis.
In binary images, each voxel has a value of 1 or
0. Here we de ne the value of each voxel B(~r) inthe ~r space as follows:2

From this de nition, the connectivity becomes
1 when all quartz rods are connecting in thespecimen.
Shape analysis
In order to analyse the 3-dimensional shape ofquartz rods, a 2-point correlation function (e.g.
Berryman, 1985; Berryman and Blair, 1986) was
calculated in three dimensions. The 2-point
correlation function represents two signi cantB(~r) =

949

when the voxel was in the quartz rods
{ 10 whenthe voxel was not in the quartz rods

S. IKEDA ETAL.

FIG. 5. (a) Histograms of the CT value for 50 slices (grey level) of 16-bit CT images shown in Fig. 4a (before the CT
value correction ). The horizonta l axis indicates the CT value. The vertical axis indicates the frequency of the voxels.
The solid circle on each histogram indicates the peak level of alkali feldspar. (b) Variation of cross-sectiona l area ofthe specimen with slice no. The areas correspond to the number of voxels with CT value >1000.

The 2-point correlation function of quartz is
expressed as
AC d

B r

B r

d

r

4

AC(0,0,0) using the method shown in the
Appendix.
Results and discussion

where the brackets r indicate a volume averageover the spatial coordinate r. In three dimensions,
AC(d) is written as AC(dx, dy, dz). AC(0,0,0),which is the value at the origin, shows the highest
value and represents the volume fraction of theobject (i.e. the quartz rods). Around the origin,
there is a spatial  eld having high AC(dx, dy, dz)
values, and the  eld is considered to represent the
average shape of the objects. In this study, anellipsoid was  tted to the  eld whose AC(dx, dy,dz) values were greater than or equal to half ofData correction and thresholding
Three slices of a 2-dimensional X-ray CT image
are shown in Fig. 1 (right). As described above,
there was a large difference in the CT value
between the slices. Figure 4a shows the 11 slices(every 5 slices) of grey-scale images (90690
voxels) before the CT value correction. If wecarry out thresholding with a  xed threshold value(grey level), we obtain the binary images shownin Fig. 4b. The slices having a larger slice number
show a higher volume fraction of quartz rods.

950

SHAPE OF CRYSTALS IN GRAPHIC GRANITE

This trend is not consistent with what we see inthin-sections. From our observations of several
thin-sections, the volume fraction of quartz rodsin the specimen is almost the same from slice no.
1 to no. 50. Figure 5a shows histograms of the CT
value for all slices. The solid circle at the top of
the each histogram indicates the CT value ofalkali feldspar. We recognize the systematic shift
in peak position of alkali feldspar. Figure 5bshows the variation of the cross-sectional area ofthe slices. The areas of the slices increase withincreasing slice number. Since the variation of thearea is almost linear with the slice no. (Fig. 5b),
we considered that the linear shift in the peak
position of the histogram (Fig. 5a) was caused by
the beam-hardening, depending on the variationof the cross-sectional sample area. Thus wecorrected the CT value of all slices by shifting
the histograms as shown in Fig. 5a, on theassumption that the peak position of alkalifeldspar must be the same in all slices. The
corrected (shift) value of each slice wasdetermined as (maximum corrected value)6(50– slice no.)/49.
The corrected grey level images are shown in
Fig. 4c and the binary images obtained by thethresholding are shown in Fig. 4d. The thresholdCT value extracting only quartz (2024) was
determined by comparing slice no. 16 with its
thin-section. After the data correction, the area ofthe quartz rods become almost the same for all
slices (Fig. 4d).
Connection analysis (cluster labelling)
Figure 6 shows the result of the 3-dimensional
connection analysis (cluster labelling) using thebinary images. The sub- gures show a bird’s-eyeview reconstructed using: (a) 1 slice; (b) 5 slices;(c) 10 slices; (d) 20 slices; (e) 30 slices; and(f) 50 slices, respectively. The left column of this gure shows all the quartz rods and the rightcolumn shows only the largest cluster ofconnected quartz rods. The sub- gures in theright column (a f) indicate the results ofconnection analysis calculated individuallyusing the data from 1, 5, 10, 20, 30 and 50
slices, respectively. Results (a e) are not part ofresult f analysed using 50 slices. We havede ned the ‘connectivity ’ in equation 2 as thevolume fraction of the largest cluster of thequartz rods (Fig. 6, right) to all quartz rods
(Fig. 6, left). The calculated connectivities were(a) 0.133, (b) 0. 160, (c) 0.337, (d) 0.379,

FIG. 6. ‘Bird’s eye view’ representin g 3-dimensiona lconnection analysis (cluster labelling ) of quartz rods.
The left column shows all the quartz rods in the analysedarea while only the largest cluster of quartz rods isshown in the right column. Each  gure shows the image
constructe d using slices of (a) no. 1, (b) no. 1 5, (c) no.
1 10, (d) no. 1 20, (e) no. 1 30 and (f) no. 1 50,
respectively. Numerical values indicate the connectivit yof quartz rods de ned by equation 2. In (a f), the top
slice is the CT image no. 1. The XYZ-axes are identicalto those of Fig. 1. The spatial extent of (f) is 28 mm(X)628 mm (Y)625 mm (Z).

951(e) 0.794 and (f) 0.899, respectively. With up
to  ve sl ic es (Fi g. 6b; 2. 5 mm in t he
Z-direction), connectivity did not increase.

S. IKEDA ETAL.

Between (b) and (c), connectivity started toincrease. With 50 slices (Fig. 6f; 25 mm in the Zdirection), almost all quartz rods (89.9%) were
found to be connected. In Fig. 6f, the separate
quartz rods (10.1%), which did not belong to thelargest cluster, were located at the corners of theanalysed space (28628625 mm3 ). These separate rods are also shown in Fig. 4d in grey, andwe can con rm that they are located at thecorners. Thus, such quartz rods would be
included in the largest cluster if the analysed
space became larger. This result is consistent
with the hypothesis of Simpson (1962) that allquartz rods are parts of the same crystal. By
image analysis we were able to assess theinterconnection quantitatively.
We could also investigate the structures of thejunctions of connected quartz rods in XZ-planes
synthesized by image processing (Fig. 7) though
we could not clearly con rm those structures in
XY-planes (X-ray CT images). Figure 7 shows
how the quartz rods join or separate. Suchinformation will be important in revealing the
crystal growth mechanism of the quartz rods andthe origin of graphic granite. These subjects (e.g.
Fenn, 1986; Lentz and Fowler, 1992; Stel, 1992)
will be discussed in a separate paper. By image
analysis, we can obtain tomographic (2-dimensional) images of sections cut in any orientationand it becomes easy to investigate the structures.
This is one of the signi cant points that advanced
the technique compared to the serial surface
grinding of Simpson (1962).
Thus, the combination of X-ray CT and
3-dimensional image analysis is a powerful toolin analysing the 3-dimensional interconnection of amineral. However, we cannot identify the crystallographic orientation by these techniques. For
example, even the quartz rods that did not show
simultaneous extinction in Fig. 2b were also
thought to be components of the largest clustershown in Fig. 6f. However, this result is not a
mistake of the above connection analysis. The
quartz rod surrounded by the broken line in Fig. 2b
consists of two domains some of which show

simultaneous extinction with some of the quartz
rods and some of which do not. Therefore, thequartz rods which do not show simultaneous
extinction also touch the largest cluster, and theconclusion on the connection was true. The
connectivity obtained by X-ray CT carries noinformation on the crystallographic orientations
because the X-ray CT technique is not based on
XRD but on X-ray absorption. Conversely, this
aspect will be a good point to analyse theinterconnection of poly-crystals having different
crystallographic orientations, e.g. chains (network
structure) of plagioclase phenocrysts in basalt
discovered by Philpotts et al. (1998) in a partial
melting experiment.
Shape analysis
We calculated a 2-point correlation function
AC(dx, dy, dz) using the binary images in orderto analyse the 3-dimensional shape of the quartz
rods (Fig. 8). The value at the origin AC(0,0,0) was
0.232, which is the highest value and represents thevolume fraction of quartz. In the dz = 0 slice, we
can recognize the high AC(dx, dy, dz) area aroundthe origin (bright area). With decreasing dz, thearea having the high AC value moves away fromthe centre point, and the AC value decreases
gradually. Since the 2-point correlation function issymmetrical with respect to the origin (dx, dy, dz)= (0,0,0), the bright area moves in the opposite
direction in the slices for +1 4 dz 4 +10, thoughthese slices are not shown in Fig. 8. We carried out
the ellipsoid  tting to the area with high correlation
value 50.116 (bold contours in Fig. 8), which was
half of AC(0,0,0). The  tting was performed byminimizing or maximizing the 2nd moments
around the axes (Appendix).
The result of the ellipsoid  tting is shown in
Fig. 9. The shape of the quartz rods could beapproximated to an elongated ellipsoid with large
aspect ratio. The elongation axis of the ellipsoid isconsidered to correspond with the elongation axesof the quartz rods with preferred orientation. Thiscan be con rmed by comparing Figs 9 and 6f.

FIG. 7. Junction s of quartz rods recognize d in binary images on the XZ-plane. The images were synthesize d by acomputer using 50 slices (XY-plane) obtained by X-ray CT. The numbers above the images indicate slice numbers
along the Y-axis.

952

SHAPE OF CRYSTALS IN GRAPHIC GRANITE

FIG. 8. Three-dimensiona l distribution of two-point correlatio n function AC(dx, dy, dz) of the quartz rods. Eachdiagram shows a slice of AC(dx, dy, dz) for  xed dz. The ranges of dx, dy and dz are 18 4 dx 4 +18, 18 4 dy 4
+18 and 10 4 dz 4 0, respectively , and the centre point of each diagram is (dx, dy) = (0,0). Since the function issymmetrical with respect to the origin (dx, dy, dz) = (0,0,0), the slices for +1 4 dz 4 +10 can be derived from theslices shown in the  gure. Contour lines are drawn at an interval of 0.0232 which is 1/10 of AC(0,0,0). The bold line
indicates AC(dx, dy, dz) = 0.116 which is half of AC(0,0,0).

953

S. IKEDA ETAL.

values based on the cross-sectional areas of thespecimen.
(2) The cluster labelling technique showed atleast 89.9% of the quartz rods to be connectedwith each other in the specimen. This result ofinterconnection of quartz rods was consistent withthe previous study by Simpson (1962).
(3) The 2-point correlation function was
calculated in order to investigate the 3-dimensional shape of the quartz rods. The ellipsoid tted to the high-value area of the function
represented a very elongated shape whose
elongation axis almost corresponded to the
crystallographic c-axes of the quartz rods.
Furthermore, we con rmed the following
useful points: (1) We can obtain quantitative
information on the connection of crystals using
the ‘connectivity’ value. (2) We can observe theprecise structure of junctions of connected
crystals from various angles by image processing.
(3) We can obtain the average shape andelongation axis of crystals by 2-point correlation
function and an ellipsoid  tting. (4) All of theseprocesses can be performed 3-dimensionally andnon-destructively. (5) The techniques are nottime-consuming.

FIG. 9. The average three-dimensiona l shape of quartz
rods in the specimen derived from an ellipsoid  tting to
the high value area of two-point correlation function(AC(dx, dy, dz) 5 0.116) shown in Fig. 8. The
rectangula r cell correspond s to that of Fig. 6f.

Acknowledgements
According to the crystallographic analysis by aprecession camera, this elongation direction (i.e.
c-axis of the  tted ellipsoid) corresponds approximately to the crystallographic c-axis of the quartz
rods.
Such image analysis of crystallographic
preferred orientation using the elongation axis
has been carried out in only two dimensions (e.g.
Simigian and Starkey, 1986). Thus this shapeanalysis by 2-point correlation function andellipsoid  tting is useful not only for investigationof shape but also for the crystallographic
orientation analysis in three dimensions.
Conclusions
We developed techniques to analyse 3-dimensional connection and the shapes of crystals by
X-ray CT and image analysis. Graphic granite
was the material selected for these analyses. The
conclusions in this study are summarized asfollows:(1) The difference in CT values from one slice
to another was recognized. This difference is dueto the beam-hardening, and we corrected the CT954

We would like to thank C. Matsue for sample
preparation and M. Shimizu and Y. Yabe, whoprovided useful information on the geology ofpegmatites in Ishikawa town. We are deeply
indebted to T. Ozawa and M. Nespolo for their
help in the determination of crystallographic
orientations using a precession camera. The
analysis with the precession camera was performedat the Mineralogical Institute, University of Tokyo.
We are grateful to A. Tsuchiyama and T.
Hanamoto for their kind suggestion with respect
to the theory of X-ray absorption. Helpful advice
from K. Nakamura, T. Ohtani, M. Toriumi, H.
Nagahara, and T. Hiraga is gratefully acknowledged. Thanks are due to S. Sasaki for allowing us
to use a personal computer (CD-R writing system).
A review by an anonymous referee helped to
improve the manuscript. This research was
supported in part by a Research Fellowship of
the Japan Society for the Promotion of Science for
Young Scientists, and a Grant-in-Aid for Scienti cResearch from the Ministry of Education, Science
and Culture, Japan (to S.I.) and by the grant fromthe Agency of Industrial Science and Technology(to Y.N.).

SHAPE OF CRYSTALS IN GRAPHIC GRANITE

References
Berryman, J.G. (1985) Measurement of spatial correlation functions using image processing techniques. J.
Appl. Phys., 57, 2374 84.
Berryman, J.G. and Blair, S.C. (1986 ) Use of digital
image analysis to estimate  uid permeabilit y ofporous materials: Application of 2-point correlatio nfunctions. J. Appl. Phys., 60, 1930 8.
Bryon, D.N., Atherton, M.P. and Hunter, R.H. (1995)
The interpretatio n of granitic textures from serial
thin sectioning , image analysis and three-dimen sional reconstruction . Mineral. Mag., 59, 203 11.
Cashman, K.V. and Marsh, B.D. (1988) Crystal size
distributio n (CSD) in rocks and the kinetics anddynamics of crystallizatio n II: Makaopuhi lava lake.
Contrib. Mineral. Petrol., 99, 292 305.
Cooper, M.R. and Hunter, R.H. (1995) Precision serial
lapping, imaging and three-dimensiona l reconstruc tion of minus-cemen t and post-cementatio n intergranular pore-system s in the Penrith Sandstone ofnorth-wester n England. Mineral. Mag., 59, 213 20.
Denison, C. and Carlson, W.D. (1997) Three-dimensional quantitativ e textural analysis of metamorphic
rocks using high resolutio n computed X-ray tomography: Part II. Application to natural samples. J.
Metam. Geol., 15, 45 57.
Denison, C., Carlson, W.D. and Ketcham, A. (1997)
Three-dimensiona l quantitativ e textural analysis ofmetamorphic rocks using high-resolutio n computed
X-ray tomography: Part I. Methods and technique s.
J. Metam. Geol., 15, 29 44.
Dilks, A. and Graham, S.C. (1985 ) Quantitativ emineralogica l characterizatio n by back-scattere delectron image analysis. J. Sed. Petrol., 55, 347 55.
Fenn, P.M. (1986) On the origin of graphic granite.
Amer. Mineral., 71, 325 30.
Hubbell, J.H., Gimm, H.A. and Øverbø, I. (1980) Pair,
triplet and total atomic cross sections (and massattenuatio n coef cients) for 1 MeV 100 GeV
photons in elements z=1 to 100. J. Phys. Chem.
Ref. Data, 9, 1023 145.
Kawano, Y. and Ueda, Y. (1967) K-Ar dating on theigneous rocks in Japan (VI) – granitic rocks,
summary. J. Japan Assoc. Mineral. Petrol. Econ.
Geol., 57, 177 87 (Japanes e with English abstract ).
Koch, B. and MacGillavry, C.H. (1962 ) X-ray
absorptio n. Pp. 157 60 in: Internationa l Tables for
X-ray Crystallography , Vol. III (C.H. MacGillavry
and G.D. Rieck, editors ). The Kynoch Press, UK.
Lentz, D.R. and Fowler, A.D. (1992) A dynamic modelfor graphic quartz-feldspa r intergrowth s in granitic
pegmatites in the southwester n Grenville Province.
Canad. Mineral., 30, 571 85.
Marschallinger , R. (1998) Correction of geometricerrors associated with the 3-D reconstructio n ofgeological materials by precision serial lapping.
Mineral. Mag., 62, 783 92.
Matsubara, H. (1956) Reports on the pegmatites in
Ishikawa district, Fukushima Prefecture – One of thefundamenta l studies on uranium-thoriu m resource s
in Japan. Bull. Geol. Surv. Japan, 7, 335 48(Japanese with English abstract ).
McCullough, E.C. (1975) Photon attenuatio n incomputed tomography. Medical Physics, 2, 307 20.
Morishita, R. (1998) Statistical properties of ideal rocktextures: relationshi p between crystal size distribution and spatial correlatio n of minerals. Math. Geol.,30, 409 34.
Morishita, R. and Obata, M. (1995 ) A new statistica ldescriptio n of the spatial distribution of minerals inrocks. J. Geol., 103, 232 40.
Nakano, T. and Fujii, N. (1989a ) Softwares for digital
image processing: (1) painting and virtual screen
management. Geoinformatics, 14A, 93 107(Japanese with English abstract ).
Nakano, T. and Fujii, N. (1989b ) The multiphase grain
control percolation : its implication for a partiallymolten rock. J. Geophys. Res., 94, 15653 61.
Nakano, T. and Fujii, N. (1991 ) Softwares for digital
image processing : (3) cluster labelling and perimeter
extractio n. Geoinformatics, 2, 23 44 (Japanese withEnglish abstract).
Nakano, T., Nishizawa, O., Masuda, K., Inazumi, T. andKasama, S. (1992) Three-dimensiona l distribution of
cracks and minerals in a rock obtained by X-ray CT.
Geotomography , SEGJ Int. Pub., 2, 361 71.
Nakano, T., Nakamura, K., Someya, T. and Ohtsuka, H.
(1997 ) Observatio n of 3-dimensiona l internal structure of rock using X-ray CT: (1) density calibratio nof CT value. Geoinformatics, 8, 239 55 (Japanes e
with English abstract ).
Nakashima, Y., Hirai, H., Koishikawa, A. and Ohtani, T.
(1997 ) Three-dimensiona l imaging of arrays of  uid
inclusions in  uorite by high-resolutio n X-ray CT.
Neues Jahrb. Mineral. Mh., 559 68.
Philpotts, A.R., Shi, J. and Brustman, C. (1998) Role ofplagioclase crystal chains in the differentiatio n ofpartly crystallized basaltic magma. Nature, 395,
343 6.
Raynaud, S., Fabre, D., Mazerolle, F., Geraud, Y. andLatiere, H. (1989) Analysis of the internal structureof rocks and characterizatio n of mechanical deformation by a non-destructiv e method: X-ray tomodensitometr y. Tectonophysics, 159, 149 59.
Simigian, S. and Starkey, J. (1986) Automated grain
shape analysis. J. Struct. Geol., 8, 589 92.
Simpson, D.R. (1962) Graphic granite from the Ramonapegmatite district, California. Amer. Mineral., 47,1123 38.
Stauffer, D. (1985 ) Introductio n to Percolation Theory.
Taylor & Francis, Philadelphia , PA, USA.

955

S. IKEDA ETAL.

Stel, H. (1992) Diagnostic microstructure s for primaryand deformationa l quartz rods in graphic granite.
Amer. Mineral., 77, 329 35.
Wahlstrom, E.E. (1939) Graphic granite. Amer.

Mineral., 24, 681 98.
[Manuscript received 16 February 1999:
revised 7 July 1999]

Appendix
Ellipsoid-fitting to a three-dimensional object
We carried out ellipsoid- tting to the 3-dimensional object (i.e. high-value region of a 2-point
correlation function) and obtained the directionsof principal axes and the ratios between the radiiof the object (ellipsoid). This  tting was
performed on the assumption that the centre ofbalance of the object corresponded with the centreof the ellipsoid. Firstly, the direction the c-axis of
the ellipsoid was determined by minimizing the
2nd moment around c-axis. Next, the a-axis wasdetermined by maximizing the 2nd moment
around a-axis. The detailed calculation process
is as follows.
Expression of an ellipsoid
Here we show the fundamental equations which
express an ellipsoid. Internal space within an
ellipsoid was expressed as:
a
A

b
B

2

2

c
C

2

41

A1

dadbdc

4
pABC
3

M0

a

b

1
V A2
5

B2

2

c2

Mb

c2

a2

Mc

2

2

C2

C2

A2

A2

B2

A2

x
y
z

b

A4

ea

eb

c dadbdc
A3

ec

The 2nd moments around the a, b and c-axes
are expressed as:

for

956

aea

beb

cec

A5

where

2

C2

a
B2

dadbdcTransformation between the coordinates
Here we show the way to transform a coordinate
belonging to a coordinate system into a coordinate
belonging to other coordinate systems. If the
abc-coordinate system is not the same as thexyz-coordinate system, i.e. if the angles l, f and y
in Fig. A1b =0, we have to transform thecoordinates to express the ellipsoid in thexyz-coordinate system. The transformation of thecoordinate belonging to the abc-coordinate
system can be transformed into the coordinate
belonging to the xyz-coordinate system by theequation

The 2nd moment around the origin of thecoordinate axes (a, b and c-axes) is expressed as:
2

b2

1
V
5

where a, b and c indicate the values of acoordinate of the point in abc-coordinate
system, and A, B and C indicate the radii on thea, b and c-axes, respectively (Fig. Ala). In thisappendix we de ne A4B4C. The volume of theellipsoid is expressed as:VMasin lcos y cos lsin f sin ycos lcos y sin lsin f sin ycos f cos ysin lsin y cos lsin f cos y
cos lsin y sin lsin f cos y
cos f cos y
cos lcos fandsin lcos f
sin f
p4 l4p,

1p4f4 1p,2
2

A6

p4y4p.

SHAPE OF CRYSTALS IN GRAPHIC GRANITEab

FIG. A1. Coordinate system expressing an ellipsoid. The sign of the Z-axis is reversed. (a) An ellipsoid and its
principal axes. (b) The generalize d coordinate system of ellipsoid in the case where the abc-axes of the ellipsoid do
not correspon d to the XYZ-axes of the observatio n coordinat e system.

Image datasystem into the coordinate in abc-coordinate
systemImage data are constructed in the xyz-coordinate
system and the origin of the coordinate systemdoes not lie at the centre of the object. Thus, as asimple calculation, we moved the origin of thecoordinate system to the centre (of gravity) of theobject. The centre of gravity of the object isexpressed as a vector in the xyz-coordinate system:
N

1
N i 1

R0

Xi
Yi
Zi

xi
yi
zi

Xi
Yi
Zi

ci ec

(A9)

(A10)

V = ND
N

M0

A7

i 1

R0

bi eb

The volume of the object and the 2nd moment
around the origin of the coordinate axes are
expressed as;

where N is the number of voxels which belong tothe object, and X, Y and Z are the values of thecoordinate of the ith voxel. In the following
calculations, this point (coordinate) is regarded as
the centre of the xyz-coordinate system.
Furthermore, the centre point of the ellipsoidwas  xed at the origin of the xyz-coordinate
system. Thus we transformed all the coordinatesin the old xyz-coordinate system into the newcoordinates as follows:
ri

ai ea

ri

a2i

b2i

N

c2i D

2

ri D

L0 D

A11

i 1

where D indicates the volume of each voxel and
L0 is
N

L0
i 1

x2i

y2i

z2i

A12

From equations A3, A10 and A11, we canobtainA2 + B2 + C2 = 5M0 /V = 5L0 /N(A13)

Determination of the c-axis

A8

Using equation A6, we can transform thecoordinate of each voxel in xyz-coordinate

957

From the calculations above we established theessential values of the object. Here we determine
the c-axis of the object. Only the angles l and f
are involved in these calculations. The angle y is

S. IKEDA ETAL.

irrelevant. From equations A6, A8, A9, A10 and
A11, the 2nd moment around the c-axis isexpressed as:
N

Mc
i 1

a2i

b2i D
N

M0

c2i D

N

2

ri

c2i

i 1
N

i 1

A14

yi sin l cos fzi sin f 2 D

In the last section we determined the direction of
the c-axis (i.e. the angles l and f) and radius C.
In this section, we determine the angle y by
maximizing the 2nd moment around the a-axis.
Maximizing the 2nd moment around the a-axis(Ma ) is equivalent to minimizing the 2nd moment
around the b-axis (Mb ). Then we determine theradii B and C.
The 2nd moments around the a and b-axes are
expressed as: (see equations A18 to 21 below)i 1

We determined l and f which minimizes thevalue of Mc . From equations A4 and A13, we candetermine
A2 + B2 = 5Mc /V = 5Lc /N(A15)

Lc = Mc /D

(A16)

where

N

Ma
i 1

b2i

N

c2i D
N

2

ri

a2i D

i 1

2
ri ea D

M0

N

Mb
i 1

i 1

a2i D

ni sin y 2 D

hi cos y

A18

i 1

Lab y D
c2i

a2i D
N

M0

N

ri

2

b2i D

i 1

2
ri eb D

N

M0
i 1

N

M0

hi sin y

i 1

M0

N

M0

N

M0

i 1

M0

(A17)

Determination of the a and b-axes

i 1

xi cos lLc )/Nand then, we can determine the radius of C.

N

M0

C2 = 5(L0D

ri ec 2 D

M0

Finally, from equations A13 and A16, we candetermineb2i D

ni cos y 2 D

A19

i 1

1
pD
2

Lab y

where
N

Lab y

hi cos yni sin y 2

A20

i 1

and
hi
ni

xi sin l yi cos lxi cos l yi sin l sin f
i

1

N

958

zi cos f

A21

SHAPE OF CRYSTALS IN GRAPHIC GRANITEWe determined y (04y<p) which minimizes thevalue of La b . Then we calculated
La
Lb

Ma /D

L0

Mb /D

L0

Finally, from equations A13 and A16, we canobtain:

Lab y
Lab y

1
p
2

C2

5Ma /V

5La /N2

2

5Mb /V

5Lb /N

C

A

5 L0

La /N

2

5 L0

Lb /N

B

A22

A24

and then, we can determine the radii A and B.

From equations A4 and A13, we can obtain
B2

A2

A23

959