Computationally efficient approach to quantify 2 D particle morphological descriptors

The macro-scale behavior of granular materials is strongly influenced by grain kinematics. The mobility of the grains in turn is affected by grain morphology which needs to be comprehensively characterized. Initially, morphological descriptors were determined by manual processes which were tedious and cumbersome. But with the help of image processing techniques and computational geometry, this problem can be handled with ease. The long-established descriptors such as sphericity and roundness can be estimated by operating over the grain boundary obtained from 2D images of the particle. The main objective of this study is to quantify these descriptors in a computationally effective way. The roundness quantification in previous studies involves the sequence of processes such as removing the noise, corner identification and fitting circles. This paper details the necessary modifications to the quantification process required to reduce the cost of time for a single particle. Further, the influence of different smoothing techniques and a new corner identification method will be detailed.


Introduction
Realistic quantification of an ensemble/continuum level granular material response through grain to grain interactions is a complex problem under varied boundary conditions. The morphological features of grains plays an important role in predicting the material response. Macroscale experimental studies showed that the particle shape affects the minimum and maximum packing configuration (e max & e min ), stiffness, critical state friction angle (φ c ) and dilatancy (ψ) [1,2]. Thus in order to evaluate the particle morphology, several descriptions are available in the literature [3,4]. In this study, sphericity (particle form) and roundness (particle corners) based on Wadell [5] will be used.
The descriptors usually distinguish the morphology by elucidating how deviant it is from a perfect circle and hence these processes involves fitting circles. Earlier this process of circle fitting is either manually etched or by visual inspection from charts which can lead to an under or over prediction of the true value. However in the past decade, using image processing and computational geometry, an accurate description can be obtained within a relatively short time. Zheng & Hryciw [6] used these techniques and provided accurate quantification of particle morphology. An efficient way of estimating the morphological characteristics is needed especially when used on a granular assembly. The current study focuses on modifications to the geometric techniques in order to improve the accuracy at a reduced computational cost. * e-mail: ce19s020@smail.iitm.ac.in * * e-mail: rameshkk@iitm.ac.in

Roundness
Roundness represents the sharpness of the particle corners. Wadell [5] defines roundness as the ratio of average radius of curvature of the corners to the maximum inscribed circle. The corresponding equation 1 is as follows.
where r i is the radius of circle fitted to the i th corner, n is the total number of corners while r ins is the maximum inscribed circle radius. The computation of roundness from a 2D image of the particle involves the following steps -Extraction of the raw boundary of the particle -Removing noise from the raw boundary -Convex corners identification on the smooth boundary -Fitting circles to the identified convex corners The boundary is extracted from a Scanning Electron Microscopy (SEM) image of a sand particle (Figure 1a). The extracted boundary consists of sharp crest and trough, this could be due to the particle surface texture or electronic noise. The noise removal is essential to compute roundness. Hence, the closed raw boundary is converted to an open profile (Figure 1b) by finding the radial distances at every small θ interval from the center of the closed curve. The ordinates of the open profile are expressed as, where y i are the raw coordinates (radial distance), f (x i ) are the noise-free boundary values at x i (angle) and i is the noise. The noise-free values can be estimated by choosing an appropriate non-parametric smoothing technique. Some of these techniques and their influence over computational time are discussed in the ensuing.

Locally weighted regression smoothing
Locally weighted scatter plot (LOESS) is a weighted least square (equation 3) regression technique which is used to find the smoothed value f (x i ) at any x i , for a selected local neighborhood and fitting a two-degree polynomial. The weight w for any point x k in the neighborhood is given by equation 4. w controls the influence of the neighborhood points (x k ) over the local fitting of the polynomial. If the distance between the x k and x i is less then w(x) increases and vice versa. The weight function helps to capture the local curvature while the neighborhood is defined by the span parameter (β). β is the ratio between the number of local points that are considered for fitting a polynomial to the total data points and it ranges from 0 to 1.

smoothing splines
The smoothing spline is a piece-wise continuous nonparametric polynomial function [7] having continuous derivatives. Like LOESS, smoothing splines can also estimate the noise-free points, but the latter gives a smooth function instead of a value. Smoothing spline regression is treated as penalty based least-squares method. The induced penalty controls the level of smoothness. The smooth values are obtained from the spline function S(x) that minimizes the equation 5. The second part of the equation 5 is the noise controller and λ is called a smoothing parameter which varies from 0 to 1.

Gaussian Process Regression (GPR)
Gaussian process regression is a Bayesian approach to the regression technique. In this method, the noise i is assumed to be following a normal distribution with mean 0 and variance σ 2 noise . The smooth term f (x) is a random variable distributed as a Gaussian process. Any set of random variables following a joint Gaussian distribution is said to be a Gaussian process (equation 6). The mean function and covariance function define the Gaussian process. A mean of 0 and squared exponential covariance function (equation 7) are used in this study. The process of obtaining f (x) is explained in [8].
In equations (7 & 8), l is the length scale parameter, σ 2 f is the variance of f (x) , δ pq is the Kronecker delta and x p , x q are input values.
The parameters (β-LOESS, λ-smoothing splines, l, σ 2 f , σ 2 noise -GPR) affects the estimation of smoothed profile. Hence, these are to be optimized such that there is no under-fitting or over-fitting. The process of optimizing parameters is termed as hyperparmeter tuning. In this study, cross-validation technique [6,8] is used to for tuning the parameters. This is implemented with MATLAB's inbuilt functions for cross-validation. The raw open profile shown in figure 1b is smoothed by using these three smoothing techniques with appropriate optimum parameters ( Figure 2). LOESS and smoothing spline methods perform the task faster than Gaussian process regression as the time taken for optimization is more in the case of GPR as it has more parameters.

