Structure tensors

Structure tensors ${\bf a}^{(k)}$ are defined as the mean of the $k$-th repeated outer product of a crystallographic axis with itself, say $\bf n$:

\[ {\bf a}^{(k)} := \langle {\bf n}^k \rangle. \]

In the case of a discrete ensemble of $N$ grains with slip plane normals ${\bf n}_i$ for $i=1,\cdots,N$, the structure tensors are

\[ {\bf a}^{(k)} = \dfrac{\sum_i w_i {\bf n}_i^k}{\sum_i w_i}, \]

where $w_i$ is the weight attributed to the $i$-th grain, typically its size.

Alternatively, if the distribution function $n(\theta,\phi)$ of ${\bf n}$-axes is known, the structure tensors are

\[ {\bf a}^{(k)} = \dfrac{\int_{S^2} n \hat{{\bf r}}^k \, \mathrm{d}\Omega}{\int_{S^2} n \, \mathrm{d}\Omega} , \]

where $\mathrm{d}\Omega = \sin(\theta) \mathrm{d}\theta \mathrm{d}\phi$ is the infinitesimal solid angle and $\hat{{\bf r}}(\theta,\phi)$ is the radial unit vector.

Principal frame

Since $n(\theta,\phi)$ and $b(\theta,\phi)$ are antipodally symmetric, odd moments (odd $k$) vanish identically. Hence, ${\bf a}^{(2)}(n)$ and ${\bf a}^{(2)}(b)$ measure the variance of $\bf n$- and $\bf b$-axes, respectively, around the three coordinate axes.

Posing ${\bf a}^{(2)}$ in its principal frame

\[ {\bf a}^{(2)} = \left[\begin{matrix} \lambda_1 & 0 & 0\\ 0 & \lambda_2 & 0\\ 0 & 0 & \lambda_3\\ \end{matrix}\right] \]

therefore has a similar interpretation as in PCA: the first principal component (eigenvector ${\bf m}_1$) is the direction that maximizes the variance (eigenvalue $\lambda_1$) of the projected data (red curve), the second component is the direction orthogonal to the first component that maximizes the variance of the projected data, and so on with the third component.

Convert to spectral

Converting between spectral and tensorial representations is a linear problem in the sense that

\[ {\bf a}^{(k)} = {\bf f}(\hat{n}_2^{m}, \hat{n}_4^{m}, \cdots, \hat{n}_k^{m}) , \qquad\text{(for all $m$)} \]

where ${\bf f}$ is linear in its arguments and $\hat{n}_l^m = n_l^m/n_0^0$.

In the case of ${\bf a}^{(2)}$, the relation is simple:

\[ {\bf a}^{(2)} = \frac{{\bf I}}{3} + \sqrt{\frac{2}{15}} \left[\begin{matrix} \operatorname{Re}[\hat{n}_2^2] - \dfrac{1}{2}\sqrt{\dfrac{2}{3}} \hat{n}_2^0 & -\operatorname{Im}[\hat{n}_2^2] & -\operatorname{Re}[\hat{n}_2^1] \\ & -\operatorname{Re}[\hat{n}_2^2] - \dfrac{1}{2}\sqrt{\dfrac{2}{3}} \hat{n}_2^0 & \operatorname{Im}[\hat{n}_2^1] \\ \mathrm{sym.} & & \sqrt{\dfrac{2}{3}} \hat{n}_2^0 \end{matrix}\right] , \]

but for higher-order structure tensors the expressions are long (not shown). The above applies to $b(\theta,\phi)$ as well.

The following code example shows how to convert between the two CPO representations:

import numpy as np
from specfabpy import specfab as sf
lm, nlm_len = sf.init(8)
nlm = np.zeros((nlm_len), dtype=np.complex64) # state vector (expansion coefficients)

"""
a2 to nlm
"""
a2 = np.diag([0.0,0.25,0.75]) # some second-order structure tensor
nlm[:sf.L2len] = sf.a2_to_nlm(a2) # l<=2 expansion coefficients
a2 = sf.a2(nlm) # nlm back to a2
print('a2 is: ', a2)

"""
a4 to nlm
"""
p = np.array([0,0,1]) # unidirectional CPO
a4 = np.einsum('i,j,k,l', p,p,p,p) # = deltafunction(r-p) 
nlm[:sf.L4len] = sf.a4_to_nlm(a4) # l<=4 expansion coefficients
a4 = sf.a4(nlm) # nlm back to a4 
print('a4 is: ', a4)

"""
a6 to nlm
"""
p = np.array([0,0,1]) # unidirectional CPO
a6 = np.einsum('i,j,k,l,m,n', p,p,p,p,p,p) # = deltafunction(r-p) 
nlm[:sf.L6len] = sf.a6_to_nlm(a6) # l<=6 expansion coefficients
a6 = sf.a6(nlm) # nlm back to a6
print('a6 is: ', a6)

Construct from measurements

The harmonic expansion coefficients of any CPO can be determined from discrete measurements of its crystallographic axes. This requires constructing the corresponding structure tensors (of each crystallographic axis), from which the expansion coefficients may be derived.

In the case of glacier ice where ${\bf n} = {\bf c}$, this can be done as follows:

import numpy as np
from specfabpy import specfab as sf
from specfabpy import discrete as sfdsc
lm, nlm_len = sf.init(8) 

"""
Generate sample of N grains with random c-axis orientations and sizes
"""

N     = 1000 # number of grains
lon   = np.random.uniform(0, 2*np.pi, N)       # uniform longitude
colat = np.arccos(np.random.uniform(-1, 1, N)) # uniform colatitude
caxes = sfdsc.sph2cart(colat, lon)
weights = np.random.uniform(0.5, 1, N) # grain weights (typically grain sizes)
weights /= np.linalg.norm(weights)     # normalize weights

"""
Determine state vector nlm (array of harmonic expansion coefficients)
"""

### Do it all yourself

a6 = np.zeros((3,3,3,3,3,3))
for c in caxes:
    a6 += np.einsum('i,j,k,l,m,n', c,c,c,c,c,c) # sixth outer product of c-axis with itself
a6 /= N # normalize by number of grains if grain sizes are not known
nlm = np.zeros((nlm_len), dtype=np.complex64) # state vector (expansion coefficients)
nlm[:sf.L6len] = sf.a6_to_nlm(a6) # l<=6 expansion coefficients
print(nlm)

### Alternatively, use built-in function ri_to_nlm(caxes, weights, L)

weights = np.ones(N)/N # equal wight to each grain, should sum to 1
nlm[:sf.L6len] = sf.ri_to_nlm(caxes, weights, 6) # l<=6 expansion coefficients
print(nlm)

Note that constructing ${\bf a}^{(6)}$ is to be preferred over ${\bf a}^{(4)}$ and ${\bf a}^{(2)}$, since it contains more information on the fine-scale structure of the distribution: ${\bf a}^{(6)}$ encodes information about the expansion coefficients $l\leq 6$, ${\bf a}^{(4)}$ about $l\leq 4$, and ${\bf a}^{(2)}$ about $l\leq 2$.