The two objectives of this paper were to demonstrate use the of the discrete element method for generat ing synthetic images of spherical particle configurations, and to compare the performance of 9 classic feature extraction methods for predict ing the particle siz e distribution s (PSD) from these images . The discrete element code YADE was used to generate synthetic images of granular materials to build the dataset. Nine feature extraction methods were compared: Haralick features, Histogram s of Oriented Gradients, En tropy, Local Binary Pattern s, Local Configuration Pattern, Complete Local Binary Pattern s, the Fast Fourier transform, Gabor filter s, and Discrete Haar Wavelet s. The feature extraction methods were used to generate the inputs of neural networks to predict the PSD. The results show that feature extraction methods can predict the percentage passing with a root -mean -square error (RMSE) on the percentage passing as low as 1.7%. CLBP showed the best result for all particle sizes with a RMSE of 3.8 %. Better RMSE were obtained for the finest sieve (2.1%) compared to coarsest sieve (5.2%). Keywords: image analysis, feature extraction, particle size distribution, discrete elements, neural network * Corresponding author : Laboratory for Geotechnical and Geoenvironmental Engineering (LG2), École de technologie supérieure, 1100 Notre -Dame Ouest, Mo ntreal, Quebec, H3C 1K3, Canada. Email: mjmanashti@gmail.com 1 Introduction Sieving, the current standard method for PSD determination in geotechnical engineering [1], dates back to the first half of the 20th century. With this method, the specimen is poured in a series of sieves and shaken for a predetermined duration [2]. The cu mulative mass passing each sieve is presented against the sieve opening on a logarithmic scale. Several useful parameters can be extracted from the PSD, such as the characteristic sizes corresponding to specific percentages passing (e.g. D50, the particle size for which 50 % of the granular material mass is smaller), and parameters describing the shape of the PSD (e.g. Cu = D60/D10, the coefficient of uniformity). These parameters are used for soil classification and the prediction of soil properties (e.g. [3]). Sieving yields accurate results, but is time -consuming. Several methods have been developed to determine the PSD based on soil photographs and image analysis techniques. These methods can be classified in two main groups. Methods of the first group use image segmentation to define the contour of each particle in the photograph. The size of each segmented zone is then used to estimate the PSD. Methods of the second group extract textural features from soil photographs to quantify the PSD [4]. Commerc ial codes such as WipFrag [5] and VisioRock [6] use image segmentation to determine the particle size distribution of granular materials from photographs. The relationship between the size of the segmented zone in the photographs and the real particle size is not straightforward as segregation, capturing and overlapping errors will cause the segmented zone and particle size distributions to differ [7, 8]. Segregation refers to the reorganization of particles according to their size when they are displaced, for instance when they are pushed down a slope by machinery or when they are shaken. Capturing refers to the unequal probability of particles of different sizes of being observed on a surface. For example, a particle that has a diameter of the same order a s the sample or pile is much more likely to be visible at the surface compared to a small particle. Overlapping errors are due to particle superposition. Particles are hiding each other. Even if some methods are available to take these errors into account, they remain largely unresolved [8]. The accuracy of segmentation methods differs depending on the parameters that are compared, the reference method, and the granular material used in the comparison. The parameter that is most often compared is the D50. For 10 digital images of blasted sandstone, Sudhakar et al. [9] compared automatic segmentation using commercial codes WipFrag and Fragalyst [10], with manual segmentation on tracing paper. They did not compare their PSD with sieving results. Their comparison focused on segmentation errors and avoided the influence of segregation, capturing and overlapping errors. The difference between the D50 of Wi pFrag, Fragalyst, and the manual segmentation was below 10 % for half of the digital images. Errors of up to a factor 2 were observed on some of the other images. Liu and Tran [11] compared the PSD obtained with commercial segmentation codes FragScan [12], WipFrag [5] and Split [13], with sieving results for dry waste rock from a mining operation. The ratio between the D50 for the three commercial codes and the reference D50 obtained by sieving varied between 1.3 and 3.2. On average, Split and WipFrag overe stimated the D50 by around 55%, while the overestimation by FragScan was around 100%. The PSD given by Sudhakar et al. [9] and Liu and Tran [11] also allow the root mean square error (RMSE) on the percentage passing for each size class to be estimated and compared. RMSE values between 13 and 36 % were obtained when comparing FragScan, WipFrag and Split for dry waste rock [11]. When comparing only the segmentation, RMSE of 15 and 20 % were obtained for WipFrag and Fragalyst [9]. Textural features are functions that describe the spatial arrangement of pixel intensities [14, 4]. Methods for the extraction of textural features can be divided in three main categories [15]. Methods of the first type use statistical parameters to descri be the spatial distribution of pixel intensities. For example, the Haralick features [16] use statistics (e.g. average, variance, correlation) of the global co -occurrence matrix that describes the number of pixels of given intensities separated by a given distance. Methods of the second type are based on local patterns found in the image (e.g. local binary pattern (LBP) [17]). The third type includes transform -based feature extraction methods, such as Gabor wavelets [18]. Statistical textural parameters hav e been used for PSD determination or rock classification. Itoh et al. [19] compared the performance of 32 textural and colour features, including Haralick features and run - length statistics, for the prediction of particle size. Five materials with differen t uniform PSD were photographed under 22 different lighting levels for a total of 3960 images. An analysis of variance showed that color properties, like saturation and hue, are sensitive to changing illumination. Haralick features and run -length statistic s applied on edge -enhanced images were found to be better predictor of the particle size. Local pattern methods describe the relationship between the intensity of the gray levels in the neighborhood of each pixel of an image. The Local Binary Pattern (LB P) method is a good example [17]. With this method, the relative intensity of each pixel in the neighborhood of a center pixel is encoded as a binary number. Statistics are then computed on the occurrence of each binary pattern in the image. Local pattern methods can be rotation invariant and independent of illumination [17]. These methods are generally used for classification purposes. For example, Paclík et al. [20] used LBP, among other methods, in a rock classification system. To our best knowledge, lo cal pattern methods have never been used to determine the PSD of granular materials. Transform -based methods allow features to be extracted from images in the frequency domain. The most common method for PSD determination is the Harr wavelet transform [21, 22]. Wavelet transforms give local information on the frequency content of an image. The wavelet transformation of an image involves multiplying the local pixel intensities in an image by a mother wavelet. The parts of the image where the frequency conten t corresponds to the frequency of the wavelet will show a stronger response. Shin and Hryciw [22] compared the accuracy of two -dimensional Haar wavelet and Haralick features for PSD determination. Both wavelet transforms and Haralick features were dependen t on the illumination level and particle colour. However, the normalization of the wavelet energy was shown to eliminate the effect of illumination and particle color. Hryciw et al . [21] also presented the theoretical basis behind the use of wavelet analys es for PSD determinations. Yaghoobi et al. [23] determined the particle size of fragmented rock using features extracted from Fourier transforms, wavelet transforms, Gabor filter and their combinations. With Fourier transforms, the mean and standard deviat ion of the magnitude spectrum for 20 rings of equal thickness was used for each image. The best performances were obtained with Fourier transforms, followed by Gabor filters, and a combination of Fourier and wavelet transforms. Size estimation was better f or fine to medium particles. Several researchers have combined image processing methods with Artificial Neural Networks (ANN) for PSD determination [24, 25, 26] . Because of the structure of classic ANN, it is generally not possible to use the intensity of each pixel directly as ANN inputs. Consequently, preprocessing techniques have been used to extract parameters, such as textural parameters, from the images to be used as ANN inputs. Ghalib et al. [25] used a selection of eight Haralick features as the input for an ANN to predict the average size of sub -angular sand particles. Their dataset included 220 images. The particle size RMSE varied between 1.4 and 6.1 %. In another study, Ko and Shang [27] defined a uniformity coefficient that was combined with an initial PSD estimate obtained with commercial code WipFrag as inputs for an ANN. The ANN output was the particle sizes corresponding to weight percentages passi ng of 50, 75, and 90 % (D50, D75 and D90). The uniformity coefficient was obtained through the thresholding of 5 areas in each photograph (center, top left, right and bottom). The uniformity coefficient was defined as the difference in the proportion of da rker pixels in the center area with respect to the four surrounding areas. Only 9 images were used for the ANN training. The ratio between the D50, D75 and D90 estimated from the neural network and the target values obtained from sieve analyses varied betw een 0.7 and 1.4. Hamzeloo et al. [28] used image processing in combination with an ANN to predict the PSD of crushed rocks on a conveyor belt. They segmented their images manually after several image processing steps (e.g. image sharpening, edge detection, thresholding). They then calculated different geometrical features for the segmented particles (e.g. Ferret diameter). A Principal Component Analysis (PCA) was then used to select the best features and to use them as the input to a series of ANN. They com pared their results with sieve analyses. Their dataset included 21 sieve analyses. Yaghoobi et al. [23] used Fourier transforms, Gabor filters, wavelet transforms, and their combinations to extract image features that were later used as inputs for neural n etworks to determine the PSD of fragmented rock. Their database included 226 photographs. Manual segmentation results obtained with commercial code Split were used as the target for the neural network. The training of neural networks requires large datase ts that are not readily available. The previously cited examples are based on a limited number of images, generally much less than 1000, and they involve a limited number of input parameters. Because they involve different granular materials and different targets for the neural network (i.e. sieving results or manual segmentation results), it is also difficult to compare their results. The main objectives of this paper are 1) to compare the accuracy of various textural parameters in predicting the PSD of ma terials composed of spherical particles and 2) to introduce a new dataset of synthetic images of granular materials generated using the discrete element method. The textural parameters are used as the input for a series of neural networks. The selected par ameters include some that have never been tested for PSD analyses, notably LBP and the related methods, Entropy, and Histograms of Oriented Gradients (HOG). These new methods are compared with more established methods like Haralick parameters, Haar wavelet transforms, Fourier transforms, and Gabor filters. The database presented in this paper contains 53 003 pairs of images (top and bottom viewpoints) of spherical particle configurations created with the DEM. The paper first describes the synthetic dataset, the feature extraction methods and the neural network parameters. The paper then presents the results and compares the accuracy of each method with previously published results. The synthetic dataset is show n to facilitate the comparison of a different image analysis methods. 2 Methodology 2.1 Synthetic granular image dataset A large image dataset was prepared by Pirnia et al. [29] to compare the feature extraction methods and for training, validation, and testing of the neural networks. The discrete element (DE) code YADE was used to generate a large number of synthetic photographs of granular materials. A text file containing a list of 53003 PSD with particle sizes between 75 and 1180 µm was created. The list was obtained by considering all possible weight percentage combinations with 5 % increments for 5 sieves (106, 150, 250, 425 and 710 µm). COMSOL’s Java interface was used to re ad the data and control YADE’s python interface to generate the DE specimen [30]. The particles were created in a transparent virtual box. They were given random colors. The particles were allowed to settle in the box. Pictures were taken using two virtual cameras from the top and bottom of the box ( Fig. 1). For this paper, all pictures were converted to grayscale images. A total of 53003 pairs of top and bottom images sized 128 x 128 pixels were produced. Image pairs were merged into 256 x 128 -pixel images , which were used as the input of the model. Because each particle has a finite mass, the real PSD of each image can differ slightly from the specified PSD. However, the RMSE of the percentage passing for each sieve was found to be smaller than 0.2 %. Fig. 1 Top (upper row, a -c) and bottom (lower row, d -f) views of the synthetic particles in the virtual transparent box. Left (a, d), middle (b, e) and right (c, f) images show examples of coarse, widely graded, and fine PSD 2.2 Feature Extraction methods The following sections describe the feature extraction methods that were applied on the dataset. These feature extraction methods were used to generate features that were fed to a series of artificial neural networks presented later in the paper. 2.2.1 Haralick Features Haralick features are statistics of the gray level cooccurrence matrix (GLCM) [31, 16, 19] . Item ( i, j) in the GLCM gives the number of pixels with intensity level i that have a specified spatial relationship with other pixels with intensity level j. The GLCM has two important parameters: offset and number of gray levels ( Ng). The offset defines the spatial relationship of the two pixels. Horizontal offsets of 1 and 10 pixels ([0 1] and [0 10]) were tested in this study. The number of gray levels determines the size of the GLCM. GLCM were calculated for both 8 and 64 gray levels. The 14 textural features defined by Haralick et al. [16] were calculated for each GLCM (t otal of 56 features) using the MATLAB script of Monzel [32]. Preliminary tests showed that some Haralick features are better correlated with the percentage passing. This is the case of correlation and contrast: 𝐶𝑜𝑟𝑟𝑒𝑙𝑎𝑡𝑖𝑜𝑛 =∑ ∑ 𝑖 𝑗 𝐺𝐿𝐶𝑀 (𝑖,𝑗)−𝜇𝑥𝜇𝑦𝑁𝑔 𝑗=1𝑁𝑔 𝑖=1 𝜎𝑥𝜎𝑦 (1) 𝐶𝑜𝑛𝑡𝑟𝑎𝑠𝑡 =∑ 𝑛2𝑁𝑔−1 𝑛=0∑ ∑ 𝐺𝐿𝐶𝑀 (𝑖,𝑗)𝑁𝑔 𝑖=1𝑁𝑔 𝑗=1 |𝑖−𝑗|=𝑛 (2) where µ𝑥 and σ x are respectively the mean and standard deviation of px(i), and µ𝑦 and σ y are the mean and standard deviation of py(j), respectively : 𝒑𝒙(𝒊)=∑ 𝑮𝑳𝑪𝑴 (𝒊,𝒋)𝑵𝒈 𝒋=𝟏 (3) 𝐩𝐲(𝐣)=∑ 𝐆𝐋𝐂𝐌 (𝐢,𝐣)𝐍𝐠 𝐢=𝟏 (4) Correlation describes the linear dependency in pixel intensity and contrast increases when the image contains abrupt changes in pixel intensity [31]. 2.2.2 Histogram of Oriented Gradients (HOG) Histograms of Oriented Gradients (HOG) describe the image by counting the occurrence of different orientations of the gray intensity gradient for a series of predefined cells [33]. The cells can be square or rectangular. They are defined by their size. The gradients are assigned to bins. The gradient magnitude is used to weight the gradient when they are assigned to the bins: higher gradients have a higher weight. Larger blocks tha t include several cells are also defined to normalize the content of each bin by the norm of all bins in the block. Normalizing allows local contrast variations in an image to be taken into account. Cells that are smaller, or similar in size to the particl es will tend to be dominated by a small number of particle edges that will give a preferential orientation to the gradient. Cells that are larger than the particle size will result in a more random gradient orientation [34]. Fig. 2 shows a comparison of th e HOG presented as rose diagrams for the coarsest PSD and a finer PSD for a cell size of 8 pixels. The gradient orientation is more random for the fine PSD. Fig. 2 Examples of gradient orientations calculated for a cell size of 8 pixels for a) coarse and b) fine PSD In this study, square HOG cells with sizes of 2, 4, 8, 16, 32, and 64 pixels were used. MATLAB function extractHOGFeatures [35] was used to generate the HOG features. The mean and standard deviation of all HOG for each size were calculated to produce 12 features that were used as the ANN input. The normalizing blocks were defined to include 4 cells with an overlap of one cell width for contiguous blocks. Nine orientations bins were used. 2.2.3 Local entropy of grayscale image Entropy [36] is a measurement of the randomness of the gray level intensity in the neighborhood of a pixel. Entropy is a statistical parameter related to the image texture. The entropy of a grayscale image can be calculated as follows: 𝐻(𝑋)= −∑ 𝑃𝑗𝑁𝑔 𝑗=1 𝑙𝑜𝑔2𝑃𝑗 (5) where H(X) is the entropy of the neighborhood X around a center pixel, Ng is the number of gray level bins, and Pj is the probability of occurrence of gray level j in the neighborhood X [37]. The neighborhood around the center pixel is defined usin g a filter. Fig. 3 shows the filters that were used in this study. These filters were defined using the MATLAB functions getnhood and strel. The filter diameters vary between 1 and 9 pixels. The entropy was calculated with MATLAB function entropyfilt with bin numbers varying between 2 and 256. Figure 4 presents the effect of the disk filter size and number of bins on H(X). Coarser particles correspond to lower values of the local entropy. The mean and standard deviation of H(x) were calcula ted for the 9 disk filters (total of 18 features) and used as the ANN input. Fig. 3 Nine disk filter sizes used to calculate entropy Fig. 4 Effect of disk filter size and number of bins on the local entropy for the widely graded PSD (Fig. 1 b and 1e). Coarser particles are associated with lower values of the local entropy (darker gray levels) 2.2.4 Local Binary Pattern (LBP) The Local Binary Pattern (LBP) sets the center pixel of a neighborhood as a threshold [17]. Neighboring pixels take a value of 0 or 1 based on the gray level difference between them and the center pixel [38]. This pattern is encoded into a binary value by multiplying the 0 or 1 values by 2p, where p is the neighboring pixel number. For example, the LBP for 8 neighboring pixels corresponds to an 8 -bit number. LBP is calculated for each pixel of an image with this relationship: 𝐿𝐵𝑃 =∑ 𝑠(𝑥𝑝−𝑥𝑐)𝑝−1 𝑝=0 2𝑝 (6) where xc is the center pixel gray level and xp are the gray levels for the neighboring pixels. s( xp – xc) is the step function: 𝑠(𝑥𝑝−𝑥𝑐)={1,𝑥𝑝−𝑥𝑐≥0 0,𝑥𝑝−𝑥𝑐<0 (7) In its original form, LBP uses the closest 8 pixels around the center pixel. In later LBP implementations, different numbers of neighboring pixels and distances between the center and neighboring pixels can be chosen [39]. To achieve rotation invariance, t he encodings that only differ by a rotation can be combined into bins, thus limiting the number of features. In this project, the LBP features were calculated for radius values of 1, 2, 3, 4, 5, 6, 8, and 12 pixels, and for 4, 6, 8, 12, and 16 neighbors to find the best filters for the fine and coarse particles. 2.2.5 Local configuration pattern (LCP) The Local Configuration Pattern (LCP) combines LBP with microscopic configurations (MiC) [40, 41]. LBP encodes the structure of the neighborhood around a cen ter pixel based on a step function that does not include any information on the variability of the gray level intensity. MiC adds information on this variability. With MiC, the value of the center pixels for each rotation invariant LBP bin is defined as a linear combination of the neighborhood pixel intensities multiplied by a set of weights: [𝑥𝑐,0 ⋮ 𝑥𝑐,𝑁]=[𝑥0,0… 𝑥𝑝,0 ⋮ ⋱ ⋮ 𝑥𝑝,𝑁… 𝑥𝑝,𝑁][𝑎0 ⋮ 𝑎 𝑝] (8) where xc,0 through xc,N are the center pixels in the current LBP bin, x0,0 through xp,N are the neighborhood pixels and a0 through ap are the weight for the neighborhood pixels. The weight are determined using a least -square method. A Fourier transform is used to make the weight s rotation invariant. The LCP features combines the LBP histogram with the MiC weights. In this project, LCP features were calculated with radius of 1, 2, or 3 pixels and neighborhoods of 1, 2, 3, 4, 5, 6, 7, 8, or 12 pixels. 2.2.6 Completed Local Binary P attern (CLBP) The Completed Local Binary Pattern (CLBP) has three outputs: CLBP -Center, CLBP -Sign, and CLBP -Magnitude. CLBP -Sign corresponds to LBP [42]. CLBP -Center describes the gray level of the center pixel. It is converted into a binary value via global thresholdin g using the mean gray level for the image. CLBP -magnitude is calculated using the absolute value of the difference between the gray level of the neighboring pixels and the center pixel ( |xp−xc|). The absolute difference for each neighboring pixel is converte d to a binary value by thresholding using the mean absolute difference for the image. These binary values are encoded in a number similarly to LBP. In this project, the radius was varied between 1 and 3. The number of neighbors ranged from 1 to 6. The MAT LAB codes used for the calculation of LBP [17], LCP [41], and CLBP [43] were obtained from the Center for Machine Vision Research, Department of Computer Science and Engineering, University of Oulu, Finland [44]. 2.2.7 Fourier Transform The two -dimensional discrete Fourier transform 𝐹(𝑢,𝑣) of an image with pixel intensity f( x, y) is defined as follows [45, 23]: 𝐹(𝑢,𝑣)=1 𝑀𝑁∑ ∑ 𝑓(𝑥,𝑦)𝑒−2𝜋√−1(𝑢𝑥 𝑀+𝑣𝑦 𝑁)𝑛−1 𝑦=0𝑀−1 𝑥=0 (9) where M×N is the image size, u and x=0, 1, 2, …, M-1 and v and y=0, 1, 2, …, N-1. The magnitude and standard deviation of the 2D FFT was used to find the PSD. Figure 5 shows the magnitude spectrum of the Fourier Transform for the images shown in Figs. 1a, 1b and 1c. The center of the magnitude spectrum is associated w ith the coarse particles while the surrounding area is associated with the fine particles. The brighter center of the magnitude spectrum in Fig. 5a is associated with coarser particles in the original image. This brighter center can be compared with the d arker center in Fig. 5c associated with fine particles. Fig. 5 Magnitude spectrum of the Fourier transform for images corresponding to the coarse (a), widely graded (b), and fine (c) PSD in Figures 1a, 1b, and 1c, respectively The particle size distribution can be related to the mean value and standard deviation of the Fourier transform magnitude for a set of ring filter s [23]. In this study, the mean value and the standard deviation of the magnitude were calculated for different sets of ring filters. Two parametric studies were conducted. In the first, the number of rings was set to five and the thickness of consecutive rings was multiplied by a constant increment X between 0.8 and 2.4. Two examples are shown in Fig. 6 for X=1.0 (first row) and X=2.0 (second row). In the second method, the thickness of each ring was constant, and the number of ring filters was varied between 5 and 20. The third row in Fig. 6 shows eight rings with a constant thickness. Fig. 6 Examples of concentric circular filte rs for which the mean and standard deviation of the Fourier transform magnitude spectrum is calculated. The first row shows five filters with equal thickness (X=1.0). The second row illustrates five filters in which the thickness of every ring is twice th e previous ring thickness ( X=2.0). The third row shows eight filters with an equal thickness The mean value and standard deviation of the magnitude were used as the input for an ANN. Average values were calculated for the top and bottom images using the same ring filters. The total number of inputs for the ANN is equal to twice the number of rings. It varies between 10 and 40 for the two parametric studies. 2.2.8 Gabor Filters Gabor filters with different wavelengths and orientations were applied to the images to extract features to determine the PSD. Gabor filters combine a Gaussian kernel function with a sinusoidal plane wave [4, 23] . They allow the local frequency content of an image to be evaluated. A two - dimensional Gabor function ( F(x, y)) and its two -dimensional Fourier transform ( F(u, v)) can be defined as follows: 𝐹(𝑥,𝑦)=(1 2𝜋𝑆𝑥𝑆𝑦)𝑒𝑥𝑝 (−1 2(𝑋2 𝑆𝑥2+𝑌2 𝑆𝑦2))𝑐𝑜𝑠(2𝜋𝑢0𝑋) (10) 𝐹(𝑢,𝑣)=1 2{𝑒𝑥𝑝 {−1 2[(𝑢−𝑢0) 𝑆𝑢2+𝑣2 𝑆𝑣2]}+𝑒𝑥𝑝 {−1 2[(𝑢+𝑢0) 𝑆𝑢2+𝑣2 𝑆𝑣2]}} (11) 𝑆𝑢=1 2𝜋𝑆𝑥 ,𝑆𝑣=1 2𝜋𝑆𝑦 (12) where Sx and Sy are the Gaussian function standard deviations in the x and y directions, respectively. Su and Sv are the Gaussian function standard deviation in the u and v directions, respectively, and u0 is the sinusoidal wavelet central frequency in the x direction. Fig. 7 shows the Gabor filter magnitude for the wavelength of 3,9, and 15 pixels and the orientat ion of 0°, 90°, and 150°, for the top images corresponding to the coarse (top row) and fine (second row) PSD. A filter bank corresponding to different wavelengths (3 to 15 pixels) and orientations (0 to 150°) was prepared using MATLAB function Gabor. The m ean and standard deviation of a set of 42 filters was used as the input of the ANN (84 parameters). MATLAB function imgaborfilt was used to apply the filter bank to the images. Fig. 7 Magnitude of 3 selected Gabor filters with different wavelengths and o rientations for the coarse (top row) and fine (second row) PSD 2.2.9 Haar wavelet transforms Wavelet transforms use mother wavelets (Ψ) for the decomposition of a signal. Mother wavelets are functions that are localized in space and in the frequency domain. The Haar mother wavelet can be described as follows [36]: 𝜓𝑢,𝑣={√𝑢 −√𝑢 0 0≤𝑢𝑥−𝑣≤0.5 0.5≤𝑢𝑥−𝑣≤1.0 𝑒𝑙𝑠𝑒𝑤 ℎ𝑒𝑟𝑒 (13) where u is a scaling factor and v is a translation of the mother wavelet on the x axis. Detailed calculations of 2D Haar wavelet transforms of images in the context of PSD determinations are presented by Hryciw et al. (2015). Transforms can be applied successively to 4×4 pixel regions ( I) that are defined as follows [36]: 𝐼= [𝑎𝑏 𝑐𝑑] (14) Applying the transform results in the following four coefficients: 𝐴= 𝑎+𝑏+𝑐+𝑑 2 (15) 𝐻= 𝑎+𝑏−𝑐−𝑑 2 (16) 𝑉= 𝑎−𝑏+𝑐−𝑑 2 (17) 𝐷= −𝑎+𝑏+𝑐−𝑑 2 (18) where A, H, V and D are respectively the approximation, horizontal, vertical and diagonal coefficients. Each application of the wavelet transform decreases the resolution of the image by a factor 2. Transforms can thus be applied 7 times for 128×256 pixel images. At each new level, the transform is applied to the previous approximation matrix. The square of each coefficient can be used to describe the energy of the image. The sum of the energy associated with the A, H, V and D matrices is preserved after each transform. Consequently, each new level results in three independent energy values [46]. For each level, the mean and standard deviation of the energy associated with the H, V and D coefficients were used as features. This results in 49 features for the 7 levels. MATLAB function haart2 was used to compute the coefficients. Figure 8 presents one exam ple of the successive application of Haar wavelet transforms and the resulting approximations. Fig. 8 Seven levels of wavelet approximations for horizontal, diagonal, and vertical orientations 2.3 Artificial Neural Networks The features that were extracted for the 53 003 images were used as the input of a series of neural networks to predict the percentage passing for sieves of 106, 150, 250, 425 and 710 µm ( Fig. 9). Function fitting neural networks were trained with MATLAB to predict the particle size based on the extracted features. For all methods in this study, one hidden layer with 10 neurons was used. A hyperbolic tangent sigmoid was used as the transfer function. Fig. 9 Artificial Neural Network to predict the percentage passing for a series of sieves For optimal training of the neural networks, 70% of the 53 003 images were used for training, 15% for validation, and 15% for testing. The root mean square error (RMSE) on the percentage passing was used to calculate the error associated with the neural network. The same parameter is used in the presentation of the results to compare the performance of ea ch type of features. The Scaled Conjugate Gradient method was used for the network training. Figure 10 illustrates the performance of the neural network during training for a selection of the 618 best features for all methods. The RMSE decreases for the training, test, and validation datasets during the network training until the maximum number of epochs i s reached (5000). Fig. 10 Neural network training performance for a selection of the 618 best features Larger networks were also tested for the wavelet features and a combination of the 618 best features. Similar performances were obtained for larger net works. For example, an ANN with four hidden layers with 20, 11, 15 and 7 neurons improved the RMSE by less than 0.1 % for wavelet transforms. It was decided to use the same architecture for each type of features to facilitate the comparison. Compared to ob ject classification, good performances for smaller networks have been reported for other applications involving textural features [47, 14] . 3 Results 3.1 Entropy Different number of bins were used to calculate the local entropy. Fig. 11 shows the RMSE on the percentage passing for number of bins between 2 and 256. When the RMSE for all sieves is considered, the best performances were obtained for 8 bins (RMSE = 4.8 %). The optimal number of bins can vary for specific sieves. For example, for the finest sieve, the lowest RMSE (3.5 %) was obtained with 4 bins. Except for 256 bins, better performances were obtained systematically for the finest sieve. This was also obs erved for the other types of features. Fig. 11 RMSE on the percentage passing for the fine sieve, coarse sieve, and all sieves using local entropy with different bin sizes 3.2 LBP To find the best LBP parameters, different radius values (1 to 12 pixels) and number of neighbors (4, 6, 8, 12, 16 pixels) were tested. By changing these parameters, the number of LBP features varied between 15 and 243. Figure 12 shows some examples of RSME for the percentage passing for all sieves combined. The two curves show the influence of the number of neighbors and the radius . The radius was set to 3 for the neighbor curve and the number of neighbors was set to 16 for the radius curve. It can be seen that better performances were obtained for higher numbers of neighbors. The radius curve shows that the minimum RMSE was obtaine d for an optimal radius of 3 pixels. Fig. 12 RMSE on the percentage passing for LBP with different radius values and number of neighbor pixels for all sieve sizes. Increasing the number of features through the number of neighbors or a combination of features for different LBP parameters had a relatively small influence on the RMSE. For example, the RMSE decreased from 6.5 to 4.1 % when the number of features was incr eased from 15 (radius of 2 pixels with 4 neighbors) to 3880 (all features combined). On the other hand, a selection of the best features coud lead to optimized results with a limited number of features. For instance, combining the features corresponding to 4 neighbors and radius values from 1 to 3 pixels, a total of only 45 features , led to a RMSE of 5.1 %. This accuracy is similar to the results for 12 neighbors and a radius of 3 pixels (133 features), and 16 neighbors and a radius of 3 pixels (243 feature s). 3.3 LCP and CLBP As with LBP, the RMSE for LCP and CLBP decreased with an increasing number of neighbors. Figure 13 shows the RMSE for all sieve sizes as a function of the number of neighbors ( N) and the radius ( R) for LCP and CLBP. As with LBP, the best results for LCP when considering all sieve sizes were obtained for a radius of 3 pixels. The lowest RMSE (4.8 %) for LCP was obtained for 12 neighbors and a radius of 3 pixels. A RMSE of 3.8 % was obtained for a combination of all features 024681012141618 024681012141618 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 Neighbor in pixelRadius in pixel RMSE (%) Radius Neighbor with 1 to 12 neighbors and radius values of 1 to 3 pixels. For CLBP, the best result were obtained for 5 neighbors and a radius of 2 pixels (RMSE of 5.3 %). A RMSE of 3.8 % was obtained when combining all CLBP features. Fig. 13 RMSE of LCP and CLBP for different numbers of neighbors and radius for all sieve sizes. N refers to the number of neighbors and R to the radius 3.4 Fourier Transform Figure 14 presents the results obtained using the mean and standard deviation of the Fourier transform magnitude spectrum for 5 rings with size multipliers between 0.5 and 2.4. The size multiplier corresponds to the ratio between the radius es of each consecutive ring. When all sieves are considered, t he lowest RMSE was 5.8 % for a size multiplier of 1.6. This ratio is roughly equal to th e ratio between the opening size of consecutive sieves (e.g. 710 µm / 425 µm = 1.67). The best results for the finest sieve (RMSE = 1.9 %) were obtained for a size multiplier of 1. Fig. 14 RMSE results for the mean and standard deviation of the Fourier transform magnitude spectrum of 5 rings with size multipliers from 0.5 to 2.4 The influence of the number of rings was verified for a size multiplier of 1. The number of rings was found to have relatively little influence on the predicted PSD. For instance, when varying the number of rings from 5 to 20, the RMSE only varied between 6.9 and 7.1 %. Using a large number of rings, as Yaghoobi et al. [23] who used 20 equal rings to predict the size of fragmented rock, does not significantly increase the accuracy of the PSD determination. 3.5 Other methods and combination of all features Figure 15 compares the RMSE of each type of features when predicting the percentage passing for the finest and coarsest sieves, and all sieves. For LBP, LCP, and CLBP, the methods with the best results, RMSE are presented for both the best set of parameters and a combination of all features. For each method, better results were obtained for fine particles. The combination of all LCP features (LCP ALL) provided the best results for fine particles with a RM SE of 1.9 %. The combination of all CLBP features (CLBP ALL) provided the best results for coarse particles with a RMSE of 5.2 %. LCP ALL and CLBP ALL both obtained the lowest RMSE for all sieves (3.8 %). The RMSE for the coarsest sieve, finest sieve, and all sieves were respectively 4.7, 1.7, and 3.4 % for a combination of the 618 best features from all methods. This set of features combined Fourier Equal Ring (68), Wavelet (49) , HOG (12), Fourier Ring size Multiplier (20), Haralick (56), CLBP (7), Gabor ( 84), LBP (135), LCP (169), and Entropy (18). Fig. 15 Comparison of the RMSE on the percentage passing for each type of extracted features and some selected combinations. The number of features for each method is given between parentheses Irrespective of the type of textural feature and sieve size, the RMSE of 1.7 to 10 % in Fig. 16 are low compared to the RMSE of the experimental results presented in the literature. For instance, the comparison of FragScan, WipFrag and Split by Liu and Tran [11] resulted in RMSE varying between 13 and 36 % for dry waste rock. RMSE of 15 and 20 % w ere achieved by WipFrag and Fragalyst when comparing automatic segmentation with manual particle contouring [9]. On the other hand, it should be realized that the synthetic soil photographs ( Fig. 1) are idealized and probably much easier to handle than rea l soil photographs with varying illumination, particle colour, and particle shape . Figure 16 compares the real and predicted percentages passing for each image and each sieve for the comb ination of the 618 best features. The overall coefficient of determination ( R2) is 0.99. As noted previously, larger sieve sizes are associated with poorer predictions. The RMSE gradually increases from 1.7 to 4.7 % with increasing sieve sizes. The R2 coefficient also decreases from 0.99 to 0.91 between the 106 and 710 µm sieves. 6.96.66.25.8 5.6 5.35.2 5.14.8 4.8 4.13.8 3.83.410.0 8.8 7.28.1 6.8 6.87.26.66.35.95.4 5.55.24.7 2.43.54.4 2.9 3.05.0 2.94.03.7 3.7 2.41.9 2.11.7 0.02.04.06.08.010.012.0FOURIER EQUAL RING(68)WAVELET (49)HOG (12)FOURIER RING SIZE MULTIPLIER (20)HARALICK (56)CLBP (7)GABOR (84)LBP (135)LCP (169)ENTROPY (18)LBP ALL (3880)LCP ALL (731)CLBP ALL (90)SELECTION OF BEST FEATURES (618) RMSE (%) All Particles Coarsest Sieve (710 µm) Finest Sieve (106 µm) Fig. 16 Comparison of the real and predicted percentages passing for a selection of the 618 best features Figure 17 compares some randomly selected examples of real and predicted PSD for the 618 best features. The percentage passing is plotted on the usual semi log plot. Blue diamonds and lines show the real PSD and red circles and lines indicate the predictions. The RMSE for these examples varied between 1.4 and 5.2 %. Overall, the real and predicted PSD are similar, even for a RMSE of 5.2%. Fig. 17 Examples of real and predicted PSD for the synthetic images obtained with the DEM and the 618 best features The error percentage on D50 was calculated for the combination of 618 features. The D50 for both the predicted and observed PSD were calculated based on the method presented by Thyabat et al. [48]. The D50 error for the combination of all methods is 6.1 %. The predicted and observed D50 values are compared in Fig. 18 . As with the RMSE on the percentages passing, the error on the D50 is low compared to the values presented in the literature for real photographs. For example, Liu and Tran (1996) obtained D50 errors of approximately 55 % with Split and WipFrag. Again, the good performances for this dataset were expected considering its idealized nature. Fig. 18 Prediction and Observation of D50 for the 618 best features In some cases, the predicted D50 were not accurate because of the shape of the PSD. An example of such PSD is shown in Fig. 17 i. The flat part of the curve corresponding to a percentage passing of 50 % implies that the DEM specimen does not contain par ticles between 150 and 425 µm. H alf of the specimen mass corresponds to particles that are larger than 425 µm than while the other half comprises particles that are smaller than 150 µm. This renders the D50 undefined. These cases correspond to the vertical lines in Fig. 18 at observed D50 values corresponding to sieves sizes of 150, 250 and 425 µm. 4 Conclusion This paper introduced a dataset containing 53003 pairs of synthetic photographs of granular material. This dataset was used to train a series of ANN to predict the PSD. A series of textural features were used as the ANN inputs. All features were able to pr edict the PSD separately or in combinations with a good accuracy. Methods based on the concept of Local Binary Pattern (LBP), such as LCP and CLBP, could predict the PSD better than the other methods, and could be used independently. RMSE of 4.1%, 3.8%, an d 3.8% were obtained for LBP, LCP and CLBP, respectively. To our best knowledge, this paper presented the first example of PSD determination with these methods. The best results were obtained with a combination of 618 features. The RMSE for the coarsest si eve, finest sieve, and all sieves were respectively 4.7 %, 1.7 %, and 3.4 %. The relative error on the D50 was 6.1 %. These errors are small compared to previously published results due to the idealized nature of the images in the dataset. The next step i s to test these methods with real datasets. Other applications can also be envisioned. For example, similar methods could be used for soil classification (e.g., ASTM D2487, ASTM International 2020), and the estimation of soil properties (e.g., shear streng th, compressibility, compaction, hydraulic conductivity). Synthetic image datasets like the one that was introduced in this paper offer several advantages. As they are easier to generate than real photograph datasets, they allow various parameters to be t ested rapidly. For example, the influence of lighting on the performance of different methods could be tested by adding a randomly positioned light source before the image rendering. Synthetic images could also be used to increase the size of real photogra ph datasets. Acknowledgements The authors gratefully acknowledge the support of Hydro -Québec and NSERC for this project.