A Novel Multi-Scale Particle Morphology Descriptor with the Application of SPHERICAL Harmonics

Particle morphology is of great significance to the grain- and macro-scale behaviors of granular soils. Most existing traditional morphology descriptors have three perennial limitations, i.e., dissensus of definition, inter-scale effect, and surface roughness heterogeneity, which limit the accurate representation of particle morphology. The inter-scale effect refers to the inaccurate representation of the morphological features at the target relative length scale (RLS, i.e., length scale with respective to particle size) caused by the inclusion of additional morphological details existing at other RLS. To effectively eliminate the inter-scale effect and reflect surface roughness heterogeneity, a novel spherical harmonic-based multi-scale morphology descriptor Rinc is proposed to depict the incremental morphology variation (IMV) at different RLS. The following conclusions were drawn: (1) the IMV at each RLS decreases with decreasing RLS while the corresponding particle surface is, in general, getting rougher; (2) artificial neural network (ANN)-based mean impact values (MIVs) of Rinc at different RLS are calculated and the results prove the effective elimination of inter-scale effects by using Rinc; (3) Rinc shows a positive correlation with the rate of increase of surface area RSA at all RLS; (4) Rinc can be utilized to quantify the irregularity and roughness; (5) the surface morphology of a given particle shows different morphology variation in different sections, as well as different variation trends at different RLS. With the capability of eliminating the existing limitations of traditional morphology descriptors, the novel multi-scale descriptor proposed in this paper is very suitable for acting as a morphological gene to represent the multi-scale feature of particle morphology.


Introduction
Particle morphology, an intrinsic characteristic for granular soils, plays a significant role in the grain-scale and, consequently, macro-scale mechanical behaviors of granular soils. Many research findings, either experimental or numerical, proved that the fabric evolution and fabric anisotropy in grain-scale, as well as the compressibility, crushability, shear strength, dilatancy and small strain stiffness in macro-scale, can be highly influenced by the morphology features of granular particles [1][2][3][4][5][6][7][8]. Therefore, the precise and quantitative representation of particle morphology is the prerequisite of further geological and geo-mechanical investigations of granular soils.
To estimate particle morphology like sphericity and roundness, researchers, in early years, compared the microscopic view of particles with the standard charts developed by Krumbein and Sloss [9]. Although quite convenient for a small number of particles, this chart comparing method is very subjective and is not suitable for a population of particles. In recent years, the development of optical equipment and imaging techniques boosted the accuracy and efficiency of representation

Data Acquisition by X-ray μCT
X-ray μCT as a powerful tool to visualize and characterize the grain-scale mechanical behavior of granular soils was widely adopted in many researches, including grain morphology [12,22,25], grain-scale kinematics [26][27][28][29], and fabric evolution [30,31]. Due to its three-dimensional, highresolution, and non-destructive merits, X-ray μCT is utilized in this paper to acquire data of particle morphology. Specifically, LBS with a diameter range of 400-800 μm are utilized to form a dry cylindrical sample with a diameter and a height of 8 mm and 16 mm, respectively. Then, the sample is scanned within a synchrotron-based μCT scanner at the BL13W beamline of the Shanghai Synchrotron Radiation Facility (SSRF). An X-ray energy of 25 KeV and voxel size of 6.5 μm are used for the CT scanning, since this setting can provide a good contrast between solid particles and void spaces [31].
Once a raw CT image of the sample is acquired, a series of image processing steps and analyses are generally required to extract the surface of each particle. The image processing procedure, in general, is carried out in five steps. Firstly, a series of raw projection images obtained from synchrotron μCT scans at different rotation angles is transformed to gray-scale CT slices on free software PITRE [32] from SSRF [33]. Secondly, an anisotropic diffusion filter method [34] is conducted on those gray-scale CT slices obtained to remove unexpected noise. This filter technique can maintain the original boundaries and enhance the contrast among different phases while dramatically eliminating the noise in the background. Thirdly, an intensity value threshold [35] is utilized on filtered CT slices to convert them into binary images. Fourthly, a marker-based watershed algorithm is adopted on binary images to extract individual particles and store them in a three-dimensional labeled image (Figure 1a). Lastly, the intrinsic MATLAB function bwprim is applied to extract particle boundary voxels for each particle in the LBS sample (Figure 1b), while regionprops3 is utilized to calculate volume and surface area of each particle. A more detailed description of the experimental device and the image processing can be found in Cheng and Wang [31], in which the sample data were used to study the fabric evolution of granular soils under shear.
In order to eliminate the over-segmentation effect, particles with average values of the major and the minor principal axis lengths less than 0.4 mm are removed, and the remaining 4155 LBS particles are extracted for further analysis.

