- Email: [email protected]
Spherical harmonics in texture analysis
Tectonophysics 370 (2003) 253 – 268 www.elsevier.com/locate/tecto
Spherical harmonics in texture analysis Helmut Schaeben *, K. Gerald van den Boogaart Geoscience Mathematics and Informatics, Freiberg University of Mining and Technology, Gustav-Zeuner-Str. 12, Freiberg D-09596, Germany Accepted 31 March 2003
Abstract The objective of this contribution is to emphasize the fundamental role of spherical harmonics in constructive approximation on the sphere in general and in texture analysis in particular. The specific purpose is to present some methods of texture analysis and pole-to-orientation probability density inversion in a unifying approach, i.e. to show that the classic harmonic method, the pole density component fit method initially introduced as a distinct alternative, and the spherical wavelet method for highresolution texture analysis share a common mathematical basis provided by spherical harmonics. Since pole probability density functions and orientation probability density functions are probability density functions defined on the sphere X3oR3 or hypersphere X4oR4, respectively, they belong at least to the space of measurable and integrable functions L1(Xd), d = 3, 4, respectively. Therefore, first a basic and simplified method to derive real symmetrized spherical harmonics with the mathematical property of providing a representation of rotations or orientations, respectively, is presented. Then, standard orientation or pole probability density functions, respectively, are introduced by summation processes of harmonic series expansions of L1(Xd) functions, thus avoiding resorting to intuition and heuristics. Eventually, it is shown how a rearrangement of the harmonics leads quite canonically to spherical wavelets, which provide a method for high-resolution texture analysis. This unified point of view clarifies how these methods, e.g. standard functions, apply to texture analysis of EBSD orientation measurements. D 2003 Elsevier B.V. All rights reserved. Keywords: Real-valued spherical harmonics; Real-valued harmonics for the rotation group; Harmonic series expansions; Spectral decomposition of the X-ray transform; Summability; Spherical singular integral; Spherical radial basis functions; Spherical wavelets; Methods of texture analysis
1. Introduction The subject of texture analysis is the experimental determination (collecting and processing data) and interpretation of the statistical distribution of orienta-
* Corresponding author. Fax: +49-3731-394-067. E-mail address: [email protected] (H. Schaeben).
tions of crystals within a specimen of polycrystalline materials, which could be metals, rocks or bio-materials. It is often referred to as analysis of crystallographic lattice preferred orientation (LPO), and the objective of the interpretation is to relate an observed pattern of preferred orientation to its generating processes or vice versa. The orientation of an individual crystal is assumed to be unique and given by the rotation g 1aSO(3) which brings a coordinate system Ks fixed to the specimen into coincidence with a coordinate
0040-1951/03/$ - see front matter D 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0040-1951(03)00190-2
254
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
system Kc fixed to the crystal, g 1: KsiKc. Thus, in the context of texture analysis, rotation always means passive or ‘‘frame’’ rotation. The coordinates of a unique direction represented by haX3oR3 with respect to the crystal coordinate system Kc (referred to a crystallographic direction) and by raX3 with respect to the specimen coordinate system Ks (referred to as specimen direction) are related to each other by rotations in the set of rotations Cðh; rÞ ¼ fqaSOð3Þ j h ¼ grg which constitutes a one-dimensional circle on X4oR4 with center O parametrized by h and r when rotations are represented as unit vectors (‘‘quaternions’’) on X4. More specifically, conventional texture analysis is the determination and interpretation of the orientation probability density function f of a polycrystalline specimen by volume without individual orientation measurements; it is emphasized that size, shape, and spatial location of the contributing grains are not considered. In X-ray diffraction experiments, the orientation probability density function f cannot be directly measured, but with a texture goniometer only pole probability density functions P˜(h, r) can be sampled, which represent the probability that an experimentally fixed crystal direction h or its antiparallel h statistically coincide with a specimen direction r. Then, a pole probability density function is basically the mean of two spherical X-ray transforms 1 ˜ Pðh; rÞ ¼ ððXf Þðh; rÞ þ ðXf Þðh; rÞÞ 2 of the orientation probability density function f, where the spherical X-ray transform is defined as mean of f along the one-dimensional circle C(h, r), i.e. ðXf Þðh; rÞ ¼
1 2p
Z
f ðgÞdg
ð1Þ
fgaSOð3Þjh¼grg
In mathematical tomography, the transform assigning to a function defined on a d-dimensional manifold its mean values with respect to the family of dVdimensional submanifolds with 1 V dVV d 1 is referred to as dV-plane transform. If dV= 1, it is explicitly called X-ray transform; if dV= d 1, it is explicitly called Radon transform. For d = 2, the two cases coincide. The start of spherical mathematical tomog-
raphy dates back to the pioneering papers by Funk (1913, 1916), while mathematical tomography in an Euclidean setting commenced later with the classic paper by Radon (1917). In this sense, the term ‘‘X-ray transform’’ applies to texture analysis and is used by us (cf. Cerejeiras et al., 2002), even though it does not refer to the actual radiation, which could be c, neutron, or synchrotron as well. The invariant Haar measure dg in Eq. (1) is uniquely defined by the postulate that the measurements should essentially be independent of the choices of the involved coordinate systems Ks, Kc, which in turn requires an invariance of the form ðX ½f ðBÞ Þðh; rÞ ¼ ðX ½f ðgc1 Bgs Þ Þðgc h; gs rÞ Corrected experimental X-ray, neutron or synchrotron diffraction intensity data are thought of as being discretely sampled with a texture goniometer from continuous pole probability density functions P˜(h, r), i.e. even probability density functions defined on the cross-product X3 X3 of two unit spheres in R3. With respect to the diffraction experiment, the feasible crystal directions are the normals of the crystallographic lattice planes. Several representations of pole probability density functions exist, among others (i) truncated series expansion into spherical harmonics, (ii) finite series expansion into different standard functions, (iii) series expansion into spherical wavelets. These representations lead to practical methods to resolve the inverse problem of texture goniometry to reconstruct a reasonable orientation probability density function f defined on SO(3), or equivalently on the hypersphere X4oR4, from the given corrected intensity data, which are mean values of f along onedimensional great circles of X4 parametrized by h and r. Each method exploits that the X-ray transform of the harmonics for the group SO(3) are the harmonics for the sphere X3.
2. Spherical harmonics To compute orthonormal real-valued and symmetrized harmonics is not only of special interest in texture analysis but of general mathematical interest (cf. Zheng and Doerschuk, 2000).
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
The purpose of this section is to clarify the simple and elegant structure of the harmonic functions for the rotation group SO(3) and their straightforward relation to the spherical harmonics for the sphere X3 mediated by the spherical X-ray transform (Eq. (1)). In contrast to more conventional presentations of this subject (cf. Bunge, 1982; Matthies et al., 1987), harmonics are introduced here directly as real functions. Introduced in this way, they apply most readily to texture analysis, and unnecessary efforts with complex functions are avoided (cf. Boogaart, 2001). A fundamental theorem of representation theory assures that a realvalued representation exists (cf. Curtis and Reiner, 1962). An orientation may given by three Euler angles providing the angles of an sequence of frame rotations about some axes of the coordinate system ordered according to some convention. These rotations are given by 1 0 1 0 0 C B C B C Rx ðuÞ ¼ B 0 cosðuÞ sinðuÞ C B A @ 0 sinðuÞ cosðuÞ
255
Thus, a rotation about the z-axis can be achieved by first rotating the z-axis with Rx( p/2) = Rtx (p/2) onto the y-axis, then applying the rotation Ry(u) by the given angle u about y, and eventually rotating the yaxis. Analogously with Qy ¼ Ry
p 2
;
Qz ¼ R z
p 2
it holds Rx ðuÞ ¼ Qtz Ry ðuÞQz ¼ Qy Rz ðuÞQty Ry ðuÞ ¼ Qy Rx ðuÞQty ¼ Qtx Rz ðuÞQx Rz ðuÞ ¼ Qty Rx ðuÞQy Then any rotation gaSO(3) may be decomposed into g ¼ Rz ðu1 ÞRx ðUÞRz ðu2 Þ ¼ Rz ðu1 ÞQy Rz ðUÞQty Rz ðu2 Þ ¼ Rz ðaÞRy ðbÞRz ðcÞ
0 B B Ry ðuÞ ¼ B B @
cosðuÞ
0
sinðuÞ
0
1
0
sinðuÞ 0
C C C C A
0 cosðuÞ
sinðuÞ
cosðuÞ
1
0
1
B C B C B Rz ðuÞ ¼ B sinðuÞ cosðuÞ 0 C C @ A 0 0 1 p p ¼ Rx Ry ðuÞRtx ¼ Qx Ry ðuÞQtx 2 2 with
Qx ¼ R x
0
1
B B ¼B B0 2 @ 0
p
0 0 1
0
1
C C 1C C A 0
¼ Rz ðaÞQtx Rz ðbÞQx Rz ðcÞ Harmonics for SO(3) can generally be introduced by application of group representation theory (cf. Curtis and Reiner, 1962), i.e. the theory of all possible, essentially different representations of elements of SO(3) by matrices of arbitrary order. Harmonics provide a set of representatives of all irreducible representations, which in turn are the characteristic representations. This set is not unique; an additional specification is required to impose uniqueness. Following the complex approach of Gel’fand et al. (1963) and Vilenkin (1968), the common application in texture analysis is to transform initial complex harmonics into real harmonics (cf. Bunge, 1969, 1982), which, however, are no longer harmonic representations of rotations, i.e. compared with the complex harmonics, they do not possess the single most important property of being a representation. For the proper transform preserving the representation property, the reader is referred to Boogaart (2001).
256
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
The real and complex harmonic functions for the rotation group have a very similar structure. For real harmonics, we start with the harmonic representation of a simple rotation about the z-axis S RS ðuÞ ¼ ðrmn Þm¼S ;...;S ;n¼S ;...;S 1 0 cosðS uÞ 0 sinðS uÞ C B C B C B O q C B C B C B C ¼B 0 1 0 C B C B C B C B q O C B A @ sinðS uÞ 0 cosðS uÞ
The same object looks a little less intuitive for complex harmonics like C RS C ðuÞ ¼ ðrmnS Þm¼S ;...S ;n¼S ;...;S 0 iS u e B B B O B B B ¼B 1 B B B B O B @
0
1 C C C C C C C C C C C C A
eiS u
0
Together with a rotation Q by p/2 about another axis, harmonics can be rewritten similarly to rotation matrices as TS ðgÞ ¼ ðTSmn ðgÞÞm¼S ;...;S ;n¼S ;...;S ¼ RS ðu1 ÞQSx RS ðUÞQSx t RS ðu2 Þ
and QSy wRSy
QSx wRSx
2 p 2
As with rotation matrices, harmonics in matrix notation possess a distinguished property also known as the representation property (cf. Gel’fand et al., 1963; Vilenkin, 1968) TS ðggVÞ ¼ TS ðgÞTS ðgVÞ In texture analysis, the orientation probability density function f is expanded into a Fourier series of harmonic functions. With the trace of a matrix trA= P trðaij Þij ¼ i aii , this expansion can be written in condensed matrix notation as
f ðgÞ ¼
l X S S X X
CSmn TSmn ðgÞ ¼
S ¼0 m¼S n¼S
l X
trðCtS TS ðgÞÞ
S ¼0
with matrices CS of CSmn coefficients CS ðCSmn Þm¼S ;...;S ;n¼S ;...;S The representation property of the harmonics allows to calculate the harmonic expansion of the orientation probability density function with respect to a different sample or crystal coordinate system by simple matrix calculations as follows
f˜ ðgÞwf ðggVÞ ¼
l X
trðCtS TS ðgÞTS ðgVÞÞ
S ¼0
¼ RS ðaÞQSy RS ðbÞQSy t RS ðcÞ ¼
with
p
l X
trðTS ðgVÞCtS TS ðgÞÞ
S ¼0
RSz ðuÞwTS
ðRz ðuÞÞ ¼ R ðuÞ S
¼
trððCS TtS ðgVÞÞt TSmn ðgÞÞ |fflfflfflfflffl{zfflfflfflfflffl} S ¼0 CSmn of f ðgVÞ
RSx ðuÞwTS ðRx ðuÞÞ ¼ QSy RS ðuÞQSy t RSy ðuÞwTS ðRy ðuÞÞ ¼ QSx t RS ðuÞQSy
l X
¼
l X S ¼0
˜ t TS ðgÞÞ trðC S
ð2Þ
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
Let us now focus on X-ray transforms and pole probability density functions, respectively. A crystallographic direction h may be presented by two Euler angles u, U corresponding to the rotations Rz(u) and Rx(U) which rotate the z-axis of Kc onto h.
Analogously, a specimen direction r can be defined by two Euler angles uV, UV, providing the rotations Rz(uV) and Rx(U)Vwhich rotate the z-axis of Ks onto r. r ¼ Rz ðuVÞRx ðUVÞez Calculating (Xf)(h, r) is most simple for h = r = ez, where ez denotes the unique vector representing the direction of the coinciding z-axis of the coordinate systems. Then 1 2p
Z
˜ S e0 C˜ S00 ¼ e0 C ¼ tret0 RSz ðuVÞRSx ðUVÞCS ðRSz ðuÞRSx ðUÞÞt e0 ¼ trCtS ðRSx ðUVÞÞt ðRSz ðuVÞÞt e0 et0 Rz ðuÞRz ðUÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffi t 4p 4p 2S þ1 YS
h ¼ Rz ðuÞRx ðUÞez
f ðRz ð/ÞÞd/ Z 2p l X tr CS RS ð/Þd/ ¼ S ¼0
l X S 1 X þ CSmm cosðm/Þd/ ¼ 2p S ¼0 m¼S Z 2p Z 2p CSmm sinðm/Þd/ þCSmm cosðm/Þd/ 0 0 |fflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} l 1 X ¼ C 00 2p S ¼0 S
0
l X
4p trCSt YtS ðhÞYS ðrÞ 2S þ 1 S ¼0
l X S S X X CSmn YS ;m ðhÞYS ;n ðrÞ 2S þ 1 S ¼0 m¼S n¼S
ð3Þ
is retrieved; for a complementary interpretation and application the reader is referred to (Boogaart and Schaeben, 1999). Harmonics, being real-valued right from their introduction, may then be symmetrized in much the same way as their complex-valued harmonic relatives.
0 when m p 0
Now using u, U, uV, UV to rotate the coordinate systems of the specimen and the crystal, respectively, such that h = ezaKc and r = ezaKs, the orientation probability density function in the rotated space is given by f˜ ðgÞ ¼ f ððRz ðuVÞRx ðUVÞÞt gRz ðuÞRx ðUÞÞ According to formula (2), the C-coefficients of ˜f are given by ˜ S ¼ RS ðuVÞRS ðUVÞCS ðRS ðuÞRS ðUÞÞt CS RS C z x z x ðuÞQSy RS ðUÞQSy
ðhÞ
which could be used as definition of spherical harmonics, the well-known spectral decomposition of the X-ray transform (Bunge, 1969, 1982; Matthies et al., 1987)
¼ 4p
Z 2p l X S 1 X ¼ CSmm sinðm/Þ 2p S ¼0 m¼S 0
2S þ1 YS
rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2S þ 1 S YS ðhÞ ¼ ðRx ðUVÞÞt ðRSz ðuVÞÞt e0 4p
ðXf Þðh; rÞ ¼
0
ðrÞ
Thus, setting
2p 0
257
3. Standard functions in texture analysis According to the intuitive introduction of standard orientation probability density functions into texture analysis by the Rossendorf school (cf. Matthies, 1980; Matthies et al., 1987, 1988, 1990; Eschner, 1993, 1994), it should satisfy the requirements (i) to be specified by a small number of parameters, (ii) mathematically and numerically tractable, (iii) generally appealing from the point of physical interpretation, and (iv) possess analytical expressions for corresponding pole probability density functions and all normalization constants, thus avoiding truncation errors. This informal characterization includes that their harmonic series expansion
258
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
is known. Later, all Rossendorf standard orientation probability density functions widely applied in pole probability density (or ambitiously: texture) component fit methods have been shown to actually reduce to special cases of the von Mises– Fisher matrix distribution on SO(3) or, equivalently, to the Bingham distribution of axes on X+4oR4, particularly to the symmetrically unimodal and the symmetrically circular ‘‘fibre’’ case (Schaeben, 1996a). Thus, the von Mises– Fisher matrix and the Bingham distribution, respectively, provide a general mathematical model orientation probability density function. Similarly, other standard functions like the spherical Brownian (Savelova, 1984; Nikolayev and Savyolova, 1997; Schaeben and Nikolayev, 1998) and the spherical de la Valle´e Poussin (Schaeben, 1997, 1999) have been introduced, which mark with respect to harmonic series expansion and analytical representation of the limiting function the ‘‘extreme ends’’ on the scale of standard functions: There is no analytical representation for the Brownian, while the harmonic series expansion of the de la Valle´e Poussin is finite. Recently, another standard probability density function proportional to exp(j tan2 (x/2)) has been introduced (Ivanova and Nikolayev, 2001). Complementary to their initial introduction, standard functions are derived here by summability or summation processes, respectively, of harmonic series expansion revealing their intricate relation with them and avoiding resorting to intuition and heuristics. 3.1. Notation Basically following Freeden et al. (1998) and Hochstadt (1986), some basic notation is recalled and fundamental results are summarized. The unit sphere in d dimensions is denoted Xd and defined as Xd ¼ fraRd j ArA ¼ 1g If e1, . . ., ed are orthonormal vectors in Rd, and raXd, then pffiffiffiffiffiffiffiffiffiffiffiffi r ¼ ted þ 1 t 2 rd1
where t = rted and rd 1 is a unit vector in the space spanned by e1, . . ., ed 1, and dxd ¼ ð1 t 2 Þðd3Þ=2 dtdxd1 (cf. Watson, 1983). Thus NXd N ¼
Z Xd
dxd ¼
2pd=2 Cðd=2Þ
e.g. NX2N = 2p, NX3N = 4p, NX4N = 2p2. The class L p (Xd) consists of scalar functions f: d X iR1, which are measurable and for which Nf NLp ðXd Þ ¼
Z Xd
Af ðxÞAp dxp ðxÞ
1=p
Thus, for d = 3, p = 1, it holds pN1N ffiffiffiffiffiffi L1(X3) = 4p, for d = 3, p = 2 it holds N1NL2 ðX3 Þ ¼ 4p. Since our main interest is in results for X3oR3, which apply to experimentally accessible pole probability density functions, omitting the superscript and writing in short X always refers to X3. Most results analogously generalize to d >3. Real spherical harmonics of different orders are orthogonal in the sense of the L2 inner product Z ðYn ; Ym ÞL2 ðX3 Þ ¼ Yn ðrÞYm ðrÞdx3 ðrÞ ¼ 0 X3
There are 2n + 1 linearly independent harmonics of order n for X3; therefore, let { Yn, 1, . . ., Yn, 2n + 1} denote an orthonormal sequence of harmonics of order n with respect to the inner product of L2 (X3), and let Yn,j, j = 1, . . ., 2n + 1 denote one of its members. Then it holds n X j¼n
Yn;j ðr1 ÞYn;j ðr2 Þ ¼
2n þ 1 Pn ðr1 r2 Þ 4p
where Pn denotes the Legendre polynomial (cf. Freeden et al., 1998; p. 32 (2.3.12), p. 38). It can be written as 2 1 d ðt 2 1Þn ; Pn ðtÞ ¼ n 2 n! dt
ta½1; 1
(cf. Freeden et al., 1998; p. 41 (3.2.14)).
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
As defined above, a pole probability density function is actually a probability density function P˜(h, r) defined on X3 X3. Mathematically, haX3 and raX3 are any two-unit vectors in R3. With respect to the diffraction experiment, however, h refers to the unit normal direction of a crystallographic lattice plane and therefore cannot vary arbitrarily as can raX3. To emphasize the role of h as a parameter of the experiment, pole probability density functions used to be denoted P˜h(r) or P˜h(r) to indicate symmetry of the crystallographic form h. In the following, harmonic analysis is applied to functions F(h, r) defined on X3 X3 read as a function of r parametrized by h. Then, any function F(h, r)aX (X3 X3), where X denotes either Lp, 1 V p < l, or C, can formally be associated with its Fourier series, i.e. Fðh; rÞf
l X n X
ðFðh;BÞ; Yn;j ðBÞÞYn;j ðrÞ
integrals, i.e. to convolutions of the function with a zonal kernel function approximating identity. The nth partial sum of the harmonic series expansion is defined as Sn ðFðh; BÞ; rÞ ¼
with constants (independent of r)
^ Fn;j ðhÞYn;j ðrÞ
S ¼0 j¼S
¼
n X S X
ðFðh; BÞ; Yn;j ðBÞÞYn;j ðrÞ
S ¼0 j¼S
¼
Fðh; BÞ;
n X S X
! Yn;j ðBÞYn;j ðrÞ
S ¼0 j¼S n X 2S þ 1 PS ðr BÞ Fðh; BÞ; 4p S ¼0
¼
Gn ðtÞ ¼
^ ðhÞ ¼ ðFðh; BÞ; Yn;j ðBÞÞ Fn;j
! ð4Þ
n X ð2S þ 1ÞPS ðtÞ S ¼0
referred to as its spherical Fourier-transform or Fourier-coefficients. Thus, Fðh; rÞf
n X S X
Setting
n¼0 j¼n
l X n X
259
¼ ðn þ 1Þ
Pnþ1 ðtÞ Pn ðtÞ t1
ð5Þ
which may be addressed as the spherical Dirichlet kernel, it is
^ Fn;j ðhÞYn;j ðrÞ
n¼0 j¼n 2
3
It can be shown that any function F(h, B)aL (X ) may be approximated in the L2 (X3) sense, i.e. in the norm of L2 (X3), by finite truncations of its Fourier orthogonal series expansion in terms of any L2 (X3) orthonormal system of spherical harmonics Yn,j. It should be noted that the above statement does not hold true for every continuous F(h, B)aC (X3) or F(h, B)aLp (X3) with pa[1, 4/3]v[4, l) (Freeden et al., 1998, p. 58). Thus, finite truncations of Fourier series expansions by means of L2 (X3)-orthonormal spherical harmonics do not generally form an ‘‘approximation process’’ on the sphere. For convergence in C(X3) norm it is required that F is Lipschitzcontinuous (Freeden et al., 1998, p. 59). Instead of imposing conditions on the function to be approximated, it may be more efficient to augment the Fourier series with ‘‘convergence factor’’ to ensure convergence. As it will be shown, the augmented series expansions correspond to spherical singular
1 Sn ðFðh; BÞ; rÞ ¼ 4p
Z X3
Fðh; rVÞGn ðr rVÞdx3 ðrVÞ
Of mathematical interest is the existence and uniqueness of the limit lim Sn ðFðh; BÞ; rÞ
n!l
As already mentioned above, the series of partial sums does not generally converge for FaC(X3) or FaL1 (X3) in the sense of the corresponding norm (cf. Freeden et al., 1998). Therefore, the linear space of all absolutely convergent sequences {hq (n, j)} of real numbers hq (n, j), (n, j)aJ ={(n, j) | n = 0, 1, . . ., j = 1, . . ., 2n + 1} for q in some parameter set IA is introduced. The sequence {hq (n, j)} furnished with appropriate properties is referred to as convergence factor eventually constituting a convergent summation process of the ‘‘augmented’’ harmonic series expansion as follows.
260
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
the set of functions g: [ 1,1]iR1, and endowed with norms NgAX ½1;1 ¼ Nqðe BÞNXðXd Þ
Instead of Sn( F(h, B); r), the focus is on Uq ðFðh; BÞ; rÞ ¼
l X n X
^ hq ðn; jÞFn;j ðhÞYn;j ðrÞ
n¼0 j¼n
¼
l X n X
¼
hq ðn; jÞðFðhÞ; BÞ; Yn;j ðBÞYn;j ðrÞ
¼
Fðh;BÞ;
n X
! hq ðn; jÞYn;j ðBÞYn;j ðrÞ
ð6Þ
n¼0 j¼n
which leads to different summation processes for different choices of hq (n, j), qaIA. If hq (n, j) = hq (n) for all j = 1, . . ., 2n + 1, then Uq ( F(h, B); r) simplifies to the special case Uq ðFðh; BÞ; rÞ ¼
l X n X
^ hq ðnÞFn;j ðhÞYn;j ðrÞ
¼
Fðh; BÞ
n¼0
!
Uq ðFðh; BÞ; rÞ ¼
1 4p
hq ðnÞPn ðr BÞ
4p Z
ð7Þ
Fðh; rVÞ
X3 l X
p
2p
1=p
AgðtÞA dt
they form subspaces of X (Xd) where X (Xd) denotes either Lp (Xd) or C (Xd), and X [ 1, 1] analogously. For zonal functions, their harmonic series expansions reduce to an expansion in terms of Legendre polynomials Pn,d for Rd with corresponding Fourier– Legendre coefficients. The spherical convolution of a zonal function gaL1 [ 1, 1] and a function faX (Xd) is defined by Z ðg*f ÞðyÞ ¼ gðy xÞf ðxÞdxd ðxÞ Xd
n¼0 j¼n l X 2n þ 1
1 1
n¼0 j¼n l X
Z
!
ð2n þ 1Þhq ðnÞPn ðr rVÞ
n¼0
dxðrVÞ Z 1 ¼ Fðh; rVÞKq ðr rVÞdx3 ðrVÞ 4p X3 ð8Þ If more specifically IA = IN, then hq (n, j) = hm (n), and if additionally 8 < 1 if nVm hm ðnÞ ¼ : 0 otherwise then the summation process Um( F(h, B); r) degenerates to the partial sums Sm ( F(h, B); r). 3.2. Summation processes and singular integrals on the sphere Any scalar function gy: XdiR1, xigy(x) = g(y x), x, yaXd is called zonal (synonymously central, rotationally invariant). Zonal functions are isomorphic to
A convolution has the nicer property of the corresponding property of g and f, respectively. In particular, the Fourier coefficients of the convolution are the product of the Fourier( – Legendre) coefficients of the functions being convoluted. A set of functions {vq | qaIA}oL1 [ 1, 1] with q in some parameter set IA is called a kernel if ðvq Þ^ ð0Þ ¼ 2p
Z ½1;1
vq ðtÞdt ¼ 1
holds. The convolution of a function faX (Xd) with a kernel vqaL1 [ 1, 1] is called a spherical singular integral vq ðf ; xÞ ¼ ðvq *f ÞðxÞ ¼
Z Xd
vðx yÞf ðyÞdxd ðyÞ
If the kernel function vaX [ 1, 1] is an approximate identity, i.e. if it is uniformly bounded for all qaIA, and if it possess the spatial localization property Z
d
lim
q!q0
1
vq ðtÞdt ¼ 0;
dað1; 1Þ
then the singular integral converges in X (Xd)-norm towards f for q ! q0. A nonnegative kernel is an approximate identity if lim ðvq Þ^ ð0Þ ¼ 1
q!q0
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
With respect to texture analysis, it is well known that if the orientation probability density function f is zonal with respect to g0aSO(3), then its X-ray transform is zonal, too, and it holds ðXf Þðh; rÞ ¼ ðXf Þðh g0 rÞ ¼ ðXf Þðg01 h rÞ Correspondingly, the familiar form of the harmonic series expansion of an orientation probability density function f ðgÞf
l X S S X X S ¼0 n¼S m¼S
for any zonal function f defined on X4, where the last equation is an application of the identity CS(1) (cos w) = US (cos w) with the Chebychev polynomials of the second kind US (cos w). Further, following an argument introduced by Arnold (1941), the corresponding series expansion for an even function f is 1 ½ f ðx; jÞ þ f ðp x; jÞ
2 l 1 X ¼ 2 ðS þ 1Þ2 f ^ ðS ÞPS ;4 ðcosxÞ 2p S ¼0ð2Þ
f ðx; jÞ ¼
4p C mn T mn ðgÞ 2S þ 1 S S
l 1 X ð2S þ 1Þ2 f ^ ð2S ÞP2S ;4 ðcosxÞ 2p2 S ¼0 l 1 X sin½ð2S þ 1Þx
¼ 2 ð2S þ 1Þf ^ ð2S Þ 2p S ¼0 sinx
¼
and its X-ray transform ðXf Þðh; rÞf
l X S S X X
CSmn YS ;m ðhÞYS ;n ðrÞ
ð12Þ
S ¼0 n¼S m¼S
is turned into f ðg; g0 ; jÞf
261
with the even-order Chebychev polynomials of second kind U2n (x/2)=((sin (2n + 1)x/2)/(sin x/2)). Formally, a spherical singular integral may be associated with its harmonic series
l X
2S þ 1 ^ f ð2S Þ 2p2 S ¼0
sin½ð2S þ 1Þxðgg01 Þ=2
sinðxðgg01 Þ=2Þ
ð9Þ vq ðF; rÞf
l X n X
ðvq ðF; BÞ; Yn;j ðBÞÞYn;j ðrÞ
n¼0 j¼n
and l X 2S þ 1 ðXf Þ^ ðX ½f ðB; g0 ; jÞ Þðh; rÞf 4p S ¼0
ðS ÞPS ðh g0 rÞ ¼ ðXf Þðh g0 r; jÞ
¼
l X n X
ððvq *FÞðBÞ; Yn;j ðBÞÞYn;j ðrÞ
n¼0 j¼n
ð10Þ
¼
l X n X
^ v^q ðnÞFn;j ðhÞYn;j ðrÞ
ð13Þ
n¼0 j¼n
(cf. Matthies et al., 1987; Schaeben, 1996b). Eq. (9) follows from f ðx; jÞ ¼
l 1 X ðS þ 1Þ2 f ^ ðS ÞPS ;4 ðcosxÞ 2 2p S ¼0
¼
l 1 X ð1Þ ðS þ 1Þf ^ ðS ÞCS ðcosxÞ 2p2 S ¼0
¼
l 1 X ðS þ 1Þf ^ ðS ÞUS ðcosxÞ 2p2 S ¼0
¼
with constants (independent of r)
1 2p2
l X
ðS þ 1Þf ^ ðS Þ
S ¼0
sin½ðS þ 1Þx
sinx
v^q ðnÞ ¼ 2p
Z
1
1
vq ðtÞPn ðtÞdt ¼ ðG; Pn Þ
and ^ Fn;j ðhÞ ¼ ðFðh; BÞ; Yn;j ðBÞÞ
ð11Þ
referred to as Fourier( –Legendre) coefficients. Comparing Eqs. (7) and (13), at least the formal analogy is obvious that hq (n, k) is a convergence factor and provides a summation process if it is equal
262
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
to the Fourier – Legendre coefficients of a (zonal) approximate identity vq (t). For detailed results and proofs the reader is referred to Dunkl (1966), Berens et al. (1968), Butzer and Nessel (1971), Butzer et al. (1979), and Freeden et al. (1998). 3.3. Examples For reasons of notational simplicity and practical applications, we confine ourselves again to the case X3, i.e. to X-ray transforms. The most simple choice for a spherical kernel is the Dirichlet kernel (Eq. (5)). Choosing more sophisticatedly as spherical kernel the (i) spherical Brownian, (ii) spherical de la Valle´ e Poussin, (iii) spherical Cauchy, (iv) von Mises – Fisher and (v) X-ray transformed von Mises-Fisher probability density function l 1 X ð2S þ 1Þexp½S ðS þ 1Þq PS ðtÞ (i) vq ðtÞ ¼ 4p S ¼0
l X
1 I S ðqÞ I S þ1 ðqÞ PS ðtÞ ð2S þ 1Þ 4p S ¼0 2B 32 ; q þ 12 1 1þt q ¼ ¼ vn ðtÞ Bð1; q þ 1Þ 2
(iii) vq ðtÞ ¼
l 1 X ð2S þ 1Þq2S PS ðtÞ 4p S ¼0 2
1 1q ¼ pq ðtÞ 4p ð1 þ q2 2qtÞ3=2 l I S þ1=2 ðqÞ 1 X PS ðtÞ ð2S þ 1Þ (iv) vq ðtÞ ¼ I 1=2 ðqÞ 4p S ¼0 ¼
¼
3.3.1. Singular integral of Gauss –Weierstrass, spherical Brownian standard probability density function Choosing hq ðn; jÞ ¼ hq ðnÞ ¼ expðnðn þ 1ÞqÞ with qaIA=(0, l), q0 = 0, it follows that l 1 X ð2n þ 1Þhq ðnÞPn ðtÞ 4p n¼0 ¼
l 1 X ð2n þ 1Þexpðnðn þ 1ÞqÞPn ðtÞ 4p n¼0
¼ bðt; qÞ
¼ bðt; qÞ (ii) vq ðtÞ ¼
leads to the singular integrals of (i) Gauss– Weierstrass, (ii) de la Valle´e Poussin, (iii) Abel – Poisson, (iv) von Mises– Fisher and (v) X-ray transformed von Mises– Fisher. All of them provide summation processes and thus accomplish convergence in the X (Xp) norm. In case of C (Xd), this is uniform convergence which in turn implies pointwise convergence.
q expðqtÞ ¼ mðt; qÞ 4psinhq
which is the spherical Brownian distribution for X3 (cf. Mardia and Jupp, 2000). Moreover, it is the X-ray transform of the spherical Brownian for X4 (cf. Schaeben, 1996b). Then Z Wq ðFðh; BÞ; rÞ ¼ Fðh; rVÞbðr rV; qÞdxðrVÞ ð15Þ X3
is the spherical singular integral of Gauss– Weierstrass (cf. Freeden et al., 1998; Schaeben, 1996b). It converges for q ! 0 to F(h, r) in the corresponding norm. A texture component fit method using spherically Brownian distributed components evolves as follows. In case of a finite data set F(h, rk) according to some spherical sampling scheme rkaSh adjusted to q wq ðFðh; BÞ; rÞ ¼
(v) vq ðtÞ ¼ ¼
l 1 X I S ðqÞ I S þ1 ðqÞ PS t ð2S þ 1Þ 4p S ¼0 I 0 ðqÞ I 1 ðqÞ
q q 1 I0 ð1 þ tÞ exp ðt 1Þ 4psinhq 2 2
¼ ðXM Þðt; qÞ
ð14Þ
NS hN X
ak ðhÞbðr rk ; qÞ
k¼1
provides a practical approximation of F(h, r) in the sense that NS hN X k¼1
ak ðhÞbðr rk ; qÞcFðh; rk Þ
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
The iterative fit of ak (h) may start with a(0) k (h) = F(h, rk). In practice, the series expansion may be truncated by the user to a rough approximation by K major texture components with q = q(k), k = 1, . . ., K uðFðh; BÞ; rÞ ¼
K X
ak ðhÞbðr rk ; qðkÞÞ
where a small number K bNShN denotes the total number of spherically Brownian distributed components to be considered. Usually, these major components are placed in some rk where F(h, r) displays strong local maximum. Strictly speaking, this usercontrolled method turns texture approximation by finite series expansion into texture modeling by a very small number of texture components, which may be very successful and helpful in practice when the pattern of lattice preferred orientation is known or sufficiently simple. 3.3.2. Singular integral of de la Valle´e Poussin, spherical de la Valle´e Poussin standard probability density function Choosing 8 ðjÞ > < I n ðjÞI 3 nþ1 1 if nVj 2B ; j þ hq ðn; jÞ ¼ hj ðnÞ ¼ 2 2 > : 0 if n > j where
0 1 2S S X B C I S ðjÞ ¼ ð1Þk @ A k¼0 2k 1 1 B k þ ;j þ S k þ ; 2 2
Z X3
Fðh; rVÞvn ððr rV; qÞdxðrVÞÞ ð17Þ
is the spherical singular integral of de la Valle´e Poussin (cf. Schaeben, 1996b). It converges for n ! l to F(h, r) in the corresponding norm. A texture component fit method using spherically de la Valle´e Poussin distributed components evolves analogously to the case of Gauss– Weierstrass. Since the infinite series representing the kernel reduces to a finite series, truncation is not required in practical numerical applications involving CS coefficients, and errors due to truncation are avoided. This property essentially distinguishes the de la Valle´e Poussin standard orientation probability density function from other known standard functions in an advantageous way, as it provides a computationally efficient model. 3.3.3. Singular integral of Abel – Poisson, spherical Cauchy- or Lorentz-type standard probability density function Choosing hq ðn; jÞ ¼ hh ðnÞ ¼ hn with ha[0,1), h0 = 1, it follows that l 1 X ð2n þ 1Þhh ðnÞPn ðtÞ 4p n¼0
S ¼ 0; 1; . . . with jaIA = N, j0 = l, it follows that n 1 X ð2S þ 1Þhj ðS ÞPS ðtÞ 4p S ¼0 n 1 X I S ðjÞ I S þ1 ðjÞ PS ðtÞ ð2S þ 1Þ ¼ 4p S ¼0 2B 32 ; j þ 12 1 1þt j ¼ ¼ vn ðtÞ Bð1; j þ 1Þ 2
which is the spherical de la Valle´e Poussin kernel for X3 and the X-ray transform of the spherical de la Valle´e Poussin kernel for X4 (cf. Schaeben, 1997, 1999). Then Vn ðFðh; BÞ; rÞ ¼
k¼1
263
ð16Þ
¼
l 1 X ð2n þ 1Þhn Pn ðtÞ 4p n¼0
¼
1 1 h2 ¼ ph ðtÞ 4p ð1 þ h2 2htÞ3=2
ð18Þ
which is the spherical Cauchy distribution for X3 (cf. Mardia and Jupp, 2000) or Lorentz-type distribution (Matthies et al., 1987), respectively; it is also the X-
264
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
ray transform of the spherical Cauchy distribution for X4. Then Z Fðh; rVÞph ðr rVÞdxðrVÞ ð19Þ Uh ðFðh; BÞ; rÞ ¼ X3
is the spherical singular integral of Abel –Poisson (cf. Freeden et al., 1998; Schaeben, 1996b). A texture component fit method using spherically Cauchy distributed components evolves analogously to the case of Gauss – Weierstrass or de la Valle´e Poussin. Owing to the poor localization, i.e. the fat tails, of the Cauchy distribution, it is rarely used in practical texture component fit procedures. 3.3.4. Standard von Mises – Fisher probability density function Choosing Inþ1=2 ðjÞ I1=2 ðjÞ
hq ðn; jÞ ¼ hj ðnÞ ¼
with jaIA = R+, j0 = l, it follows that
¼
Fðh; rVÞMðr rV; jÞdxðrVÞ
l 1 X ð2n þ 1Þhq ðnÞPn ðtÞ 4p n¼0 l 1 X Il ðjÞ Ilþ1 ðjÞ Pn ðtÞ ð2n þ 1Þ 4p n¼0 I0 ðjÞ I1 ðjÞ j j ¼ CM ðjÞI0 ð1 þ tÞ exp ðt 1Þ 2 2
¼
ð20Þ
ð22Þ
Z
Fðh; rVÞXM ððr rV; jÞdxðrVÞÞ
X3
ð23Þ
which is the von Mises– Fisher distribution for X (cf. Mardia and Jupp, 2000). It is not the X-ray transform of the rotationally symmetric von Mises – Fisher matrix probability density function on SO(3) or equivalently the bimodal Bingham, i.e. Watson probability density function for X4. Then Uj ðFðh; BÞ; rÞ ¼
with jaIA = R+, j0 = l, it follows that
Uj ðFðh; BÞ; rÞ ¼
3
Z
In ðjÞ Inþ1 ðjÞ I0 ðjÞ I1 ðjÞ
hq ðn; jÞ ¼ hj ðnÞ ¼
which is the X-ray transform of the von Mises –Fisher distribution for X4. Then
l Inþ1=2 ðjÞ 1 X Pn ðtÞ ð2n þ 1Þ I1=2 ðjÞ 4p n¼0
j expðjtÞ ¼ mðt; jÞ 4psinhj
3.3.5. X-ray-transformed standard von Mises –Fisher probability density function Choosing
¼ ðXM Þðt; jÞ
l 1 X ð2n þ 1Þhj ðnÞPn ðtÞ 4p n¼0
¼
Fisher pole probability density function seems to be unknown, it is not applied in texture component fit methods.
ð21Þ
is a spherical singular integral that converges for j ! l to F(h, r) in the norm. A texture component fit method using projected von Mises– Fisher distributed pole probability density components evolves analogously to the previous cases and corresponds to von Mises– Fisher distributed orientation probability density components. 3.4. Remark
X3
which is the spherical singular with respect to the von Mises –Fisher kernel (Schaeben, 1996b). Since the form of the orientation probability density function which is projected on the von Mises –
This first class of radial spherical basis functions (Gauss – Weiertrass, de la Valle´e –Poussin, von Mises– Fisher, Abel – Poisson) is characterized by the fact that all Legendre coefficients of the corresponding kernel are nonvanishing (Freeden et al., 1998). Note
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
that this is actually true for the de la Valle´e Poussin kernel for n ! l. The second class consists of locally supported kernel functions. Depending on the diameter of their support, some of their Legendre coefficients vanish (Freeden et al., 1998). Application of the Riemann– Lebesgue kernel, which is the indicator with respect to spherical caps was briefly considered in Schaeben (1996b); it is not considered here as its presentation is beyond the scope of this contribution. Unimodal radial basis functions can be transformed and adapted to represent multimodal fibre textures; formally, the required adaptation is accomplished by interpreting the angle x = arc cos t as orientation distance xf with respect to the fibre C(h0, r0) defined by h0, r0aX3.
As shown above, any function F(h, B) thought of as being parametrized by haX defined on X may be associated with its harmonic series expansion Fðh; rÞf
l X S X
FS^;j ðhÞYS ;j ðrÞ
S ¼0 j¼S
¼
l X S X
ðFðh; BÞ; YS ;j ðBÞÞYS ;j ðrÞ
ð24Þ
S ¼0 j¼S
The space of all spherical harmonics of degree S with S = 0, 1, 2, . . ., i.e. the restriction to the unit sphere of all homogeneous and harmonic polynomials in three variables, is denoted by Y S ; then dim Y S = 2S + 1. That any function in L2(X3) may be represented by its harmonic series expansion may be rewritten as
4. Spherical wavelets
L2
l
Another representation of spherical distributions in texture analysis is provided in terms of spherical wavelets (Schaeben et al., 2001). The main idea of wavelet analysis is to obtain a multiscale representation of the data or functions, which allows localization in space and frequency. Initially, a spherical probability density function is sampled on a coarse, ideally almost equidistributed, grid on the sphere. These coarse grid measurements are approximated by a spherical polynomial of low degree. This polynomial is clearly a sufficiently good approximation in regions of the sphere where the underlying function does not vary too much, i.e. where the underlying function consists of low frequencies only. Where the data show large variations, the initial approximation requires a local improvement, which is accomplished by a local refinement of the grid of measurements. The crucial point then is to construct a high-degree polynomial from the global coarse grid and the locally refined grid. This can be seen as adding adaptively and locally a wavelet part to the global approximation of low degree. Thus, wavelets provide a digital device to zoom into areas of special interest. Moreover, wavelets will provide the means to control the scanning process with a texture goniometer adaptively to a required local refinement of the spatial resolution.
265
L2 ðXÞ ¼
Y
S
S ¼0
Further, Pn (X) = PSn = 0 Y S with dim Pn (X)= (n + 1)2. Let Nk be a sequence of strictly monotone increasing positive integers. Then, scaling spaces Vk are introduced as polynomial spaces Vk ðX3 Þ ¼ PNk ðX3 Þ with dim Vk (X3)=(Nk + 1)2. Hence, there is a chain V0 ðX3 ÞoV1 ðX3 ÞoV2 ðX3 Þo . . . and it is sensible to define the corresponding orthogonal complements Nkþ1
Wk ðX3 Þ ¼ Vkþ1 ðX3 ÞOVk ðX3 Þ ¼
S ¼Nk þ1
YS
with dim Wk=(Nk + 1 Nk)(Nk + 1 + Nk + 2). Finally, l
L ðX3 Þ ¼ V0 ðX3 ÞP 2
W ðX Þ k
L2
3
k¼0
The classical approach of wavelet theory uses Nk = 2k. However, in case of spherical probability density functions, this would imply very large dimensions of the wavelet spaces. Therefore, the choice maybe more like Nk = 20 k.
266
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
with Ik={k=(k1,k2)TaN02: k1 < 2 Nk, k2 < 2 Nk} and
Using again 2X S þ1 j¼0
2S þ 1 PS ðr1 r2 Þ; YS ;j ðr1 ÞYS ;j ðr2 Þ ¼ 4p
rk;k
r1 ; r2 aX3
T
by using the Clenshaw –Curtis weights wkV(k). Analogously, with
with Legendre polynomials PS of degree S with PS (1) = 1 results in Fðh; rÞ ¼
pk1 pk2 pk1 pk2 pk2 cos cos ; cos sin ; sin 2Nk Nk 2Nk Nk Nk
l X
Fðh; BÞ;
S ¼0
2S þ 1 PS ðr BÞ 4p
SVk ðFðh; BÞ; rÞ ¼
Fðh; BÞ;
S ¼0
2S þ 1 PS ðr BÞ 4p
Nkþ1 X
ð2S þ 1ÞPS ðtÞ
S ¼Nk þ1
ð25Þ
Considering harmonics only up to order Nk, the approximation SVk ( F(h, B); r) with respect to Vk is given by Nk X
HV k ðtÞ ¼
and SWk ( F(h, B); r) c SWI k ( F(h, B); r) by numerical integration, it holds I SVI kþ1 ðFðh; BÞ; rÞ ¼ SVI k ðFðh; BÞ; rÞ þ SW ðFðh; BÞ; rÞ k
ð30Þ
Substituting SWI k ( F(h, B); r) by ð26Þ
I Sˆ W ðFðh; BÞ; rÞ ¼ k
X
wV k þ1 ðkÞLj Fðh; rV k þ1;k Þ
kaIkþ1
HV k ðr rkVþ1;k Þ
and the corresponding next wavelet part by
SWk ðFðh; BÞ; rÞ ¼
N k þ1 X
S ¼Nk þ1
2S þ 1 PS ðr BÞ Fðh; BÞ; 4p
ð27Þ Introducing the kernel GkVðtÞ ¼
ð31Þ
with Lk Fðh; rkþ1;k Þ 8 < Fðh; rkþ1;k Þ ¼ : I SV ðFðh; BÞ; rkþ1;k Þ
if FðhÞ is large ð32Þ otherwise
Nk X ð2S þ 1ÞPS ðtÞ
eventually yields
S ¼0
I Sˆ VI kþ1 ðFðh; BÞ; rÞ ¼ SVI k ðFðh; BÞ; rÞ þ Sˆ W ðFðh; BÞ; rÞ k
it holds 1 SVk ðFðh; BÞ; rÞ ¼ 4p
Z
ð33Þ Fðh; rVÞGkVðr rVÞdxðrVÞ
X3
ð28Þ The right-hand side can be numerically evaluated by a sampling theorem (Potts et al., 1996) such that SVk ( F(h, B); r) c SVIk ( F(h, B); r), i.e. SVI k ðFðh; BÞ; rÞ ¼
X kaIk
wkVðkÞFðh; rVk;k ÞGkVðr rVk ;k Þ
Summarizing, the proposed procedure is as follows. The function F(h, r) is known on a coarse r-grid and SVI ( F(h, B); r) is computed. By inspection or algorithms (Mhaskar et al., in press), the regions where F(h, r) and SVI ( F(h, B); r) have large absolute values and large variations, respectively, are detected. In the next step, the approximation SVI F is improved by adding the next wavelet part SWI F. Since F(h, r) should not be totally sampled on a refined r-grid and since it is known that SWI F is almost zero in regions
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
where F(h, r) does not oscillate too much, SVI F + SWI F is replaced by SVI F + SˆWI F. The definition of Lk shows that the function F(h, r) must be resampled for some additional points in the regions where large variations of F(h, r) are observed or expected. It can be shown that the spherical X-ray transform of the wavelet representation of f is the wavelet representation of the X-ray transform Af(h, B), and that the wavelet representation can be applied to solve the inverse problem to determine numerically the crystallographic orientation probability density function f from corrected experimental X-ray, neutron or synchrotron intensities measured in diffraction experiments with a texture goniometer.
267
analysis. Wavelet texture analysis combines the best elements of both methods, in particular the frequency localization of the harmonic and the spatial localization of the component fit method. It generalizes both with mathematical rigor and provides a wellfounded means for high-resolution texture analysis.
Acknowledgements The authors would like to thank the guest editors of this special issue of Tectonophysics for the invitation to contribute and the opportunity to bring a unified view of the fundamentals of mathematical texture analysis to the attention of a broader audience.
5. Conclusions References Harmonics for the group of rotations SO(3) and the sphere X3, respectively, provide a major mathematical prerequisite for texture analysis. Their application in texture analysis can be largely clarified and turned instructive if they are introduced in the way they are used-as real-valued functions avoiding the classical complex approach featured in textbooks of texture analysis. Their single most important property is that of being a group-theoretical representation. This property is lost if, e.g. the conversion of complex into real harmonics (Bunge, 1969, 1982) is done with less care than actually required. Convergence of the associated harmonic series expansions has never been an issue in texture analysis; it has been taken for granted. However, considerations of their summability provide a complementary introduction of standard distributions as required in texture component fitting. Initially introduced on heuristic grounds to deal with the ‘‘ghost problem’’ and the nonnegativity constraint, the focus here is on their common mathematical basis in terms of harmonics. The results presented here in the context of texture analysis are also fundamental for the analysis of individual crystallographic orientation measurements, e.g. the approximate identities are used in kernel density estimations of the orientation probability density function. The duality of summation processes and spherical singular integrals provides the key of an appropriate re-arrangement of the series leading to wavelet texture
Arnold, K.J., 1941. On spherical probability distributions. PhD thesis, Massachusetts Institute of Technology, MA, USA. Berens, H., Butzer, P.L., Pawelke, S., 1968. Limitierungsverfahren von Reihen mehrdimensionaler Kugelfunktionen und deren Saturationsverhalten. Publ. Res. Inst. Math. Sci., Kyoto Univ., Ser. A 4, 201 – 268. Boogaart, K.G.V.D., 2001. Statistics for individual crystallographic orientation measurements. PhD Thesis, Freiberg. Boogaart, K.G.v.d., Schaeben, H., 1999. Singular value decomposition of the odf – pdf projection operator. In: Szpunar, J.A. (Ed.), Proceedings of the 12th International Conference on Textures of Materials (ICOTOM12). NRC Research Press, Ottawa, Ont., Canada, pp. 1363 – 1368. Bunge, H.J., 1969. Mathematische Methoden der Texturanalyse. Akademie Verlag, Berlin. Bunge, H.J., 1982. Texture Analysis in Materials Science. Butterworths, London. Butzer, P.L., Nessel, R.J., 1971. Fourier Analysis and Approximation. Birkha¨user Verlag, Basel. Butzer, P.L., Stens, R.L., Wehrens, H., 1979. Approximation by algebraic convolution integrals. In: Prolla, J.B. (Ed.), Approximation Theory and Functional Analysis. North Holland Publishing, Amsterdam, pp. 71 – 120. Cerejeiras, P., Schaeben, H., Sommen, F., 2002. The spherical X-ray transform. Math. Methods Appl. Sci. 25, 1493 – 1507. Curtis, C.W., Reiner, J., 1962. Representation Theory of Finite Groups and Associative Algebras. Interscience, New York. Dunkl, C.F., 1966. Operators and harmonic analysis on the sphere. Trans. Am. Math. Soc. 125, 250 – 263. Eschner, T., 1993. Texture analysis by means of model functions. Textures Microstruct. 21, 139 – 146. Eschner, T., 1994. Generalized model function for quantitative texture analysis. In: Bunge, H.J., Siegesmund, S., Skrotzki, W., Weber, K. (Eds.), Textures of Geological Materials. DGM Informationsgesellschaft, Oberursel, pp. 15 – 28.
268
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
Freeden, W., Gervens, T., Schreiner, M., 1998. Constructive Approximation on the Sphere with Applications to Geomathematics. Clarendon Press, Oxford. ¨ ber Fla¨chen mit lauter geschlossenen geoda¨tischen Funk, P., 1913. U Linien. Math. Ann. 74, 278 – 300. ¨ Uber eine geometrische Anwendung der Funk, P., 1916. U Abel’schen Integralgleichung. Math. Ann. 77, 129 – 135. Gel’fand, I.M., Minlos, R.A., Shapiro, Z.Ya., 1963. Representations of the Rotation and Lorentz Groups and Their Applications. Pergamon, Oxford. Hochstadt, H., 1986. The Functions of Mathematical Physics. Wiley, New York. Ivanova, T.M., Nikolayev, D.I., 2001. New standard function for quantitative texture analysis. Phys. Stat. Solidi, B 228, 825 – 836. Mardia, K.V., Jupp, P.E., 2000. Directional Statistics. Wiley, New York. Matthies, S., 1980. Standard functions in texture analysis. Phys. Stat. Solidi, B 101, K111 – K115. Matthies, S., Vinel, G.W., Helming, K., 1987. Standard Distributions in Texture Analysis, vol. I. Akademie Verlag, Berlin. Matthies, S., Helming, K., Steinkopf, T., Kunze, K., 1988. Standard distributions for the case of fibre textures. Phys. Stat. Solidi, B 150, K1 – K5. Matthies, S., Vinel, G.W., Helming, K., 1990. Standard Distributions in Texture Analysis, vol. III. Akademie Verlag, Berlin. Mhaskar, H.N., Narcowich, F.J., Prestin, J., Ward, J.D., 2001. Polynomial frames on the sphere. Adv. Comput. Math. (in press). Nikolayev, D.I., Savyolova, T.I., 1997. Normal distribution on the rotation group SO(3). Textures Microstruct. 29, 201 – 233. Potts, D., Steidl, G., Tasche, M., 1996. Kernels of spherical harmonics and spherical frames. In: Fontanella, F., Jetter, K., Laurent, P.-J. (Eds.), Advanced Topics in Multivariate Approximation. World Scientific Publications, Singapore, pp. 287 – 301. ¨ ber die Bestimmung von Funktionen durch Radon, J., 1917. U ihre Integralwerte la¨ngs gewisser Mannigfaltigkeiten. Berichte
u¨ber die Verhandlungen der Ko¨niglich Sa¨chsischen Gesellschaft der Wissenschaften zu Leipzig. Math.-Phys. Klasse 69, 262 – 277. Savelova, T.I., [identical with Savyolova], 1984. Distribution functions of grains with respect to orientation in polycrystals and their Gaussian approximations. Ind. Lab. 50, 468 – 474 (translated from Zavod. Lab. 50, 48, 52). Schaeben, H., 1996a. Texture modeling or texture approximation with components represented by the von Mises – Fisher matrix distribution on SO(3) and the Bingham distribution on S4+. J. Appl. Crystallogr. 29, 516 – 525. Schaeben, H., 1996b. A unified view of methods to resolve the inverse problem of texture goniometry. Textures Microstruct. 25, 171 – 181. Schaeben, H., 1997. A simple standard orientation density function: the hyperspherical de la Valle¨e Poussin kernel. Phys. Stat. Solidi, B 200, 367 – 376. Schaeben, H., 1999. The de la Valle´e Poussin standard orientation density function. Proceedings of the International Conference on Neutron Texture and Stress Analysis, Dubna, Russia, June 23 – 26, 1997. Textures Microstruct. 33, 365 – 373 (Special Issue). Schaeben, H., Nikolayev, D.I., 1998. The central limit theorem in texture component fit methods. Acta Appl. Math. 53, 59 – 87. Schaeben, H., Potts, D., Prestin, J., 2001. Spherical wavelets with an application in preferred crystallographic orientation. Proceedings of the Annual Conference of the International Association for Mathematical Geology, Cancun, Mexico, Sep 6 – 12, 2001. Vilenkin, N.J., 1968. Special Functions and the Theory of Group Representations. Am. Math. Soc. Transl. of Math. Monographs 22, Providence, 613 pp. Watson, G.S., 1983. Statistics on Spheres. Wiley, New York. Zheng, Y., Doerschuk, P.C., 2000. Explicit computation of orthonormal symmetrized harmonics with applications to the identity representation of the icosahedral group. SIAM J. Math. Anal. 32 (3), 538 – 554.
Spherical harmonics in texture analysis Helmut Schaeben *, K. Gerald van den Boogaart Geoscience Mathematics and Informatics, Freiberg University of Mining and Technology, Gustav-Zeuner-Str. 12, Freiberg D-09596, Germany Accepted 31 March 2003
Abstract The objective of this contribution is to emphasize the fundamental role of spherical harmonics in constructive approximation on the sphere in general and in texture analysis in particular. The specific purpose is to present some methods of texture analysis and pole-to-orientation probability density inversion in a unifying approach, i.e. to show that the classic harmonic method, the pole density component fit method initially introduced as a distinct alternative, and the spherical wavelet method for highresolution texture analysis share a common mathematical basis provided by spherical harmonics. Since pole probability density functions and orientation probability density functions are probability density functions defined on the sphere X3oR3 or hypersphere X4oR4, respectively, they belong at least to the space of measurable and integrable functions L1(Xd), d = 3, 4, respectively. Therefore, first a basic and simplified method to derive real symmetrized spherical harmonics with the mathematical property of providing a representation of rotations or orientations, respectively, is presented. Then, standard orientation or pole probability density functions, respectively, are introduced by summation processes of harmonic series expansions of L1(Xd) functions, thus avoiding resorting to intuition and heuristics. Eventually, it is shown how a rearrangement of the harmonics leads quite canonically to spherical wavelets, which provide a method for high-resolution texture analysis. This unified point of view clarifies how these methods, e.g. standard functions, apply to texture analysis of EBSD orientation measurements. D 2003 Elsevier B.V. All rights reserved. Keywords: Real-valued spherical harmonics; Real-valued harmonics for the rotation group; Harmonic series expansions; Spectral decomposition of the X-ray transform; Summability; Spherical singular integral; Spherical radial basis functions; Spherical wavelets; Methods of texture analysis
1. Introduction The subject of texture analysis is the experimental determination (collecting and processing data) and interpretation of the statistical distribution of orienta-
* Corresponding author. Fax: +49-3731-394-067. E-mail address: [email protected] (H. Schaeben).
tions of crystals within a specimen of polycrystalline materials, which could be metals, rocks or bio-materials. It is often referred to as analysis of crystallographic lattice preferred orientation (LPO), and the objective of the interpretation is to relate an observed pattern of preferred orientation to its generating processes or vice versa. The orientation of an individual crystal is assumed to be unique and given by the rotation g 1aSO(3) which brings a coordinate system Ks fixed to the specimen into coincidence with a coordinate
0040-1951/03/$ - see front matter D 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0040-1951(03)00190-2
254
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
system Kc fixed to the crystal, g 1: KsiKc. Thus, in the context of texture analysis, rotation always means passive or ‘‘frame’’ rotation. The coordinates of a unique direction represented by haX3oR3 with respect to the crystal coordinate system Kc (referred to a crystallographic direction) and by raX3 with respect to the specimen coordinate system Ks (referred to as specimen direction) are related to each other by rotations in the set of rotations Cðh; rÞ ¼ fqaSOð3Þ j h ¼ grg which constitutes a one-dimensional circle on X4oR4 with center O parametrized by h and r when rotations are represented as unit vectors (‘‘quaternions’’) on X4. More specifically, conventional texture analysis is the determination and interpretation of the orientation probability density function f of a polycrystalline specimen by volume without individual orientation measurements; it is emphasized that size, shape, and spatial location of the contributing grains are not considered. In X-ray diffraction experiments, the orientation probability density function f cannot be directly measured, but with a texture goniometer only pole probability density functions P˜(h, r) can be sampled, which represent the probability that an experimentally fixed crystal direction h or its antiparallel h statistically coincide with a specimen direction r. Then, a pole probability density function is basically the mean of two spherical X-ray transforms 1 ˜ Pðh; rÞ ¼ ððXf Þðh; rÞ þ ðXf Þðh; rÞÞ 2 of the orientation probability density function f, where the spherical X-ray transform is defined as mean of f along the one-dimensional circle C(h, r), i.e. ðXf Þðh; rÞ ¼
1 2p
Z
f ðgÞdg
ð1Þ
fgaSOð3Þjh¼grg
In mathematical tomography, the transform assigning to a function defined on a d-dimensional manifold its mean values with respect to the family of dVdimensional submanifolds with 1 V dVV d 1 is referred to as dV-plane transform. If dV= 1, it is explicitly called X-ray transform; if dV= d 1, it is explicitly called Radon transform. For d = 2, the two cases coincide. The start of spherical mathematical tomog-
raphy dates back to the pioneering papers by Funk (1913, 1916), while mathematical tomography in an Euclidean setting commenced later with the classic paper by Radon (1917). In this sense, the term ‘‘X-ray transform’’ applies to texture analysis and is used by us (cf. Cerejeiras et al., 2002), even though it does not refer to the actual radiation, which could be c, neutron, or synchrotron as well. The invariant Haar measure dg in Eq. (1) is uniquely defined by the postulate that the measurements should essentially be independent of the choices of the involved coordinate systems Ks, Kc, which in turn requires an invariance of the form ðX ½f ðBÞ Þðh; rÞ ¼ ðX ½f ðgc1 Bgs Þ Þðgc h; gs rÞ Corrected experimental X-ray, neutron or synchrotron diffraction intensity data are thought of as being discretely sampled with a texture goniometer from continuous pole probability density functions P˜(h, r), i.e. even probability density functions defined on the cross-product X3 X3 of two unit spheres in R3. With respect to the diffraction experiment, the feasible crystal directions are the normals of the crystallographic lattice planes. Several representations of pole probability density functions exist, among others (i) truncated series expansion into spherical harmonics, (ii) finite series expansion into different standard functions, (iii) series expansion into spherical wavelets. These representations lead to practical methods to resolve the inverse problem of texture goniometry to reconstruct a reasonable orientation probability density function f defined on SO(3), or equivalently on the hypersphere X4oR4, from the given corrected intensity data, which are mean values of f along onedimensional great circles of X4 parametrized by h and r. Each method exploits that the X-ray transform of the harmonics for the group SO(3) are the harmonics for the sphere X3.
2. Spherical harmonics To compute orthonormal real-valued and symmetrized harmonics is not only of special interest in texture analysis but of general mathematical interest (cf. Zheng and Doerschuk, 2000).
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
The purpose of this section is to clarify the simple and elegant structure of the harmonic functions for the rotation group SO(3) and their straightforward relation to the spherical harmonics for the sphere X3 mediated by the spherical X-ray transform (Eq. (1)). In contrast to more conventional presentations of this subject (cf. Bunge, 1982; Matthies et al., 1987), harmonics are introduced here directly as real functions. Introduced in this way, they apply most readily to texture analysis, and unnecessary efforts with complex functions are avoided (cf. Boogaart, 2001). A fundamental theorem of representation theory assures that a realvalued representation exists (cf. Curtis and Reiner, 1962). An orientation may given by three Euler angles providing the angles of an sequence of frame rotations about some axes of the coordinate system ordered according to some convention. These rotations are given by 1 0 1 0 0 C B C B C Rx ðuÞ ¼ B 0 cosðuÞ sinðuÞ C B A @ 0 sinðuÞ cosðuÞ
255
Thus, a rotation about the z-axis can be achieved by first rotating the z-axis with Rx( p/2) = Rtx (p/2) onto the y-axis, then applying the rotation Ry(u) by the given angle u about y, and eventually rotating the yaxis. Analogously with Qy ¼ Ry
p 2
;
Qz ¼ R z
p 2
it holds Rx ðuÞ ¼ Qtz Ry ðuÞQz ¼ Qy Rz ðuÞQty Ry ðuÞ ¼ Qy Rx ðuÞQty ¼ Qtx Rz ðuÞQx Rz ðuÞ ¼ Qty Rx ðuÞQy Then any rotation gaSO(3) may be decomposed into g ¼ Rz ðu1 ÞRx ðUÞRz ðu2 Þ ¼ Rz ðu1 ÞQy Rz ðUÞQty Rz ðu2 Þ ¼ Rz ðaÞRy ðbÞRz ðcÞ
0 B B Ry ðuÞ ¼ B B @
cosðuÞ
0
sinðuÞ
0
1
0
sinðuÞ 0
C C C C A
0 cosðuÞ
sinðuÞ
cosðuÞ
1
0
1
B C B C B Rz ðuÞ ¼ B sinðuÞ cosðuÞ 0 C C @ A 0 0 1 p p ¼ Rx Ry ðuÞRtx ¼ Qx Ry ðuÞQtx 2 2 with
Qx ¼ R x
0
1
B B ¼B B0 2 @ 0
p
0 0 1
0
1
C C 1C C A 0
¼ Rz ðaÞQtx Rz ðbÞQx Rz ðcÞ Harmonics for SO(3) can generally be introduced by application of group representation theory (cf. Curtis and Reiner, 1962), i.e. the theory of all possible, essentially different representations of elements of SO(3) by matrices of arbitrary order. Harmonics provide a set of representatives of all irreducible representations, which in turn are the characteristic representations. This set is not unique; an additional specification is required to impose uniqueness. Following the complex approach of Gel’fand et al. (1963) and Vilenkin (1968), the common application in texture analysis is to transform initial complex harmonics into real harmonics (cf. Bunge, 1969, 1982), which, however, are no longer harmonic representations of rotations, i.e. compared with the complex harmonics, they do not possess the single most important property of being a representation. For the proper transform preserving the representation property, the reader is referred to Boogaart (2001).
256
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
The real and complex harmonic functions for the rotation group have a very similar structure. For real harmonics, we start with the harmonic representation of a simple rotation about the z-axis S RS ðuÞ ¼ ðrmn Þm¼S ;...;S ;n¼S ;...;S 1 0 cosðS uÞ 0 sinðS uÞ C B C B C B O q C B C B C B C ¼B 0 1 0 C B C B C B C B q O C B A @ sinðS uÞ 0 cosðS uÞ
The same object looks a little less intuitive for complex harmonics like C RS C ðuÞ ¼ ðrmnS Þm¼S ;...S ;n¼S ;...;S 0 iS u e B B B O B B B ¼B 1 B B B B O B @
0
1 C C C C C C C C C C C C A
eiS u
0
Together with a rotation Q by p/2 about another axis, harmonics can be rewritten similarly to rotation matrices as TS ðgÞ ¼ ðTSmn ðgÞÞm¼S ;...;S ;n¼S ;...;S ¼ RS ðu1 ÞQSx RS ðUÞQSx t RS ðu2 Þ
and QSy wRSy
QSx wRSx
2 p 2
As with rotation matrices, harmonics in matrix notation possess a distinguished property also known as the representation property (cf. Gel’fand et al., 1963; Vilenkin, 1968) TS ðggVÞ ¼ TS ðgÞTS ðgVÞ In texture analysis, the orientation probability density function f is expanded into a Fourier series of harmonic functions. With the trace of a matrix trA= P trðaij Þij ¼ i aii , this expansion can be written in condensed matrix notation as
f ðgÞ ¼
l X S S X X
CSmn TSmn ðgÞ ¼
S ¼0 m¼S n¼S
l X
trðCtS TS ðgÞÞ
S ¼0
with matrices CS of CSmn coefficients CS ðCSmn Þm¼S ;...;S ;n¼S ;...;S The representation property of the harmonics allows to calculate the harmonic expansion of the orientation probability density function with respect to a different sample or crystal coordinate system by simple matrix calculations as follows
f˜ ðgÞwf ðggVÞ ¼
l X
trðCtS TS ðgÞTS ðgVÞÞ
S ¼0
¼ RS ðaÞQSy RS ðbÞQSy t RS ðcÞ ¼
with
p
l X
trðTS ðgVÞCtS TS ðgÞÞ
S ¼0
RSz ðuÞwTS
ðRz ðuÞÞ ¼ R ðuÞ S
¼
trððCS TtS ðgVÞÞt TSmn ðgÞÞ |fflfflfflfflffl{zfflfflfflfflffl} S ¼0 CSmn of f ðgVÞ
RSx ðuÞwTS ðRx ðuÞÞ ¼ QSy RS ðuÞQSy t RSy ðuÞwTS ðRy ðuÞÞ ¼ QSx t RS ðuÞQSy
l X
¼
l X S ¼0
˜ t TS ðgÞÞ trðC S
ð2Þ
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
Let us now focus on X-ray transforms and pole probability density functions, respectively. A crystallographic direction h may be presented by two Euler angles u, U corresponding to the rotations Rz(u) and Rx(U) which rotate the z-axis of Kc onto h.
Analogously, a specimen direction r can be defined by two Euler angles uV, UV, providing the rotations Rz(uV) and Rx(U)Vwhich rotate the z-axis of Ks onto r. r ¼ Rz ðuVÞRx ðUVÞez Calculating (Xf)(h, r) is most simple for h = r = ez, where ez denotes the unique vector representing the direction of the coinciding z-axis of the coordinate systems. Then 1 2p
Z
˜ S e0 C˜ S00 ¼ e0 C ¼ tret0 RSz ðuVÞRSx ðUVÞCS ðRSz ðuÞRSx ðUÞÞt e0 ¼ trCtS ðRSx ðUVÞÞt ðRSz ðuVÞÞt e0 et0 Rz ðuÞRz ðUÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffi t 4p 4p 2S þ1 YS
h ¼ Rz ðuÞRx ðUÞez
f ðRz ð/ÞÞd/ Z 2p l X tr CS RS ð/Þd/ ¼ S ¼0
l X S 1 X þ CSmm cosðm/Þd/ ¼ 2p S ¼0 m¼S Z 2p Z 2p CSmm sinðm/Þd/ þCSmm cosðm/Þd/ 0 0 |fflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} l 1 X ¼ C 00 2p S ¼0 S
0
l X
4p trCSt YtS ðhÞYS ðrÞ 2S þ 1 S ¼0
l X S S X X CSmn YS ;m ðhÞYS ;n ðrÞ 2S þ 1 S ¼0 m¼S n¼S
ð3Þ
is retrieved; for a complementary interpretation and application the reader is referred to (Boogaart and Schaeben, 1999). Harmonics, being real-valued right from their introduction, may then be symmetrized in much the same way as their complex-valued harmonic relatives.
0 when m p 0
Now using u, U, uV, UV to rotate the coordinate systems of the specimen and the crystal, respectively, such that h = ezaKc and r = ezaKs, the orientation probability density function in the rotated space is given by f˜ ðgÞ ¼ f ððRz ðuVÞRx ðUVÞÞt gRz ðuÞRx ðUÞÞ According to formula (2), the C-coefficients of ˜f are given by ˜ S ¼ RS ðuVÞRS ðUVÞCS ðRS ðuÞRS ðUÞÞt CS RS C z x z x ðuÞQSy RS ðUÞQSy
ðhÞ
which could be used as definition of spherical harmonics, the well-known spectral decomposition of the X-ray transform (Bunge, 1969, 1982; Matthies et al., 1987)
¼ 4p
Z 2p l X S 1 X ¼ CSmm sinðm/Þ 2p S ¼0 m¼S 0
2S þ1 YS
rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2S þ 1 S YS ðhÞ ¼ ðRx ðUVÞÞt ðRSz ðuVÞÞt e0 4p
ðXf Þðh; rÞ ¼
0
ðrÞ
Thus, setting
2p 0
257
3. Standard functions in texture analysis According to the intuitive introduction of standard orientation probability density functions into texture analysis by the Rossendorf school (cf. Matthies, 1980; Matthies et al., 1987, 1988, 1990; Eschner, 1993, 1994), it should satisfy the requirements (i) to be specified by a small number of parameters, (ii) mathematically and numerically tractable, (iii) generally appealing from the point of physical interpretation, and (iv) possess analytical expressions for corresponding pole probability density functions and all normalization constants, thus avoiding truncation errors. This informal characterization includes that their harmonic series expansion
258
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
is known. Later, all Rossendorf standard orientation probability density functions widely applied in pole probability density (or ambitiously: texture) component fit methods have been shown to actually reduce to special cases of the von Mises– Fisher matrix distribution on SO(3) or, equivalently, to the Bingham distribution of axes on X+4oR4, particularly to the symmetrically unimodal and the symmetrically circular ‘‘fibre’’ case (Schaeben, 1996a). Thus, the von Mises– Fisher matrix and the Bingham distribution, respectively, provide a general mathematical model orientation probability density function. Similarly, other standard functions like the spherical Brownian (Savelova, 1984; Nikolayev and Savyolova, 1997; Schaeben and Nikolayev, 1998) and the spherical de la Valle´e Poussin (Schaeben, 1997, 1999) have been introduced, which mark with respect to harmonic series expansion and analytical representation of the limiting function the ‘‘extreme ends’’ on the scale of standard functions: There is no analytical representation for the Brownian, while the harmonic series expansion of the de la Valle´e Poussin is finite. Recently, another standard probability density function proportional to exp(j tan2 (x/2)) has been introduced (Ivanova and Nikolayev, 2001). Complementary to their initial introduction, standard functions are derived here by summability or summation processes, respectively, of harmonic series expansion revealing their intricate relation with them and avoiding resorting to intuition and heuristics. 3.1. Notation Basically following Freeden et al. (1998) and Hochstadt (1986), some basic notation is recalled and fundamental results are summarized. The unit sphere in d dimensions is denoted Xd and defined as Xd ¼ fraRd j ArA ¼ 1g If e1, . . ., ed are orthonormal vectors in Rd, and raXd, then pffiffiffiffiffiffiffiffiffiffiffiffi r ¼ ted þ 1 t 2 rd1
where t = rted and rd 1 is a unit vector in the space spanned by e1, . . ., ed 1, and dxd ¼ ð1 t 2 Þðd3Þ=2 dtdxd1 (cf. Watson, 1983). Thus NXd N ¼
Z Xd
dxd ¼
2pd=2 Cðd=2Þ
e.g. NX2N = 2p, NX3N = 4p, NX4N = 2p2. The class L p (Xd) consists of scalar functions f: d X iR1, which are measurable and for which Nf NLp ðXd Þ ¼
Z Xd
Af ðxÞAp dxp ðxÞ
1=p
Thus, for d = 3, p = 1, it holds pN1N ffiffiffiffiffiffi L1(X3) = 4p, for d = 3, p = 2 it holds N1NL2 ðX3 Þ ¼ 4p. Since our main interest is in results for X3oR3, which apply to experimentally accessible pole probability density functions, omitting the superscript and writing in short X always refers to X3. Most results analogously generalize to d >3. Real spherical harmonics of different orders are orthogonal in the sense of the L2 inner product Z ðYn ; Ym ÞL2 ðX3 Þ ¼ Yn ðrÞYm ðrÞdx3 ðrÞ ¼ 0 X3
There are 2n + 1 linearly independent harmonics of order n for X3; therefore, let { Yn, 1, . . ., Yn, 2n + 1} denote an orthonormal sequence of harmonics of order n with respect to the inner product of L2 (X3), and let Yn,j, j = 1, . . ., 2n + 1 denote one of its members. Then it holds n X j¼n
Yn;j ðr1 ÞYn;j ðr2 Þ ¼
2n þ 1 Pn ðr1 r2 Þ 4p
where Pn denotes the Legendre polynomial (cf. Freeden et al., 1998; p. 32 (2.3.12), p. 38). It can be written as 2 1 d ðt 2 1Þn ; Pn ðtÞ ¼ n 2 n! dt
ta½1; 1
(cf. Freeden et al., 1998; p. 41 (3.2.14)).
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
As defined above, a pole probability density function is actually a probability density function P˜(h, r) defined on X3 X3. Mathematically, haX3 and raX3 are any two-unit vectors in R3. With respect to the diffraction experiment, however, h refers to the unit normal direction of a crystallographic lattice plane and therefore cannot vary arbitrarily as can raX3. To emphasize the role of h as a parameter of the experiment, pole probability density functions used to be denoted P˜h(r) or P˜h(r) to indicate symmetry of the crystallographic form h. In the following, harmonic analysis is applied to functions F(h, r) defined on X3 X3 read as a function of r parametrized by h. Then, any function F(h, r)aX (X3 X3), where X denotes either Lp, 1 V p < l, or C, can formally be associated with its Fourier series, i.e. Fðh; rÞf
l X n X
ðFðh;BÞ; Yn;j ðBÞÞYn;j ðrÞ
integrals, i.e. to convolutions of the function with a zonal kernel function approximating identity. The nth partial sum of the harmonic series expansion is defined as Sn ðFðh; BÞ; rÞ ¼
with constants (independent of r)
^ Fn;j ðhÞYn;j ðrÞ
S ¼0 j¼S
¼
n X S X
ðFðh; BÞ; Yn;j ðBÞÞYn;j ðrÞ
S ¼0 j¼S
¼
Fðh; BÞ;
n X S X
! Yn;j ðBÞYn;j ðrÞ
S ¼0 j¼S n X 2S þ 1 PS ðr BÞ Fðh; BÞ; 4p S ¼0
¼
Gn ðtÞ ¼
^ ðhÞ ¼ ðFðh; BÞ; Yn;j ðBÞÞ Fn;j
! ð4Þ
n X ð2S þ 1ÞPS ðtÞ S ¼0
referred to as its spherical Fourier-transform or Fourier-coefficients. Thus, Fðh; rÞf
n X S X
Setting
n¼0 j¼n
l X n X
259
¼ ðn þ 1Þ
Pnþ1 ðtÞ Pn ðtÞ t1
ð5Þ
which may be addressed as the spherical Dirichlet kernel, it is
^ Fn;j ðhÞYn;j ðrÞ
n¼0 j¼n 2
3
It can be shown that any function F(h, B)aL (X ) may be approximated in the L2 (X3) sense, i.e. in the norm of L2 (X3), by finite truncations of its Fourier orthogonal series expansion in terms of any L2 (X3) orthonormal system of spherical harmonics Yn,j. It should be noted that the above statement does not hold true for every continuous F(h, B)aC (X3) or F(h, B)aLp (X3) with pa[1, 4/3]v[4, l) (Freeden et al., 1998, p. 58). Thus, finite truncations of Fourier series expansions by means of L2 (X3)-orthonormal spherical harmonics do not generally form an ‘‘approximation process’’ on the sphere. For convergence in C(X3) norm it is required that F is Lipschitzcontinuous (Freeden et al., 1998, p. 59). Instead of imposing conditions on the function to be approximated, it may be more efficient to augment the Fourier series with ‘‘convergence factor’’ to ensure convergence. As it will be shown, the augmented series expansions correspond to spherical singular
1 Sn ðFðh; BÞ; rÞ ¼ 4p
Z X3
Fðh; rVÞGn ðr rVÞdx3 ðrVÞ
Of mathematical interest is the existence and uniqueness of the limit lim Sn ðFðh; BÞ; rÞ
n!l
As already mentioned above, the series of partial sums does not generally converge for FaC(X3) or FaL1 (X3) in the sense of the corresponding norm (cf. Freeden et al., 1998). Therefore, the linear space of all absolutely convergent sequences {hq (n, j)} of real numbers hq (n, j), (n, j)aJ ={(n, j) | n = 0, 1, . . ., j = 1, . . ., 2n + 1} for q in some parameter set IA is introduced. The sequence {hq (n, j)} furnished with appropriate properties is referred to as convergence factor eventually constituting a convergent summation process of the ‘‘augmented’’ harmonic series expansion as follows.
260
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
the set of functions g: [ 1,1]iR1, and endowed with norms NgAX ½1;1 ¼ Nqðe BÞNXðXd Þ
Instead of Sn( F(h, B); r), the focus is on Uq ðFðh; BÞ; rÞ ¼
l X n X
^ hq ðn; jÞFn;j ðhÞYn;j ðrÞ
n¼0 j¼n
¼
l X n X
¼
hq ðn; jÞðFðhÞ; BÞ; Yn;j ðBÞYn;j ðrÞ
¼
Fðh;BÞ;
n X
! hq ðn; jÞYn;j ðBÞYn;j ðrÞ
ð6Þ
n¼0 j¼n
which leads to different summation processes for different choices of hq (n, j), qaIA. If hq (n, j) = hq (n) for all j = 1, . . ., 2n + 1, then Uq ( F(h, B); r) simplifies to the special case Uq ðFðh; BÞ; rÞ ¼
l X n X
^ hq ðnÞFn;j ðhÞYn;j ðrÞ
¼
Fðh; BÞ
n¼0
!
Uq ðFðh; BÞ; rÞ ¼
1 4p
hq ðnÞPn ðr BÞ
4p Z
ð7Þ
Fðh; rVÞ
X3 l X
p
2p
1=p
AgðtÞA dt
they form subspaces of X (Xd) where X (Xd) denotes either Lp (Xd) or C (Xd), and X [ 1, 1] analogously. For zonal functions, their harmonic series expansions reduce to an expansion in terms of Legendre polynomials Pn,d for Rd with corresponding Fourier– Legendre coefficients. The spherical convolution of a zonal function gaL1 [ 1, 1] and a function faX (Xd) is defined by Z ðg*f ÞðyÞ ¼ gðy xÞf ðxÞdxd ðxÞ Xd
n¼0 j¼n l X 2n þ 1
1 1
n¼0 j¼n l X
Z
!
ð2n þ 1Þhq ðnÞPn ðr rVÞ
n¼0
dxðrVÞ Z 1 ¼ Fðh; rVÞKq ðr rVÞdx3 ðrVÞ 4p X3 ð8Þ If more specifically IA = IN, then hq (n, j) = hm (n), and if additionally 8 < 1 if nVm hm ðnÞ ¼ : 0 otherwise then the summation process Um( F(h, B); r) degenerates to the partial sums Sm ( F(h, B); r). 3.2. Summation processes and singular integrals on the sphere Any scalar function gy: XdiR1, xigy(x) = g(y x), x, yaXd is called zonal (synonymously central, rotationally invariant). Zonal functions are isomorphic to
A convolution has the nicer property of the corresponding property of g and f, respectively. In particular, the Fourier coefficients of the convolution are the product of the Fourier( – Legendre) coefficients of the functions being convoluted. A set of functions {vq | qaIA}oL1 [ 1, 1] with q in some parameter set IA is called a kernel if ðvq Þ^ ð0Þ ¼ 2p
Z ½1;1
vq ðtÞdt ¼ 1
holds. The convolution of a function faX (Xd) with a kernel vqaL1 [ 1, 1] is called a spherical singular integral vq ðf ; xÞ ¼ ðvq *f ÞðxÞ ¼
Z Xd
vðx yÞf ðyÞdxd ðyÞ
If the kernel function vaX [ 1, 1] is an approximate identity, i.e. if it is uniformly bounded for all qaIA, and if it possess the spatial localization property Z
d
lim
q!q0
1
vq ðtÞdt ¼ 0;
dað1; 1Þ
then the singular integral converges in X (Xd)-norm towards f for q ! q0. A nonnegative kernel is an approximate identity if lim ðvq Þ^ ð0Þ ¼ 1
q!q0
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
With respect to texture analysis, it is well known that if the orientation probability density function f is zonal with respect to g0aSO(3), then its X-ray transform is zonal, too, and it holds ðXf Þðh; rÞ ¼ ðXf Þðh g0 rÞ ¼ ðXf Þðg01 h rÞ Correspondingly, the familiar form of the harmonic series expansion of an orientation probability density function f ðgÞf
l X S S X X S ¼0 n¼S m¼S
for any zonal function f defined on X4, where the last equation is an application of the identity CS(1) (cos w) = US (cos w) with the Chebychev polynomials of the second kind US (cos w). Further, following an argument introduced by Arnold (1941), the corresponding series expansion for an even function f is 1 ½ f ðx; jÞ þ f ðp x; jÞ
2 l 1 X ¼ 2 ðS þ 1Þ2 f ^ ðS ÞPS ;4 ðcosxÞ 2p S ¼0ð2Þ
f ðx; jÞ ¼
4p C mn T mn ðgÞ 2S þ 1 S S
l 1 X ð2S þ 1Þ2 f ^ ð2S ÞP2S ;4 ðcosxÞ 2p2 S ¼0 l 1 X sin½ð2S þ 1Þx
¼ 2 ð2S þ 1Þf ^ ð2S Þ 2p S ¼0 sinx
¼
and its X-ray transform ðXf Þðh; rÞf
l X S S X X
CSmn YS ;m ðhÞYS ;n ðrÞ
ð12Þ
S ¼0 n¼S m¼S
is turned into f ðg; g0 ; jÞf
261
with the even-order Chebychev polynomials of second kind U2n (x/2)=((sin (2n + 1)x/2)/(sin x/2)). Formally, a spherical singular integral may be associated with its harmonic series
l X
2S þ 1 ^ f ð2S Þ 2p2 S ¼0
sin½ð2S þ 1Þxðgg01 Þ=2
sinðxðgg01 Þ=2Þ
ð9Þ vq ðF; rÞf
l X n X
ðvq ðF; BÞ; Yn;j ðBÞÞYn;j ðrÞ
n¼0 j¼n
and l X 2S þ 1 ðXf Þ^ ðX ½f ðB; g0 ; jÞ Þðh; rÞf 4p S ¼0
ðS ÞPS ðh g0 rÞ ¼ ðXf Þðh g0 r; jÞ
¼
l X n X
ððvq *FÞðBÞ; Yn;j ðBÞÞYn;j ðrÞ
n¼0 j¼n
ð10Þ
¼
l X n X
^ v^q ðnÞFn;j ðhÞYn;j ðrÞ
ð13Þ
n¼0 j¼n
(cf. Matthies et al., 1987; Schaeben, 1996b). Eq. (9) follows from f ðx; jÞ ¼
l 1 X ðS þ 1Þ2 f ^ ðS ÞPS ;4 ðcosxÞ 2 2p S ¼0
¼
l 1 X ð1Þ ðS þ 1Þf ^ ðS ÞCS ðcosxÞ 2p2 S ¼0
¼
l 1 X ðS þ 1Þf ^ ðS ÞUS ðcosxÞ 2p2 S ¼0
¼
with constants (independent of r)
1 2p2
l X
ðS þ 1Þf ^ ðS Þ
S ¼0
sin½ðS þ 1Þx
sinx
v^q ðnÞ ¼ 2p
Z
1
1
vq ðtÞPn ðtÞdt ¼ ðG; Pn Þ
and ^ Fn;j ðhÞ ¼ ðFðh; BÞ; Yn;j ðBÞÞ
ð11Þ
referred to as Fourier( –Legendre) coefficients. Comparing Eqs. (7) and (13), at least the formal analogy is obvious that hq (n, k) is a convergence factor and provides a summation process if it is equal
262
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
to the Fourier – Legendre coefficients of a (zonal) approximate identity vq (t). For detailed results and proofs the reader is referred to Dunkl (1966), Berens et al. (1968), Butzer and Nessel (1971), Butzer et al. (1979), and Freeden et al. (1998). 3.3. Examples For reasons of notational simplicity and practical applications, we confine ourselves again to the case X3, i.e. to X-ray transforms. The most simple choice for a spherical kernel is the Dirichlet kernel (Eq. (5)). Choosing more sophisticatedly as spherical kernel the (i) spherical Brownian, (ii) spherical de la Valle´ e Poussin, (iii) spherical Cauchy, (iv) von Mises – Fisher and (v) X-ray transformed von Mises-Fisher probability density function l 1 X ð2S þ 1Þexp½S ðS þ 1Þq PS ðtÞ (i) vq ðtÞ ¼ 4p S ¼0
l X
1 I S ðqÞ I S þ1 ðqÞ PS ðtÞ ð2S þ 1Þ 4p S ¼0 2B 32 ; q þ 12 1 1þt q ¼ ¼ vn ðtÞ Bð1; q þ 1Þ 2
(iii) vq ðtÞ ¼
l 1 X ð2S þ 1Þq2S PS ðtÞ 4p S ¼0 2
1 1q ¼ pq ðtÞ 4p ð1 þ q2 2qtÞ3=2 l I S þ1=2 ðqÞ 1 X PS ðtÞ ð2S þ 1Þ (iv) vq ðtÞ ¼ I 1=2 ðqÞ 4p S ¼0 ¼
¼
3.3.1. Singular integral of Gauss –Weierstrass, spherical Brownian standard probability density function Choosing hq ðn; jÞ ¼ hq ðnÞ ¼ expðnðn þ 1ÞqÞ with qaIA=(0, l), q0 = 0, it follows that l 1 X ð2n þ 1Þhq ðnÞPn ðtÞ 4p n¼0 ¼
l 1 X ð2n þ 1Þexpðnðn þ 1ÞqÞPn ðtÞ 4p n¼0
¼ bðt; qÞ
¼ bðt; qÞ (ii) vq ðtÞ ¼
leads to the singular integrals of (i) Gauss– Weierstrass, (ii) de la Valle´e Poussin, (iii) Abel – Poisson, (iv) von Mises– Fisher and (v) X-ray transformed von Mises– Fisher. All of them provide summation processes and thus accomplish convergence in the X (Xp) norm. In case of C (Xd), this is uniform convergence which in turn implies pointwise convergence.
q expðqtÞ ¼ mðt; qÞ 4psinhq
which is the spherical Brownian distribution for X3 (cf. Mardia and Jupp, 2000). Moreover, it is the X-ray transform of the spherical Brownian for X4 (cf. Schaeben, 1996b). Then Z Wq ðFðh; BÞ; rÞ ¼ Fðh; rVÞbðr rV; qÞdxðrVÞ ð15Þ X3
is the spherical singular integral of Gauss– Weierstrass (cf. Freeden et al., 1998; Schaeben, 1996b). It converges for q ! 0 to F(h, r) in the corresponding norm. A texture component fit method using spherically Brownian distributed components evolves as follows. In case of a finite data set F(h, rk) according to some spherical sampling scheme rkaSh adjusted to q wq ðFðh; BÞ; rÞ ¼
(v) vq ðtÞ ¼ ¼
l 1 X I S ðqÞ I S þ1 ðqÞ PS t ð2S þ 1Þ 4p S ¼0 I 0 ðqÞ I 1 ðqÞ
q q 1 I0 ð1 þ tÞ exp ðt 1Þ 4psinhq 2 2
¼ ðXM Þðt; qÞ
ð14Þ
NS hN X
ak ðhÞbðr rk ; qÞ
k¼1
provides a practical approximation of F(h, r) in the sense that NS hN X k¼1
ak ðhÞbðr rk ; qÞcFðh; rk Þ
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
The iterative fit of ak (h) may start with a(0) k (h) = F(h, rk). In practice, the series expansion may be truncated by the user to a rough approximation by K major texture components with q = q(k), k = 1, . . ., K uðFðh; BÞ; rÞ ¼
K X
ak ðhÞbðr rk ; qðkÞÞ
where a small number K bNShN denotes the total number of spherically Brownian distributed components to be considered. Usually, these major components are placed in some rk where F(h, r) displays strong local maximum. Strictly speaking, this usercontrolled method turns texture approximation by finite series expansion into texture modeling by a very small number of texture components, which may be very successful and helpful in practice when the pattern of lattice preferred orientation is known or sufficiently simple. 3.3.2. Singular integral of de la Valle´e Poussin, spherical de la Valle´e Poussin standard probability density function Choosing 8 ðjÞ > < I n ðjÞI 3 nþ1 1 if nVj 2B ; j þ hq ðn; jÞ ¼ hj ðnÞ ¼ 2 2 > : 0 if n > j where
0 1 2S S X B C I S ðjÞ ¼ ð1Þk @ A k¼0 2k 1 1 B k þ ;j þ S k þ ; 2 2
Z X3
Fðh; rVÞvn ððr rV; qÞdxðrVÞÞ ð17Þ
is the spherical singular integral of de la Valle´e Poussin (cf. Schaeben, 1996b). It converges for n ! l to F(h, r) in the corresponding norm. A texture component fit method using spherically de la Valle´e Poussin distributed components evolves analogously to the case of Gauss– Weierstrass. Since the infinite series representing the kernel reduces to a finite series, truncation is not required in practical numerical applications involving CS coefficients, and errors due to truncation are avoided. This property essentially distinguishes the de la Valle´e Poussin standard orientation probability density function from other known standard functions in an advantageous way, as it provides a computationally efficient model. 3.3.3. Singular integral of Abel – Poisson, spherical Cauchy- or Lorentz-type standard probability density function Choosing hq ðn; jÞ ¼ hh ðnÞ ¼ hn with ha[0,1), h0 = 1, it follows that l 1 X ð2n þ 1Þhh ðnÞPn ðtÞ 4p n¼0
S ¼ 0; 1; . . . with jaIA = N, j0 = l, it follows that n 1 X ð2S þ 1Þhj ðS ÞPS ðtÞ 4p S ¼0 n 1 X I S ðjÞ I S þ1 ðjÞ PS ðtÞ ð2S þ 1Þ ¼ 4p S ¼0 2B 32 ; j þ 12 1 1þt j ¼ ¼ vn ðtÞ Bð1; j þ 1Þ 2
which is the spherical de la Valle´e Poussin kernel for X3 and the X-ray transform of the spherical de la Valle´e Poussin kernel for X4 (cf. Schaeben, 1997, 1999). Then Vn ðFðh; BÞ; rÞ ¼
k¼1
263
ð16Þ
¼
l 1 X ð2n þ 1Þhn Pn ðtÞ 4p n¼0
¼
1 1 h2 ¼ ph ðtÞ 4p ð1 þ h2 2htÞ3=2
ð18Þ
which is the spherical Cauchy distribution for X3 (cf. Mardia and Jupp, 2000) or Lorentz-type distribution (Matthies et al., 1987), respectively; it is also the X-
264
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
ray transform of the spherical Cauchy distribution for X4. Then Z Fðh; rVÞph ðr rVÞdxðrVÞ ð19Þ Uh ðFðh; BÞ; rÞ ¼ X3
is the spherical singular integral of Abel –Poisson (cf. Freeden et al., 1998; Schaeben, 1996b). A texture component fit method using spherically Cauchy distributed components evolves analogously to the case of Gauss – Weierstrass or de la Valle´e Poussin. Owing to the poor localization, i.e. the fat tails, of the Cauchy distribution, it is rarely used in practical texture component fit procedures. 3.3.4. Standard von Mises – Fisher probability density function Choosing Inþ1=2 ðjÞ I1=2 ðjÞ
hq ðn; jÞ ¼ hj ðnÞ ¼
with jaIA = R+, j0 = l, it follows that
¼
Fðh; rVÞMðr rV; jÞdxðrVÞ
l 1 X ð2n þ 1Þhq ðnÞPn ðtÞ 4p n¼0 l 1 X Il ðjÞ Ilþ1 ðjÞ Pn ðtÞ ð2n þ 1Þ 4p n¼0 I0 ðjÞ I1 ðjÞ j j ¼ CM ðjÞI0 ð1 þ tÞ exp ðt 1Þ 2 2
¼
ð20Þ
ð22Þ
Z
Fðh; rVÞXM ððr rV; jÞdxðrVÞÞ
X3
ð23Þ
which is the von Mises– Fisher distribution for X (cf. Mardia and Jupp, 2000). It is not the X-ray transform of the rotationally symmetric von Mises – Fisher matrix probability density function on SO(3) or equivalently the bimodal Bingham, i.e. Watson probability density function for X4. Then Uj ðFðh; BÞ; rÞ ¼
with jaIA = R+, j0 = l, it follows that
Uj ðFðh; BÞ; rÞ ¼
3
Z
In ðjÞ Inþ1 ðjÞ I0 ðjÞ I1 ðjÞ
hq ðn; jÞ ¼ hj ðnÞ ¼
which is the X-ray transform of the von Mises –Fisher distribution for X4. Then
l Inþ1=2 ðjÞ 1 X Pn ðtÞ ð2n þ 1Þ I1=2 ðjÞ 4p n¼0
j expðjtÞ ¼ mðt; jÞ 4psinhj
3.3.5. X-ray-transformed standard von Mises –Fisher probability density function Choosing
¼ ðXM Þðt; jÞ
l 1 X ð2n þ 1Þhj ðnÞPn ðtÞ 4p n¼0
¼
Fisher pole probability density function seems to be unknown, it is not applied in texture component fit methods.
ð21Þ
is a spherical singular integral that converges for j ! l to F(h, r) in the norm. A texture component fit method using projected von Mises– Fisher distributed pole probability density components evolves analogously to the previous cases and corresponds to von Mises– Fisher distributed orientation probability density components. 3.4. Remark
X3
which is the spherical singular with respect to the von Mises –Fisher kernel (Schaeben, 1996b). Since the form of the orientation probability density function which is projected on the von Mises –
This first class of radial spherical basis functions (Gauss – Weiertrass, de la Valle´e –Poussin, von Mises– Fisher, Abel – Poisson) is characterized by the fact that all Legendre coefficients of the corresponding kernel are nonvanishing (Freeden et al., 1998). Note
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
that this is actually true for the de la Valle´e Poussin kernel for n ! l. The second class consists of locally supported kernel functions. Depending on the diameter of their support, some of their Legendre coefficients vanish (Freeden et al., 1998). Application of the Riemann– Lebesgue kernel, which is the indicator with respect to spherical caps was briefly considered in Schaeben (1996b); it is not considered here as its presentation is beyond the scope of this contribution. Unimodal radial basis functions can be transformed and adapted to represent multimodal fibre textures; formally, the required adaptation is accomplished by interpreting the angle x = arc cos t as orientation distance xf with respect to the fibre C(h0, r0) defined by h0, r0aX3.
As shown above, any function F(h, B) thought of as being parametrized by haX defined on X may be associated with its harmonic series expansion Fðh; rÞf
l X S X
FS^;j ðhÞYS ;j ðrÞ
S ¼0 j¼S
¼
l X S X
ðFðh; BÞ; YS ;j ðBÞÞYS ;j ðrÞ
ð24Þ
S ¼0 j¼S
The space of all spherical harmonics of degree S with S = 0, 1, 2, . . ., i.e. the restriction to the unit sphere of all homogeneous and harmonic polynomials in three variables, is denoted by Y S ; then dim Y S = 2S + 1. That any function in L2(X3) may be represented by its harmonic series expansion may be rewritten as
4. Spherical wavelets
L2
l
Another representation of spherical distributions in texture analysis is provided in terms of spherical wavelets (Schaeben et al., 2001). The main idea of wavelet analysis is to obtain a multiscale representation of the data or functions, which allows localization in space and frequency. Initially, a spherical probability density function is sampled on a coarse, ideally almost equidistributed, grid on the sphere. These coarse grid measurements are approximated by a spherical polynomial of low degree. This polynomial is clearly a sufficiently good approximation in regions of the sphere where the underlying function does not vary too much, i.e. where the underlying function consists of low frequencies only. Where the data show large variations, the initial approximation requires a local improvement, which is accomplished by a local refinement of the grid of measurements. The crucial point then is to construct a high-degree polynomial from the global coarse grid and the locally refined grid. This can be seen as adding adaptively and locally a wavelet part to the global approximation of low degree. Thus, wavelets provide a digital device to zoom into areas of special interest. Moreover, wavelets will provide the means to control the scanning process with a texture goniometer adaptively to a required local refinement of the spatial resolution.
265
L2 ðXÞ ¼
Y
S
S ¼0
Further, Pn (X) = PSn = 0 Y S with dim Pn (X)= (n + 1)2. Let Nk be a sequence of strictly monotone increasing positive integers. Then, scaling spaces Vk are introduced as polynomial spaces Vk ðX3 Þ ¼ PNk ðX3 Þ with dim Vk (X3)=(Nk + 1)2. Hence, there is a chain V0 ðX3 ÞoV1 ðX3 ÞoV2 ðX3 Þo . . . and it is sensible to define the corresponding orthogonal complements Nkþ1
Wk ðX3 Þ ¼ Vkþ1 ðX3 ÞOVk ðX3 Þ ¼
S ¼Nk þ1
YS
with dim Wk=(Nk + 1 Nk)(Nk + 1 + Nk + 2). Finally, l
L ðX3 Þ ¼ V0 ðX3 ÞP 2
W ðX Þ k
L2
3
k¼0
The classical approach of wavelet theory uses Nk = 2k. However, in case of spherical probability density functions, this would imply very large dimensions of the wavelet spaces. Therefore, the choice maybe more like Nk = 20 k.
266
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
with Ik={k=(k1,k2)TaN02: k1 < 2 Nk, k2 < 2 Nk} and
Using again 2X S þ1 j¼0
2S þ 1 PS ðr1 r2 Þ; YS ;j ðr1 ÞYS ;j ðr2 Þ ¼ 4p
rk;k
r1 ; r2 aX3
T
by using the Clenshaw –Curtis weights wkV(k). Analogously, with
with Legendre polynomials PS of degree S with PS (1) = 1 results in Fðh; rÞ ¼
pk1 pk2 pk1 pk2 pk2 cos cos ; cos sin ; sin 2Nk Nk 2Nk Nk Nk
l X
Fðh; BÞ;
S ¼0
2S þ 1 PS ðr BÞ 4p
SVk ðFðh; BÞ; rÞ ¼
Fðh; BÞ;
S ¼0
2S þ 1 PS ðr BÞ 4p
Nkþ1 X
ð2S þ 1ÞPS ðtÞ
S ¼Nk þ1
ð25Þ
Considering harmonics only up to order Nk, the approximation SVk ( F(h, B); r) with respect to Vk is given by Nk X
HV k ðtÞ ¼
and SWk ( F(h, B); r) c SWI k ( F(h, B); r) by numerical integration, it holds I SVI kþ1 ðFðh; BÞ; rÞ ¼ SVI k ðFðh; BÞ; rÞ þ SW ðFðh; BÞ; rÞ k
ð30Þ
Substituting SWI k ( F(h, B); r) by ð26Þ
I Sˆ W ðFðh; BÞ; rÞ ¼ k
X
wV k þ1 ðkÞLj Fðh; rV k þ1;k Þ
kaIkþ1
HV k ðr rkVþ1;k Þ
and the corresponding next wavelet part by
SWk ðFðh; BÞ; rÞ ¼
N k þ1 X
S ¼Nk þ1
2S þ 1 PS ðr BÞ Fðh; BÞ; 4p
ð27Þ Introducing the kernel GkVðtÞ ¼
ð31Þ
with Lk Fðh; rkþ1;k Þ 8 < Fðh; rkþ1;k Þ ¼ : I SV ðFðh; BÞ; rkþ1;k Þ
if FðhÞ is large ð32Þ otherwise
Nk X ð2S þ 1ÞPS ðtÞ
eventually yields
S ¼0
I Sˆ VI kþ1 ðFðh; BÞ; rÞ ¼ SVI k ðFðh; BÞ; rÞ þ Sˆ W ðFðh; BÞ; rÞ k
it holds 1 SVk ðFðh; BÞ; rÞ ¼ 4p
Z
ð33Þ Fðh; rVÞGkVðr rVÞdxðrVÞ
X3
ð28Þ The right-hand side can be numerically evaluated by a sampling theorem (Potts et al., 1996) such that SVk ( F(h, B); r) c SVIk ( F(h, B); r), i.e. SVI k ðFðh; BÞ; rÞ ¼
X kaIk
wkVðkÞFðh; rVk;k ÞGkVðr rVk ;k Þ
Summarizing, the proposed procedure is as follows. The function F(h, r) is known on a coarse r-grid and SVI ( F(h, B); r) is computed. By inspection or algorithms (Mhaskar et al., in press), the regions where F(h, r) and SVI ( F(h, B); r) have large absolute values and large variations, respectively, are detected. In the next step, the approximation SVI F is improved by adding the next wavelet part SWI F. Since F(h, r) should not be totally sampled on a refined r-grid and since it is known that SWI F is almost zero in regions
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
where F(h, r) does not oscillate too much, SVI F + SWI F is replaced by SVI F + SˆWI F. The definition of Lk shows that the function F(h, r) must be resampled for some additional points in the regions where large variations of F(h, r) are observed or expected. It can be shown that the spherical X-ray transform of the wavelet representation of f is the wavelet representation of the X-ray transform Af(h, B), and that the wavelet representation can be applied to solve the inverse problem to determine numerically the crystallographic orientation probability density function f from corrected experimental X-ray, neutron or synchrotron intensities measured in diffraction experiments with a texture goniometer.
267
analysis. Wavelet texture analysis combines the best elements of both methods, in particular the frequency localization of the harmonic and the spatial localization of the component fit method. It generalizes both with mathematical rigor and provides a wellfounded means for high-resolution texture analysis.
Acknowledgements The authors would like to thank the guest editors of this special issue of Tectonophysics for the invitation to contribute and the opportunity to bring a unified view of the fundamentals of mathematical texture analysis to the attention of a broader audience.
5. Conclusions References Harmonics for the group of rotations SO(3) and the sphere X3, respectively, provide a major mathematical prerequisite for texture analysis. Their application in texture analysis can be largely clarified and turned instructive if they are introduced in the way they are used-as real-valued functions avoiding the classical complex approach featured in textbooks of texture analysis. Their single most important property is that of being a group-theoretical representation. This property is lost if, e.g. the conversion of complex into real harmonics (Bunge, 1969, 1982) is done with less care than actually required. Convergence of the associated harmonic series expansions has never been an issue in texture analysis; it has been taken for granted. However, considerations of their summability provide a complementary introduction of standard distributions as required in texture component fitting. Initially introduced on heuristic grounds to deal with the ‘‘ghost problem’’ and the nonnegativity constraint, the focus here is on their common mathematical basis in terms of harmonics. The results presented here in the context of texture analysis are also fundamental for the analysis of individual crystallographic orientation measurements, e.g. the approximate identities are used in kernel density estimations of the orientation probability density function. The duality of summation processes and spherical singular integrals provides the key of an appropriate re-arrangement of the series leading to wavelet texture
Arnold, K.J., 1941. On spherical probability distributions. PhD thesis, Massachusetts Institute of Technology, MA, USA. Berens, H., Butzer, P.L., Pawelke, S., 1968. Limitierungsverfahren von Reihen mehrdimensionaler Kugelfunktionen und deren Saturationsverhalten. Publ. Res. Inst. Math. Sci., Kyoto Univ., Ser. A 4, 201 – 268. Boogaart, K.G.V.D., 2001. Statistics for individual crystallographic orientation measurements. PhD Thesis, Freiberg. Boogaart, K.G.v.d., Schaeben, H., 1999. Singular value decomposition of the odf – pdf projection operator. In: Szpunar, J.A. (Ed.), Proceedings of the 12th International Conference on Textures of Materials (ICOTOM12). NRC Research Press, Ottawa, Ont., Canada, pp. 1363 – 1368. Bunge, H.J., 1969. Mathematische Methoden der Texturanalyse. Akademie Verlag, Berlin. Bunge, H.J., 1982. Texture Analysis in Materials Science. Butterworths, London. Butzer, P.L., Nessel, R.J., 1971. Fourier Analysis and Approximation. Birkha¨user Verlag, Basel. Butzer, P.L., Stens, R.L., Wehrens, H., 1979. Approximation by algebraic convolution integrals. In: Prolla, J.B. (Ed.), Approximation Theory and Functional Analysis. North Holland Publishing, Amsterdam, pp. 71 – 120. Cerejeiras, P., Schaeben, H., Sommen, F., 2002. The spherical X-ray transform. Math. Methods Appl. Sci. 25, 1493 – 1507. Curtis, C.W., Reiner, J., 1962. Representation Theory of Finite Groups and Associative Algebras. Interscience, New York. Dunkl, C.F., 1966. Operators and harmonic analysis on the sphere. Trans. Am. Math. Soc. 125, 250 – 263. Eschner, T., 1993. Texture analysis by means of model functions. Textures Microstruct. 21, 139 – 146. Eschner, T., 1994. Generalized model function for quantitative texture analysis. In: Bunge, H.J., Siegesmund, S., Skrotzki, W., Weber, K. (Eds.), Textures of Geological Materials. DGM Informationsgesellschaft, Oberursel, pp. 15 – 28.
268
H. Schaeben, K.G. van den Boogaart / Tectonophysics 370 (2003) 253–268
Freeden, W., Gervens, T., Schreiner, M., 1998. Constructive Approximation on the Sphere with Applications to Geomathematics. Clarendon Press, Oxford. ¨ ber Fla¨chen mit lauter geschlossenen geoda¨tischen Funk, P., 1913. U Linien. Math. Ann. 74, 278 – 300. ¨ Uber eine geometrische Anwendung der Funk, P., 1916. U Abel’schen Integralgleichung. Math. Ann. 77, 129 – 135. Gel’fand, I.M., Minlos, R.A., Shapiro, Z.Ya., 1963. Representations of the Rotation and Lorentz Groups and Their Applications. Pergamon, Oxford. Hochstadt, H., 1986. The Functions of Mathematical Physics. Wiley, New York. Ivanova, T.M., Nikolayev, D.I., 2001. New standard function for quantitative texture analysis. Phys. Stat. Solidi, B 228, 825 – 836. Mardia, K.V., Jupp, P.E., 2000. Directional Statistics. Wiley, New York. Matthies, S., 1980. Standard functions in texture analysis. Phys. Stat. Solidi, B 101, K111 – K115. Matthies, S., Vinel, G.W., Helming, K., 1987. Standard Distributions in Texture Analysis, vol. I. Akademie Verlag, Berlin. Matthies, S., Helming, K., Steinkopf, T., Kunze, K., 1988. Standard distributions for the case of fibre textures. Phys. Stat. Solidi, B 150, K1 – K5. Matthies, S., Vinel, G.W., Helming, K., 1990. Standard Distributions in Texture Analysis, vol. III. Akademie Verlag, Berlin. Mhaskar, H.N., Narcowich, F.J., Prestin, J., Ward, J.D., 2001. Polynomial frames on the sphere. Adv. Comput. Math. (in press). Nikolayev, D.I., Savyolova, T.I., 1997. Normal distribution on the rotation group SO(3). Textures Microstruct. 29, 201 – 233. Potts, D., Steidl, G., Tasche, M., 1996. Kernels of spherical harmonics and spherical frames. In: Fontanella, F., Jetter, K., Laurent, P.-J. (Eds.), Advanced Topics in Multivariate Approximation. World Scientific Publications, Singapore, pp. 287 – 301. ¨ ber die Bestimmung von Funktionen durch Radon, J., 1917. U ihre Integralwerte la¨ngs gewisser Mannigfaltigkeiten. Berichte
u¨ber die Verhandlungen der Ko¨niglich Sa¨chsischen Gesellschaft der Wissenschaften zu Leipzig. Math.-Phys. Klasse 69, 262 – 277. Savelova, T.I., [identical with Savyolova], 1984. Distribution functions of grains with respect to orientation in polycrystals and their Gaussian approximations. Ind. Lab. 50, 468 – 474 (translated from Zavod. Lab. 50, 48, 52). Schaeben, H., 1996a. Texture modeling or texture approximation with components represented by the von Mises – Fisher matrix distribution on SO(3) and the Bingham distribution on S4+. J. Appl. Crystallogr. 29, 516 – 525. Schaeben, H., 1996b. A unified view of methods to resolve the inverse problem of texture goniometry. Textures Microstruct. 25, 171 – 181. Schaeben, H., 1997. A simple standard orientation density function: the hyperspherical de la Valle¨e Poussin kernel. Phys. Stat. Solidi, B 200, 367 – 376. Schaeben, H., 1999. The de la Valle´e Poussin standard orientation density function. Proceedings of the International Conference on Neutron Texture and Stress Analysis, Dubna, Russia, June 23 – 26, 1997. Textures Microstruct. 33, 365 – 373 (Special Issue). Schaeben, H., Nikolayev, D.I., 1998. The central limit theorem in texture component fit methods. Acta Appl. Math. 53, 59 – 87. Schaeben, H., Potts, D., Prestin, J., 2001. Spherical wavelets with an application in preferred crystallographic orientation. Proceedings of the Annual Conference of the International Association for Mathematical Geology, Cancun, Mexico, Sep 6 – 12, 2001. Vilenkin, N.J., 1968. Special Functions and the Theory of Group Representations. Am. Math. Soc. Transl. of Math. Monographs 22, Providence, 613 pp. Watson, G.S., 1983. Statistics on Spheres. Wiley, New York. Zheng, Y., Doerschuk, P.C., 2000. Explicit computation of orthonormal symmetrized harmonics with applications to the identity representation of the icosahedral group. SIAM J. Math. Anal. 32 (3), 538 – 554.