Corner identification using α-shape
After obtaining the smooth open profile, it is then converted back to a closed boundary. The smooth boundary shown in figure 3 is attained from LOESS method. A corner in the closed boundary is identified as a convex corner if any two points on the corner are connected, and this line must lie inside the boundary. Zheng & Hryciw [6] identified convex corners using a two-step process. First, by discretising key points from all boundary points. Then these key points are classified as a convex corner or concave corner points.
In this study, the convex corner points are identified by constructing α-shape [9] using the boundary points. α shapes are widely used for surface reconstruction of given geometry [10]. An intuitive description of α-shape states that a huge chunk of ice-cream is lying in space, and the chocolate pieces are the set of points in it. A sphere-shaped spoon is used to scoop out the ice-cream and the shape formed is α-shape. α is the radius of the spoon used, and in 2D, this is the radius of the circle. The α-shape is a generalization of the convex hull construction technique. A convex hull is the largest convex polygon bounding a set of points with the outer points as its vertices. By α-shape method, the convex hull can be shrunk to capture some of the interior points as well. The radius which captures all the boundary points is the critical radius (α cr ). At a large radius (α ∞ ), the α-shape is the convex hull. α-shape varies from convex hull to original boundary, and this can be controlled by a shrink factor.
The shrink factor (S) is the ratio of α cr to α. When α is chosen as α ∞ then S=0 and when α is α cr then S=1. The red dotted line leaping around the smooth boundary  figure 3 is the α-shape constructed with shrink factor as 0.1. Wherever the α shape coincides with smooth boundary, those points are the convex corner points. The black circular points are the convex corner points and are 155 in number. Some of the corner points are marked in red in order to distinguish two adjacent corners.

Fitting circles to corners
Convex corners with more than two points are fitted with circles. Each corner is separated, and circles are fitted based on least squares approach [6]. The maximum inscribed circle is obtained from the binary image of the particle (akin to figure 1a). The nearest distance to the boundary from each pixel inside the boundary is calculated and maximum of these distances will be the radius of the maximum inscribed circle (r ins ). The black dotted line inside the boundary in the figure 4 is the maximum inscribed circle. Roundness is then computed from equation 1 as 0.23.

Sphericity
According to Wadell [5], 2D sphericity is the ratio of radius of maximum inscribed circle to the radius of minimum circumscribed circle. The maximum inscribed circle is already obtained as part of roundness computation. The minimum circumscribed circle is obtained by calculating the distance between boundary points and the maximum of this will be the diameter of the minimum circumscribed circle. The sphericity of the sand particle ( Figure 4) is calculated to be 0.65.

Effect of shrink factor and total computational time
The roundness is computed for a complex shaped sand particle given in [11] ( Figure 5). The particle boundary  is smoothed by LOESS, and its shown in blue. At shrink factor (S) of 0.2, some of the corners are not identified and at other corners, the identified points are not enough to fit circles and these corners are pointed out by arrows ( Figure  5a). When the shrink factor increased to 0.45, the α-shape captures all the potential corner points (Figure 5b). The roundness computed for this complex shape is 0.13. As more points will be captured when S is increased, the time taken for circle fitting also increases. If S is greater than 0.5, the α-shape starts capturing the concave parts as well.
The three smoothing techniques are applied on the particle taken from [11] (Figure 5), and their effect on computational time is given in table 1. Among the smoothing methods discussed above, GPR consumes more time due to more number of parameters that need to be optimized. The computational time for corner finding using α-shape method is more efficient compared to Zheng & Hyrciw [6] as shown in table 1. While using Zheng & Hyrciw method [6], a maximum divergence value δ 0 of 0.1 is used for finding the corner points. Using LOESS or smoothing spline method for noise removal, α-shape for corner identification and separation would reduce the computational time significantly.

Conclusion
The morphological descriptors such as roundness is calculated accurately and efficiently using an appropriate computational geometric approach. After extracting the 2D boundary of the particle, three smoothing techniques are implemented to remove the noise found on the boundaries. It was observed that the LOESS and smoothing spline methods are effective compared to GPR. Further, to identify the corners present in a closed curve, utilization of α-shape is considered less expensive in terms of computational cost than the boundary discretisation technique. The shrink factor plays an important role in capturing the interior corner points. The shrink factor (S) is usually kept under 0.5, this is because at higher values the α-shape resembles the boundary. The usage of α-shape method for computing roundness is validated using a complex particle shape provided in [11].