Spherical harmonic mathematical analysis

The procedure to generate the 2-D surface of a 3-D particle is qualitatively the same as in generating the 1-D surface of a 2-D particle. Starting from the center of mass, which is defined in 3-D similarly to 2-D, line segments are positioned from the center of mass to the surface at various angles (_i, $(\theta_i,\phi_j)$_j), where these angles are the usual spherical polar coordinates [36]. The length of these line segments, denoted R_ij(_i , _j), are found numerically at the values of the angles corresponding to the points of a double Gaussian quadrature, one for each angle, where R_ij is the distance from the center of mass to a surface point along the direction defined by ($(\theta_i,\phi_j)$_i, _j). A surface point is defined as a point on the surface of a pixel. The Cartesian coordinates of these surface points, with the origin defined at the particle center of mass, are

$$R_{ij} \sin(\theta_i)\cos(\phi_j)$$ x_ij =

$$R_{ij} \sin(\theta_i)\sin(\phi_j)$$ y_ij = (8)

$$R_{ij} \cos(\theta_i)$$ z_ij =

Now is the time to state more clearly what kinds of particles can be handled accurately by this kind of analysis. We require that the shape to be "star-like" [20,21,22]. For a shape to be star-like, any line segment whose one endpoint is the center of mass and the other endpoint is on the surface must be totally contained in the shape itself. Therefore, there can be no "overhangs" or "bubbles" in the shape, as was depicted in Fig. 2. This requirement should be satisfied by almost all aggregates, because the grinding process, whether natural, in a riverbed, or artificial, in a rock-crushing machine, should break down any "overhanging" bits, and rocks usually do not have internal porosity large enough to show up in a tomograph that uses the resolutions of 10's of micrometers per pixel. Many aggregates do have internal porosity, but as long as it is smaller than the resolution of the tomograph, it will not be seen in the particle images created by the tomograph. If there are any of these features in the x-ray tomographic images, the spherical harmonic analysis process will create a "valley" from the internal bubble to the real surface, since encountering the bubble surface first as a line segment proceeds from the center of mass will cause it to be interpreted as the real particle surface. The signature of this happening will be a significant difference between the digital volume and the volume as computed from the spherical harmonic expansion. All the examples of aggregates that will be shown in this paper are star-like, as judged from the original tomographic image, and did not have internal porosity.

Once a star-like particle has been obtained, and surface points R_ij($R_{ij}(\theta_i,\phi_j)$, _j) found, then spherical harmonic analysis (the 3-D equivalent of 2-D Fourier analysis) can be applied. The key equation in the spherical harmonic analysis process is the following, where r(, $r(\theta,\phi)$_j) is any smooth function defined on the unit sphere (0 , 0 2$0 \leq \theta \leq \pi, 0 \leq \phi \leq 2\pi$) [36,37,38]:

$$r(\theta,\phi) = \sum_{n=0}^{\infty} \sum_{m=-n}^n a_{nm} Y_n^m(\theta,\phi)$$ (9)

For these assumptions, the spherical expansion exists and converges. In our case, r(, $r(\theta,\phi)$) is given numerically by (_i , _j , R_ij ). The function is called a spherical harmonic function, and is given by:

$$Y_n^m(\theta,\phi) = \sqrt{\left (\frac{(2 n +1) (n-m)!} {4 \pi (n+m)!} \right )} P_n^m(\,\cos(\theta)\,)\, e^{im\phi}$$ (10)

The functions P_n^m(x) are called associated Legendre functions, and are a set of orthogonal polynomials found in quantum mechanics [37] and many other fields. Appendix A lists the associated Legendre functions up to order n = 8. In this case, x = cos(). Values of higher order associated Legendre functions can be found using recursion relations [36]. These recursion relations are available in user-ready Fortran programs like DXLEGF, a part of the SLATEC numerical package [39]. Using explicit formulae for the associated Legendre functions up to n = 8 helps give more accuracy to this recursion process, hence the listing herein. Using explicit formulas up to higher values of n would be still more helpful, but the algebra to calculate these quickly becomes tedious.

The computed surface points are then used to calculate the coefficients a_nm, which depend on both n and m according to the following definition:

$$a_{nm} = \int_0^{2 \pi} \int_0^{\pi} d\phi \: d\theta \: \sin(\theta) \: r(\theta,\phi) \: Y^{m*}_n$$ (11)