Spherical Harmonic-Based Particle Surface Decomposition
Based on the extracted CT data of individual particles, the particle surface can be decomposed by SHA into a series of sub-surfaces at different RLS. Here, we briefly introduce the SHA procedure for particle surface decomposition. A more detailed description can be found in Zhou et al. [20]. By In order to eliminate the over-segmentation effect, particles with average values of the major and the minor principal axis lengths less than 0.4 mm are removed, and the remaining 4155 LBS particles are extracted for further analysis.

Spherical Harmonic-Based Particle Surface Decomposition
Based on the extracted CT data of individual particles, the particle surface can be decomposed by SHA into a series of sub-surfaces at different RLS. Here, we briefly introduce the SHA procedure for particle surface decomposition. A more detailed description can be found in Zhou et al. [20]. By adopting SH functions, the polar radius of a unit sphere can be extended in different frequencies to match the particle surface points, with different frequencies relating to different RLS. The main functions are shown as follows: where r(θ, ϕ) is the polar radius from the particle center, C m n is the spherical harmonic coefficient, and Y m n (θ, ϕ) is the spherical harmonic function as given by Equation (2).
where n and m are the frequency and the order of the associated Legendre function P m n (x), Based on the spherical harmonic coefficient C m n , the given particle surface can be reconstructed using Equation (1). In addition, the second-order norm of these coefficients expressed by Equation (4) reflects the energy contained in each frequency [36] and, hence, can be utilized to partition and group RLS range.
where * 2 is the second-order norm, and * is the conjugate transpose. Note that L n is an SH invariant with respect to particle translation and rotation. Since L 0 is related to particle size, all the second-order norms are divided by L 0 to eliminate the particle size effect. L 1 is only related to particle shift and, hence, ignored in the classification. In addition, SHA was already proven to accurately describe the morphology information of LR and ST when the maximum SH frequency is 15 or greater [22,[37][38][39] and, hence, SH frequency n = 15 is chosen in this paper for further analysis. Figure 2 depicts the multi-scale features of particle morphology and the averaged second-order norms for all the 4155 particles in each frequency. Based on different L n /L 0 values in each frequency, particle surface can be decomposed into a series of sub-surfaces containing morphological genes at three different RLS, i.e., general form (GF) at n = 2 to 4, local roundness (LR) at n = 5 to 8, and surface texture (ST) n = 9 to 15 [23,[39][40][41]. In order to discuss the inter-scale effect at finer RLS, ST is further divided into first-level surface texture (1L-ST) n = 9 to 12 and second-level surface texture (2L-ST) n = 13 to 15. Therefore, with C m n calculated and Equation (1), any given particle surface can be decomposed into a series of sub-surfaces at different RLS. For instance, the corresponding LR-level sub-surface can be reconstructed by SH frequency n = 8 with the first 81 coefficients of C m n .

Verification
The precise SH reconstruction of particle surface from μCT data is the premise of spherical harmonic-based decomposition. For the verification, particle volumes and surface areas obtained by μCT and SHA are compared using the Z value (Equation (5) in a statistical approach, t-test [42]. In detail, the calculation of particle volume and surface area for μCT data is conducted using MATLAB functions regionprops3 and isosurface, while these values for SHA data are evaluated using the sum of the micro-surface areas of all faces and the sum of the micro-volumes of all the tetrahedrons, respectively (detailed calculation procedure given by Zhou et al. [20]). The Z value, as shown in Equation (5), is related to means (μ) and standard deviations (σ) of two distributions, and it can be utilized to evaluate the divergence of SHA reconstruction from μCT data.
where μ1 and μ2 are the means of two distributions, while σ1 and σ2 are the standard deviations. For all 4155 LBS particles, the mean and standard deviation of particle volume from μCT data are 0.1136 and 0.0454, while those from SHA are 0.1132 and 0.0453, respectively. Similarly, the mean and standard deviation of particle surface area are 1.5419 and 0.4375 for μCT data, and 1.5923 and 0.5203 for SHA. Figure 3 compares their cumulative frequency distributions from μCT and SHA. Solid lines show the volume and surface area distributions from μCT scans, while dotted lines are from SHA. The Z value is 0.0062 for particle volume and 0.0741 for particle surface area, and they are all within 1.96, which reflects a sufficiently close match between these two distributions with a confidence level of 95% [43].

Verification
The precise SH reconstruction of particle surface from µCT data is the premise of spherical harmonic-based decomposition. For the verification, particle volumes and surface areas obtained by µCT and SHA are compared using the Z value (Equation (5) in a statistical approach, t-test [42]. In detail, the calculation of particle volume and surface area for µCT data is conducted using MATLAB functions regionprops3 and isosurface, while these values for SHA data are evaluated using the sum of the micro-surface areas of all faces and the sum of the micro-volumes of all the tetrahedrons, respectively (detailed calculation procedure given by Zhou et al. [20]). The Z value, as shown in Equation (5), is related to means (µ) and standard deviations (σ) of two distributions, and it can be utilized to evaluate the divergence of SHA reconstruction from µCT data.
where µ 1 and µ 2 are the means of two distributions, while σ 1 and σ 2 are the standard deviations. For all 4155 LBS particles, the mean and standard deviation of particle volume from µCT data are 0.1136 and 0.0454, while those from SHA are 0.1132 and 0.0453, respectively. Similarly, the mean and standard deviation of particle surface area are 1.5419 and 0.4375 for µCT data, and 1.5923 and 0.5203 for SHA. Figure 3 compares their cumulative frequency distributions from µCT and SHA. Solid lines show the volume and surface area distributions from µCT scans, while dotted lines are from SHA. The Z value is 0.0062 for particle volume and 0.0741 for particle surface area, and they are all within 1.96, which reflects a sufficiently close match between these two distributions with a confidence level of 95% [43].

Inter-Scale Effect of Traditional Morphology Descriptors
The traditional morphology descriptors listed in Table 1 are selected at different target RLS to illustrate their sensitivity to the SH decomposition of particle surface and to further discuss their inter-scale effects. The definitions of these descriptors are based on the recommendations from Blott and Pye [21]. As mentioned above, most existing particle morphology descriptors can be classified by the target RLS into four different groups, namely, GF, LR, ST, and overall shape parameters. To

Inter-Scale Effect of Traditional Morphology Descriptors
The traditional morphology descriptors listed in Table 1 are selected at different target RLS to illustrate their sensitivity to the SH decomposition of particle surface and to further discuss their inter-scale effects. The definitions of these descriptors are based on the recommendations from Blott and Pye [21]. As mentioned above, most existing particle morphology descriptors can be classified by the target RLS into four different groups, namely, GF, LR, ST, and overall shape parameters. To be more exact, (a) GF, such as aspect ratio, is related to the three principle dimensions of a granular particle, (b) LR is considered to be a measure of the sharpness of corners and edges for a particle, (c) average texture (AT) is the averaged absolute distance of surface points to the mean surface, which can show the morphology variation of ST in a very general way, and (d) overall shape parameters are those that can be influenced independently by features at different RLS. Based on SHA, particle surface is decomposed into a series of sub-surfaces at different RLS. Then, the chosen traditional morphology descriptors are calculated on each sub-surface and this procedure is repeated for all the 4155 LBS particles. Figure 4 shows the variation of particle morphology descriptors with increasing SH frequency. It is found that all the descriptors show a similar trend of decreasing values with an increasing SH frequency except the average texture, which has an opposite trend. This result shows that the "observation length scale" has a significant effect on the estimation of traditional morphology descriptors, proving the existence of the inter-scale effect. The arithmetic average of the target surface departure from the mean surface [46] Overall shape parameter Sphericity S = 3 36πV 2

SA
Ratio of the surface area (S A ) of a sphere with the same volume (V) as the given particle to surface area of this particle [20,47] Ratio of volume to surface area V SA Ratio of particle volume to particle surface area [48] Convexity C x = V VCH Ratio of particle volume over its convex hull volume (V CH ) [49] However, a big difference in the degree of variation exists in different descriptors, with the largest and smallest variations seen in roundness and V/S ratio, respectively. To be more exact, the aspect ratio ( Figure 4a) can be affected by LR-level morphology information, and it shows little sensitivity to high-level morphologies like ST. The LR, from Figure 4b, is highly sensitive to the morphology variation of 1L-ST and 2L-ST and shows a decreasing trend with an increasing SH frequency. The average texture for the ST-level morphology information, as shown in Figure 4c, is sensitive to the smaller RLS such as 2L-ST and witnesses a positive correlation with the SH frequency. The remaining three overall shape parameters (Figure 4d,e,f) reflect the general variation of particle shape. Thus, morphology information from all the different RLS can affect the estimation of these descriptors. Therefore, the inter-scale effect exists in particle morphology representation on all length scales from GF to high-level ST.
frequency. The average texture for the ST-level morphology information, as shown in Figure 4c, is sensitive to the smaller RLS such as 2L-ST and witnesses a positive correlation with the SH frequency. The remaining three overall shape parameters (Figure 4d,e,f) reflect the general variation of particle shape. Thus, morphology information from all the different RLS can affect the estimation of these descriptors. Therefore, the inter-scale effect exists in particle morphology representation on all length scales from GF to high-level ST.

A Novel Multi-Scale Morphology Descriptor
To effectively eliminate the inter-scale effect on representing particle morphology at different RLS and to reflect the surface roughness heterogeneity, a novel spherical harmonic-based multi-scale roughness descriptor is proposed in this paper.

Definition
For any point on a given particle surface, its normal vectors on reconstructed sub-surfaces with different RLS can be calculated. Then, the angle difference of normal vectors on decomposed subsurfaces at target and preceding RLS can reflect the morphology variation at this target RLS. For example, the angle difference between normal vectors of sub-surfaces at LR and 1L-ST can depict the

A Novel Multi-Scale Morphology Descriptor
To effectively eliminate the inter-scale effect on representing particle morphology at different RLS and to reflect the surface roughness heterogeneity, a novel spherical harmonic-based multi-scale roughness descriptor is proposed in this paper.

Definition
For any point on a given particle surface, its normal vectors on reconstructed sub-surfaces with different RLS can be calculated. Then, the angle difference of normal vectors on decomposed sub-surfaces at target and preceding RLS can reflect the morphology variation at this target RLS. For example, the angle difference between normal vectors of sub-surfaces at LR and 1L-ST can depict the morphology variation of 1L-ST at a given surface point. The calculation of this angle difference at a given point, as illustrated in Figure 5, is given by Equation (6).
where n v is the normal vector at the point on the target sub-surface and n vp is the normal vector on the preceding sub-surface. n v and n vp can be determined based on the normal vectors of all the triangulated faces sharing the same vertex (i.e., the target point) on the sub-surfaces as follows: where n f i is the normal vector of face i, α i is the angle-related weight for this face i (i.e., the angle of the vertex on face i), and k is the number of faces that share the vertex (k = 4 in this case). Therefore, based on the angle difference * i Δθ between normal vectors on two sub-surfaces, the descriptor Rinc can depict the incremental morphology variation (IMV) and is defined as follows: where NS is the number of vertices after filtering, nS is the number of vertices on the target surface region, and

* min
Δθ is the filter threshold based on the instrumental resolution (i.e., 6.5 μm in this study) to eliminate the non-physical angle differences introduced by the algorithm, as given in Equation (9). Therefore, based on the angle difference ∆θ * i between normal vectors on two sub-surfaces, the descriptor R inc can depict the incremental morphology variation (IMV) and is defined as follows: where N S is the number of vertices after filtering, n S is the number of vertices on the target surface region, and ∆θ * min is the filter threshold based on the instrumental resolution (i.e., 6.5 µm in this study) to eliminate the non-physical angle differences introduced by the algorithm, as given in Equation (9).
where P Lmin is the minimum pixel length (6.5 µm in this study), and r i is the polar radius of the target point. The value of R inc varies from 0 to 1, and a larger value denotes greater morphology variation. Specifically, 0 denotes that no IMV occurs at the corresponding RLS, and 1 depicts an overall angle variation of π/2 for all the surface points, which, however, cannot be reached for LBS particles due to the surface continuity.

Variation of R inc with SH Decomposition
Since R inc is estimated by the angle difference between normal vectors of the target sub-surface and the corresponding preceding sub-surface, it can be utilized to depict the IMV of a given particle at each RLS and, hence, can act as a morphological gene for the given particle. Different parts of this gene represent morphology variation at different RLS. Table 2 illustrates the variation of R inc with SH decomposition for a sphere, a cube, and two LBS particles. For a sphere, R inc is equal to 0 at all RLS because no IMV occurs across different RLS. Unlike the sphere, R inc values for the cube and two particles show IMV at all RLS. Since R inc at SH frequency n = 4 shows the morphology variation from a sphere to sub-surface of n = 4, it can be utilized to reflect the IMV of GF. From Table 2, Particle 0036 is more visually irregular than Particle 0006, and it shows a larger R inc at n = 4 which means a larger IMV needed from a sphere to the target GF. Therefore, R inc at n = 4 can reflect the irregularity of GF for a given particle. Specifically, the cube has large corners which produce large angle differences and, hence, the cube has an overall larger R inc than these two LBS particles. Moreover, R inc at n = 8, n = 12, and n = 15 depict the IMV of LR, 1L-ST, and 2L-ST, respectively. Figure 6 illustrates the cumulative distribution of R inc at different RLS for all the 4155 LBS particles. By using R inc , each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of R inc is limited to n = 15, i.e., the 2L-ST. However, the R inc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.

Estimating Rinc Using Artificial Neural Network (ANN)
To verify the elimination of inter-scale effects, ANN is introduced to correlate Rinc with the SHbased invariants, based on which the mean impact value (MIV) is then calculated. MIV was first proposed by Dombi et al. [50] to reflect the variation of weighted matrix in ANN. The MIV of Rinc at different RLS can reflect the contribution of morphology variation at each RLS to the estimation of Rinc and, hence, can be utilized to verify the effective elimination of inter-scale effects by Rinc. To calculate the MIV of Rinc, an ANN model is established with Levenberg-Marquardt (LM) algorithm [51]. The input and output parameters are the second-order norms of SH frequencies of 1 to 15 and Rinc at different RLS, respectively. The size of the hidden layer is set to be two-thirds of the sum of input and output parameters [52] to accelerate the training and avoid over-fitting. Table 3 depicts the training performance of ANN with mean squared error (MSE) and regression R value. As can be seen  Sphere cumulative distribution of Rinc at different RLS for all the 4155 LBS particles. By using Rinc, each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.  (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.  (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is stil getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc i limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.  (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.   Figure 6 illustrates the cumulative distribution of Rinc at different RLS for all the 4155 LBS particles. By using Rinc, each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.   Figure 6 illustrates the cumulative distribution of Rinc at different RLS for all the 4155 LBS particles. By using Rinc, each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.   Figure 6 illustrates the cumulative distribution of Rinc at different RLS for all the 4155 LBS particles. By using Rinc, each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.   Figure 6 illustrates the cumulative distribution of Rinc at different RLS for all the 4155 LBS particles. By using Rinc, each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution. = 12, and n = 15 depict the IMV of LR, 1L-ST, and 2L-ST, respectively. Figure 6 illustrates the cumulative distribution of Rinc at different RLS for all the 4155 LBS particles. By using Rinc, each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.   Figure 6 illustrates the cumulative distribution of Rinc at different RLS for all the 4155 LBS particles. By using Rinc, each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.   Figure 6 illustrates the cumulative distribution of Rinc at different RLS for all the 4155 LBS particles. By using Rinc, each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.   Figure 6 illustrates the cumulative distribution of Rinc at different RLS for all the 4155 LBS particles. By using Rinc, each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution. = 12, and n = 15 depict the IMV of LR, 1L-ST, and 2L-ST, respectively. Figure 6 illustrates the cumulative distribution of Rinc at different RLS for all the 4155 LBS particles. By using Rinc, each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.   Figure 6 illustrates the cumulative distribution of Rinc at different RLS for all the 4155 LBS particles. By using Rinc, each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.   Figure 6 illustrates the cumulative distribution of Rinc at different RLS for all the 4155 LBS particles. By using Rinc, each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.   Figure 6 illustrates the cumulative distribution of Rinc at different RLS for all the 4155 LBS particles. By using Rinc, each part on the morphological gene corresponds to different target RLS with a uniform format of definition (Equation (8)). From the figure, the IMV at different RLS decreases with the increasing SH frequency, meaning that the IMV at small RLS is lower than that at large RLS, while the overall surface is still getting rougher with decreasing RLS. Due to the limitation of CT resolution, the discussion of Rinc is limited to n = 15, i.e., the 2L-ST. However, the Rinc proposed can be utilized to represent IMV at any smaller RLS, as long as particle data provided can meet that resolution.

Estimating R inc Using Artificial Neural Network (ANN)
To verify the elimination of inter-scale effects, ANN is introduced to correlate R inc with the SH-based invariants, based on which the mean impact value (MIV) is then calculated. MIV was first proposed by Dombi et al. [50] to reflect the variation of weighted matrix in ANN. The MIV of R inc at different RLS can reflect the contribution of morphology variation at each RLS to the estimation of R inc and, hence, can be utilized to verify the effective elimination of inter-scale effects by R inc . To calculate the MIV of R inc , an ANN model is established with Levenberg-Marquardt (LM) algorithm [51]. The input and output parameters are the second-order norms of SH frequencies of 1 to 15 and R inc at different RLS, respectively. The size of the hidden layer is set to be two-thirds of the sum of input and output parameters [52] to accelerate the training and avoid over-fitting. Table 3 depicts the training performance of ANN with mean squared error (MSE) and regression R value. As can be seen from the table, all MSEs are extremely close to 0 while the R values are greater than 0.97, which reflect the feasibility and accuracy of the ANN model established. In addition, another ANN model is trained with the same input parameters as but different V/S ratio to the output parameter. Since V/S ratio is one of the overall shape parameters that can be independently influenced by morphology information at all RLS, comparisons with MIV results of V/S ratio can further verify the elimination of inter-scale effects. Based on the ANN model established, each item of the input second-order norms is increased and decreased by 10% to obtain two outputs, and this procedure is repeated for all the samples. The difference between the two output datasets obtained is the impact value (IV). Then, the MIV is calculated by averaging the IV for all samples, as expressed in Equation (10). The absolute MIV represents the relative significance of independent variables with respect to dependent variables, and a greater MIV value reflects a higher influence.
where X i × 110% means the i-th item of input data increases by 10% while the others remain unchanged, is the output by ANN, and n is the number of samples. Figure 7 shows the MIV results of R inc and V/S ratio. Solid lines in Figure 7a illustrate the R inc at different RLS from GF (n = 4) to 2L-ST (n = 15). The abscissa denotes the SH frequencies which correspond to different RLS, while the ordinate depicts the MIV that can show the significance of each SH frequency in the estimation of R inc . As can be seen from the figure, R inc at different RLS shows a generally higher level of sensitivity to the corresponding RLS without inter-scale effects from other RLS. For example, R inc at n = 4, based on the angle difference between normal vectors of sub-surfaces at n = 0 and n = 4, is highly sensitive to the SH frequency range from n = 2 to n = 4 and, thus, dominates the GF of the given particle. Similarly, R inc at n = 8 dominates the LR; R inc at n = 12 and n = 15 dominate the 1L-ST and 2L-ST, respectively. Furthermore, it is found that the estimation of R inc will be slightly affected by GF, and this influence becomes lower for the reconstructed particle surface with a higher n value. Since, for a given sand particle, the GF is the same for all reconstructed particles with varying n values, the GF-related effects can be ignored. The MIVs of V/S ratio in Figure 7b, on the other hand, show an overall low sensitivity without any dominant RLS, which means that morphology information on all RLS can contribute to the evaluation of V/S ratio with similar weights. Therefore, by reflecting the incremental morphology variation at the target RLS without inter-scale effects from other RLS, the R inc proposed is very suitable for acting as a morphological gene to represent the multi-scale feature of particle morphology. where Xi × 110% means the i-th item of input data increases by 10% while the others remain unchanged, is the output by ANN, and n is the number of samples. Figure 7 shows the MIV results of Rinc and V/S ratio. Solid lines in Figure 7a illustrate the Rinc at different RLS from GF (n = 4) to 2L-ST (n = 15). The abscissa denotes the SH frequencies which correspond to different RLS, while the ordinate depicts the MIV that can show the significance of each SH frequency in the estimation of Rinc. As can be seen from the figure, Rinc at different RLS shows a generally higher level of sensitivity to the corresponding RLS without inter-scale effects from other RLS. For example, Rinc at n = 4, based on the angle difference between normal vectors of sub-surfaces at n = 0 and n = 4, is highly sensitive to the SH frequency range from n = 2 to n = 4 and, thus, dominates the GF of the given particle. Similarly, Rinc at n = 8 dominates the LR; Rinc at n = 12 and n = 15 dominate the 1L-ST and 2L-ST, respectively. Furthermore, it is found that the estimation of Rinc will be slightly affected by GF, and this influence becomes lower for the reconstructed particle surface with a higher n value. Since, for a given sand particle, the GF is the same for all reconstructed particles with varying n values, the GF-related effects can be ignored. The MIVs of V/S ratio in Figure 7b, on the other hand, show an overall low sensitivity without any dominant RLS, which means that morphology information on all RLS can contribute to the evaluation of V/S ratio with similar weights. Therefore, by reflecting the incremental morphology variation at the target RLS without inter-scale effects from other RLS, the Rinc proposed is very suitable for acting as a morphological gene to represent the multiscale feature of particle morphology.

Variation of Rinc against Incremental Surface Area
It is interesting to examine the influence of IMV represented by Rinc on the increase of particle surface area. The rate of the increase of surface area at the target RLS for a given particle can be expressed as follows:

Variation of R inc against Incremental Surface Area
It is interesting to examine the influence of IMV represented by R inc on the increase of particle surface area. The rate of the increase of surface area at the target RLS for a given particle can be expressed as follows: where S Av and S Ap are the surface area of particle at the target RLS and the preceding RLS, respectively. Figure 8 depicts the relationship of R inc and R SA at different RLS. It can be seen that R inc vs. R SA shows a similar pattern of positive correlation for all RLS, and the degree of data scatter actually decreases with the decreasing RLS (note the different scales adopted in different sub-figures of Figure 8). This result clearly suggests that the IMV represented by R inc is the cause for the increase of particle surface area. The correlations in Figure 8 can be fitted using Equation (12).
where α and β are scale-related parameters, and their values for different RLS are listed in the table of the corresponding figure. The correlation indices r 2 for all RLS are larger than 0.83, which reflects the feasibility of Equation (12).

Variation of Rinc against Traditional Descriptors at Target RLS
To help better understand the descriptor Rinc proposed in this paper, we correlate Rinc with traditional descriptors at different target RLS, as shown in Figure 9. As can be seen from the figure, Rinc shows a negative correlation with aspect ratio and sphericity while a positive one with average texture at the corresponding RLS. Specifically, Rinc at n = 4, as discussed above, reflects the IMV of GF, and a larger Rinc at n = 4 leads to a smaller aspect ratio. The intercept on the ordinate axis is around 1 which means that Rinc at n = 4 equal to 0 depicts a sphere. In addition, sphericity shows a similar variation trend with aspect ratio. This intercept on the ordinate axis is also around 1 which depicts that, for a sphere, Rinc = 0 at all RLS. Average textures of both 1L-ST and 2L-ST, on the other hand, illustrate an opposite variation trend, i.e., increasing with increasing Rinc at the corresponding RLS. Furthermore, from Figure 9b, roundness is poorly correlated to Rinc at n = 8. This is because the estimation of roundness is based on the curvature radius of all corners, as well as the maximum inscribed circle [20,53,54] and, hence, roundness focuses on the convex part for a given particle. Unlike roundness, Rinc at n = 8 reflects the IMV from GF to LR, and it contains morphology variation of both convex and concave parts for a given particle.

Variation of R inc against Traditional Descriptors at Target RLS
To help better understand the descriptor R inc proposed in this paper, we correlate R inc with traditional descriptors at different target RLS, as shown in Figure 9. As can be seen from the figure, R inc shows a negative correlation with aspect ratio and sphericity while a positive one with average texture at the corresponding RLS. Specifically, R inc at n = 4, as discussed above, reflects the IMV of GF, and a larger R inc at n = 4 leads to a smaller aspect ratio. The intercept on the ordinate axis is around 1 which means that R inc at n = 4 equal to 0 depicts a sphere. In addition, sphericity shows a similar variation trend with aspect ratio. This intercept on the ordinate axis is also around 1 which depicts that, for a sphere, R inc = 0 at all RLS. Average textures of both 1L-ST and 2L-ST, on the other hand, illustrate an opposite variation trend, i.e., increasing with increasing R inc at the corresponding RLS. Furthermore, from Figure 9b, roundness is poorly correlated to R inc at n = 8. This is because the estimation of roundness is based on the curvature radius of all corners, as well as the maximum inscribed circle [20,53,54] and, hence, roundness focuses on the convex part for a given particle. Unlike roundness, R inc at n = 8 reflects the IMV from GF to LR, and it contains morphology variation of both convex and concave parts for a given particle.

Surface Roughness Heterogeneity by Rinc
Obviously, Rinc (Equation (8)) can be utilized to evaluate morphology variation for any local region of the particle surface at a target RLS. This value represents the IMV at an individual point i when nS = 1, while it represents the IMV of the whole particle when the target surface region is equal

Surface Roughness Heterogeneity by R inc
Obviously, R inc (Equation (8)) can be utilized to evaluate morphology variation for any local region of the particle surface at a target RLS. This value represents the IMV at an individual point i when n S = 1, while it represents the IMV of the whole particle when the target surface region is equal to the entire particle surface.
Three local points P1, P2, and P3 are selected randomly from a given particle surface, and eight parts of the surface are obtained by dividing the given particle with the three Cartesian coordination planes across the particle center (i.e., the XOY, YOZ, and XOZ planes), as shown in Figure 10, are used to illustrate the surface roughness heterogeneity. Table 4 lists the surface roughness heterogeneity quantified by different R inc values at different RLS. For the whole particle, the R inc values are 0.2003 for GF, 0.1303 for LR, and 0.1051 and 0.0658 for 1L-ST and 2L-ST. However, the three local points show different R inc values at all RLS, reflecting the roughness heterogeneity of particle morphology. A similar phenomenon can be seen from the R inc of all the local surfaces at different RLS. Specifically, the local surface S1 shows a higher R inc than S2 for the GF, LR, and 2L-ST, but a lower value for the 1L-ST. This observation implies that the R inc at different RLS are independent of each other, meaning that a given surface may be smooth at one RLS but rough at other RLS. Specifically, the local surface S1 shows a higher Rinc than S2 for the GF, LR, and 2L-ST, but a lower value for the 1L-ST. This observation implies that the Rinc at different RLS are independent of each other, meaning that a given surface may be smooth at one RLS but rough at other RLS.

Conclusions
A combined X-ray μCT and SHA technique was utilized to decompose the given particle surface into a series of sub-surfaces at different RLS. A total number of 4155 LBS particles were used to create a large dataset of particle surfaces with different target RLS.
Four groups of traditional morphology descriptors were selected at different RLS to investigate the effect of SH decomposition and the inter-scale effect. It was found that the inter-scale effect, either from large scale or from small scale, affects the estimation of particle morphology at all RLS.
To effectively eliminate the inter-scale effect and precisely represent the morphology variation of a given particle at different RLS, a novel spherical harmonic-based multi-scale morphology

Conclusions
A combined X-ray µCT and SHA technique was utilized to decompose the given particle surface into a series of sub-surfaces at different RLS. A total number of 4155 LBS particles were used to create a large dataset of particle surfaces with different target RLS.
Four groups of traditional morphology descriptors were selected at different RLS to investigate the effect of SH decomposition and the inter-scale effect. It was found that the inter-scale effect, either from large scale or from small scale, affects the estimation of particle morphology at all RLS.
To effectively eliminate the inter-scale effect and precisely represent the morphology variation of a given particle at different RLS, a novel spherical harmonic-based multi-scale morphology descriptor R inc was proposed. It is concluded that, firstly, by investigating the variation of R inc against RLS, the particle surface is, in general, rougher with decreasing RLS, but the IMV at each RLS shows a decreasing trend. Secondly, ANN-based MIVs of R inc at different RLS are calculated and the results show a general high level of sensitivity to the corresponding RLS without inter-scale effects from other RLS, which proves the effective elimination of the inter-scale effect. Thirdly, by introducing the increase rate of surface area R SA , the R inc proposed shows a positive correlation with R SA at all RLS, and it can be expressed by R inc = α + β × ln(R SA ), where α and β are scale-related parameters. Fourthly, by correlating R inc with traditional descriptors at target RLS, R inc at the corresponding RLS can be utilized as an alternative method to quantify the irregularity and roughness for a given particle. Lastly, the surface roughness heterogeneity was investigated using R inc , and it was found that the surface morphology of a given particle shows different IMV in different sections, as well as different variation trends at different RLS. Therefore, with the merits of (1) having a uniform format of definition across all RLS, (2) effectively eliminating the inter-scale effects, and (3) reflecting the surface roughness heterogeneity, the R inc proposed in this paper is very suitable for acting as a morphological gene to represent the multi-scale feature of particle morphology. Artificial