where the asterisk denotes the complex conjugate. Choosing each angle to correspond to the points of a Gaussian quadrature makes evaluation of these integrals straightforward. In 3-D, a 120 point Gaussian quadrature was used for each angle, so that eq. (11) was evaluated by summing over 120²=14,400 points. In some cases, a 240 point Gaussian quadrature, with 57,600 points, was used.

A set of coefficients, once determined, then serve as a complete, within numerical error, mathematical characterization of the aggregate particle. Much of the later sections of this paper examines error analysis of how well the expansion works for simple shapes, by direct numerical comparison to analytically known quantities, and by visual and numerical comparison to the original random digital particles from the tomograph.

Many properties of the shape can be computed once the spherical harmonic expansion is known. These include volume, surface area, mean and Gaussian curvature both at a point of the surface and integrated over the surface, and the moment of inertia tensor, as defined below. The volume and moment of inertia can also be computed directly from the digital image by counting voxels.

The equation for the volume of the shape in polar coordinates is particularly simple and is given by:

$$V = \int_0^{2\pi} \int_0^{\pi} \int_0^{r(\theta,\phi)} r^2 \sin(\theta) dr \: d\theta \: d\phi$$ (12)

where the integral is over all angles and for values of r between the origin at the center of mass and the surface, r(, $r(\theta,\phi)$), which is given by the computed expansion of eq.(9). The r integral can be analytically performed, and the resulting integral is then, in terms of the function r(, $r(\theta,\phi)$),

$$V = \frac{1}{3} \int_0^{2\pi} \int_0^{\pi} r^3(\theta,\phi) \sin(\theta) d\theta d\phi.$$ (13)

The equations for the surface area and integrated curvatures involve some auxiliary terms that are defined in differential geometry. These are given below [22]. A useful way of denoting points on the surface of the particle is by the vector $\vec{X}$, which is the Cartesian coordinates of surface points. The components of $\vec{X}$ (X₁ = x, X₂ = y, X₃ = z) are similar to those given in eq. (8), and derivatives of $\vec{X}$ are denoted by a subscript.

There are 8 auxiliary quantities that are useful in analyzing surfaces, and they are all built out of components of $\vec{X}$ and the surface normal vector $\hat{n}$, which is given by:

$$\hat{n}=\frac{\vec{X}_{\theta} \times \vec{X}_{\phi}}{ \left \vert\vec{X}_{\theta} \times \vec{X}_{\phi} \right \vert}$$ (14)

The differential surface area element, d , which is the area of the patch of surface at r(, $r(\theta,\phi)$), is given by:

$$d\Lambda = S d\theta d\phi \: , \: S = \left \vert\vec{X}_{\theta} \times \vec{X}_{\phi} \right \vert$$ (15)

The parameters E, F, and G are given by

$$E=\vec{X}_{\theta} \cdot \vec{X}_{\theta}$$

$$F=\vec{X}_{\theta} \cdot \vec{X}_{\phi}$$ (16)

$$G=\vec{X}_{\phi} \cdot \vec{X}_{\phi}$$

The parameters L, M, and N are

$$L=-\vec{X}_{\theta} \cdot \hat{n}_{\theta}$$

$$M=\frac{1}{2} \left (\vec{X}_{\theta} \cdot \hat{n}_{\phi} + \vec{X}_{\phi} \cdot \hat{n}_{\theta} \right)$$ (17)

$$N=-\vec{X}_{\phi} \cdot \hat{n}_{\phi}$$

Since the function r(, $r(\theta,\phi)$) is known in terms of a spherical harmonic expansion, and the Cartesian coordinates of the surface are known in terms of r(, $r(\theta,\phi)$), all the above derivatives can be taken and the results for all the various surface parameters given in terms of derivatives of r(, $r(\theta,\phi)$). In the following, these derivatives of r(, $r(\theta,\phi)$) are denoted by the notation:

$$r = r(\theta,\phi)$$

$$\frac{\partial r(\theta,\phi)}{\partial \theta} = r_{\theta}$$

$$\frac{\partial r(\theta,\phi)}{\partial \phi} = r_{\phi}$$ (18)

$$\frac{\partial^2 r(\theta,\phi)}{\partial \phi^2} = r_{\phi\phi}$$

$$\frac{\partial^2 r(\theta,\phi)}{\partial \theta^2} = r_{\theta\theta}$$

$$\frac{\partial^2 r(\theta,\phi)}{\partial \phi \partial \theta} = r_{\phi\theta} = r_{\theta \phi}$$

The actual functional form of these derivatives, in terms of the spherical harmonic expansion (eq. (9)) are given in Appendix B.

The surface area S_A is an integral over $\theta$ and $\phi$ of the differential surface element in eq. (15). The parameter S, using the spherical harmonic expansion, is

$$S = r [r_{\phi}^2+r_{\theta}^2 \sin^2(\theta) +r^2 \sin^2(\theta)]^{1/2}$$ (19)

so that the surface area S_A is

$$S_A = \int_0^{2\pi} \int_0^{\pi} r [r_{\phi}^2+r_{\theta}^2 \sin^2(\theta) +r^2 \sin^2(\theta)]^{1/2} d\theta d\phi$$ (20)

The parameters E, F, and G are given by

$$E=r^2_{\theta}+r^2$$

$$F=r_{\theta} r_{\phi}$$ (21)

$$G=r^2_{\phi}+r^2 \sin^2{\theta}$$

The parameters L, M, and N involve vector products of the derivatives of the components of $\vec{X}$ and $\hat{n}$ (see eq. (17)). The components of $\hat{n}$ and derivatives of the components are fairly complicated, so that the explicit forms of L, M, and Nare not given in the text. The derivatives of the components of $\vec{X}$ and $\hat{n}$ are listed in Appendix B, however. Equation (17) can then be used to write out the equations for L, M, and N, and numerically evaluate them and the values of other quantities that depend on them.

The local mean curvature H is defined as the arithmetical mean of the two principal curvatures at each point on the surface [22], and is given by

$$H = \frac{EN + GL - 2FM}{2(EG - F^2)} = H(\theta,\phi)$$ (22)

The Gaussian curvature K, which is another measure of surface curvature, is the geometric mean of the two principal curvatures at each point on the surface:

$$K = \frac{LN - M^2}{EG - F^2} = K(\theta,\phi)$$ (23)

The mean curvature, averaged over the surface and weighted by the differential surface element, is defined here as h, where

$$h=\frac{1}{S_A} \: \int_0^{2\pi}\: \int_0^{\pi} \: H \: d\Lambda$$ (24)

The parameter h has units of inverse length, since the mean curvature H has units of inverse length. If h is then inverted, the resulting length depends on the size and shape of the object considered. One should note that the term "integrated mean curvature" is often referred to in the literature as eq. (24) without the normalizing factor of S_A^-1 [40].

An interesting property of the Gaussian curvature, when similarly integrated over the surface, is that

$$k = \frac{1}{4\pi} \: \int_S K d\Lambda = 1$$ (25)

when the object under consideration is topologically equivalent to a sphere. This is the case for the star-like aggregate considered here. "Topogically equivalent to a sphere" means if the object were made out of very pliant rubber, it could be deformed, without ripping or puncturing, into a sphere. Equation (25) is then a very useful quality control check for the spherical harmonic expansion procedure. Also, one should be able to make a judgement of when "enough" terms in the expansion have been computed, by when this criterion is fulfilled.

The terms of the moment of inertia tensor I_ij, where i,j = 1,2,3, are given by [7]

$$I_{ij}= \int_V \rho(\vec{r})\: d^3r \: (\delta_{ij} r^2 - x_i x_j)$$ (26)

where $\rho(\vec{r})$ is the density of the object, $\delta_{ij}$_ij is the Kronecker delta (zero for i $i \neq j$j and 1 when i = j), the integral is over the entire particle volume V, and I_ji = I_ij. If the mass density of the object is uniform, so that $\rho = M/V$ = M/V (total mass divided by total volume), then using the definition of the polar coordinates in eq. (8), and performing the r integral like in eq. (13), the I₁₁ integral, for example, transforms to

$$I_{11}=\frac{M\:}{5V} \int_0^{2\pi} \int_0^{\pi} r^5(\theta,\... ...(\theta)[1 - \sin^2(\theta) \cos^2(\phi)] \: d\theta \: d\phi$$ (27)

The remaining components of the moment of inertia tensor can be found in Appendix B.

Even though the result of eq. (25) exists for any aggregate shape that would be considered, it is important to have other checks as well, to establish the limitations and perform error analysis on the spherical harmonic expansion procedure. This error analysis is carried out next using various analytical shapes.