Comparison of Seven Texture Analysis Indices for Their Applicability to Stereological Correction of Mineral Liberation Assessment in Binary Particle Systems

Abstract

An effective correction method for stereological bias is required because of the importance of accurate assessment of mineral liberation of ore particles. Stereological bias is error caused by the estimation of a three-dimensional liberation state based on two-dimensional sectional measurements. Recent studies have proposed a stereological correction method based on sectional particle texture analysis, which employs numerical particle models. However, the applicability of this method to unfamiliar particle systems, with different shape and texture characteristics from the numerical particle model, has not been thoroughly investigated. In this study, the viability of the method for examination of the internal structure and shape of unfamiliar particles, was assessed using four types of particle systems, based on combinations of two types of internal structures (Boolean and Voronoi) and two types of particle shapes (spherical and irregular). Seven different texture analysis indices were utilized for composition distribution correction with regard to each of the four types of particle systems. The results suggested that a model based on the angular second moment and/or entropy, employed by the spatial gray level dependence method, showed the greatest viability for assessment of unfamiliar particle internal structure and/or shape.

1. Introduction

In mineral processing, comminution mainly aims to achieve mineral liberation, freeing valuable minerals from gangue to enhance the efficiency of the subsequent selection and, consequently, the total efficiency. In the overall assessment of the comminution process, it is also important to assess the mineral liberation state, which is the fraction of liberated and composite particles.

In practice, mineral liberation states are assessed through two-dimensional (2D) measurements of ore particle sections that are resin-mounted, cross-sectioned, and polished. They are typically measured using automated analyzing systems based on scanning electron microscopy/energy-dispersive X-ray spectroscopy (SEM/EDX), such as a mineral liberation analyzer (MLA) [1,2] quantitative evaluation of minerals by scanning electron microscopy [3], and a TESCAN integrated mineral analyzer [4]. These 2D assessment methods inevitably suffer from stereological bias, in which the degree of liberation apparent in 2D is overestimated when compared with the true degree of liberation in 3D.

Several means have been employed to directly measure the 3D liberation state without stereological bias, such as X-ray computed tomography (CT) [5,6,7,8,9,10], serial sectioning [11,12,13], and dense medium separation using heavy liquids [14] or magnetic fluid [15]. However, these all have disadvantages for practical use in terms of analytical speed, cost, ease of performance, safety, and/or identification of minerals with similar specific gravity.

Therefore, a stereological correction method for estimating 3D liberation states from 2D counterparts is required. A simple, elementary method estimates the 3D degree of liberation by dividing the 2D degree of liberation by an empirically obtained parameter [16,17]. Miller and Lin [13] and King and Schneider [15] developed a method for converting 2D composition distributions (corresponding to the composition distributions or liberation spectra in previous studies) into their 3D counterparts by making use of empirically obtained kernel functions. Gay and Morrison [9] estimated a kernel function from the relation between 1D and 2D measurements. Barbery proposed a composition distribution correction method based on a particle structure model incorporating Poisson mosaics and the Boolean method, and a particle shape model employing line assessment [18]. Gay [19] employed a fundamental element of particle section texture analysis by measuring sectional parameters, such as the average distance between two pixel pairs and the average area of three-pixel triangles, with a combination of two phases, and predicted the 3D composition distribution based on stereological equations. However, these methods have not become widely used, mainly because their applicability to unfamiliar particle systems has not been fully validated. The complexity of the internal mineral particle structure has a significant influence on the stereological bias in liberation assessment [11,17], and thus it is vital for any stereological correction method to possess versatile applicability with respect to various types of particles.

Recently, a newly conceptualized stereological correction method based on particle sectional texture analysis was proposed [20,21]. Essentially, the procedure consists of quantifying the mineral texture in the observed particle sections; choosing a numerical binary particle model with the best-fit texture; calculating the stereological bias of the binary particle model; and correcting the stereological bias, with respect to the real ore particles, on the basis of the model’s bias. Figure 1 summarizes the procedure, which includes the following steps:

• visualize the mineral phases in the particle sections, using existing methods such as MLA [1,2] (Figure 1(i));
• binarize the mineral phases with respect to the mineral of interest and others (Figure 1(ii));
• measure the areal fraction of the mineral of interest ($F_a$) and 2D composition distribution (${\Lambda }^{2D}$), and calculate the particle sectional texture characteristics ($f$) through image analysis (Figure 1(iii));
• estimate the stereological correction parameter (${\Lambda }^{dif}$) from the 2D parameters ($f, F_a$), using the corresponding isogram, which was initially obtained by an all-encompassing simulation using binary particle models (Figure 1(iv));
• estimate the 3D composition distribution (${\Lambda }^{3D}$) from ${\Lambda }^{dif}$ and ${\Lambda }^{2D}$ (Figure 1(v)).

A combination of spherical particle shape and spherical mineral phases (sphere-based model) was employed in previous studies [20,22]. However, since real ore never has a simply spherical shape or spherical mineral phases, it is necessary to ensure that any such model is sufficiently robust to apply to various types of particle systems. Critical to this approach is the selection of a proper texture analysis method, which efficiently extracts the texture characteristics most strongly correlated with the stereological bias. In the abovementioned studies, the fractal dimension method (FD) was employed in the texture analysis as a serviceable example [20,22]; however, systematic investigation is required, since there are numerous texture analysis methods, as described below.

Texture analysis methods are grouped into two categories, statistical and structural, with the former more suitable for natural textures and the latter for artificial textures [23]. The former were investigated in the present study, as they were considered more suitable for ore texture analysis. In preliminary testing, we found that the gray level run length method, Fourier power spectrum method, and gray level histogram method were not suitable for incorporation in the proposed stereological correction method, as they are better suited analyze stripe pattern, periodical pattern, and color tone, respectively, which are not common in the ore texture. Therefore, the following three methods were selected as candidates:

• the FD method, based on the fractal dimension ($\delta$) [24] of the surface of an imaginary 3D structure of image intensity [22,25];
• the spatial gray level dependence method (SGLDM), using a co-occurrence matrix [26,27]; the gray level difference method (GLDM), using a gray level difference vector [27].

In order to model the particle internal structure, various 3D structure models have been used for liberation analysis. For the sake of simplicity, a binary structure composed of phases A and B has been assumed, though most of the structure models described below can be extended to multi-phase structures. $F_v$ denotes the volume fraction of phase A to the total volume of the material.

• The dispersed phase model was composed of core A phases dispersed in a B phase matrix. The core phase could be spherical [20,21,22,28,29,30], regularly-shaped [16,30,31], or irregularly-shaped [32]. This differed from the Boolean model (below) because the core phases did not overlap with each other. The maximum $F_v$ did not reach 1.0, and was limited to the packing density of the core phases except when the core phase was aligned in a cubic manner [16,31].
• The Poisson mosaic model was based on Poisson tessellation [33,34]. Small polyhedra composed of flats with random location and orientation were randomly grouped into phase A or B with respect to $F_v$. These were termed Poisson polyhedra in a previous study [35,36]; however, Poisson mosaic, as defined by Gay et al. [18], was used in the present study.
• The Boolean model [34] was an extension of the dispersed phase model. The A phase was dispersed in a B phase matrix, based on a Poisson distribution, irrespective of the previous A phase position. This model is able to produce more complicated phase shapes than the dispersed phase model, though a spherical core phase is commonly used [37,38]. Barbery used the Poisson tessellations generated by the Poisson mosaic model as primary grains, and called the resulting configuration Boolean scheme texture with Poisson polyhedra as primary grains [35,36]. The magnitude of stereological bias in Barbery’s model is reported to be similar to that in the Poisson mosaic model [18].
• The Voronoi model was based on Voronoi tessellation [34,39]. The seeds were randomly distributed, and the space was divided into groups of seeds and polyhedral forms with the greatest relative proximity to one another. The small polyhedra were randomly grouped as phase A or B with respect to $F_v$. Vassiliev et al. used the Voronoi model to simulate ore structure [39]. There are various types of Voronoi model [40], but the most simple and classical Voronoi method was used in this study.
• The Markov random field technique, used to model random textures in 2D [41,42], can be expanded to 3D mineral structures.

In preliminary simulations, we employed the Markov random field technique to model 3D binary structures. However, this resulted in extremely complex structures that could not be compared to the other models. In this study, two texture models (Boolean and Voronoi) were employed as representatives of spherical and irregularly shaped mineral phase models, respectively.

Four types of binary particle systems, combining two types of particle shape with two types of mineral phase, were prepared. One was used for assessing the proposed correction method, and the rest were used as numerical validation samples, one for unfamiliar particle shape and internal structure, one for unfamiliar particle shape alone, and one for unfamiliar internal structure alone. As discussed above, stereological correction was conducted using various types of texture analysis techniques, and their respective viability was assessed. Note that Boolean and Voronoi models are compared to assess the viability of the stereological correction method for unfamiliar particles; however, assessment of the validity of these models themselves, for modeling real ore structure, is beyond the scope of this study.

2. Materials and Methods

2.1. Methodology

2.1.1. Modeling of Binary Particle Systems

Multi-phased particles were the primary targets of the present study. However, given the mineral phase binarization process shown in Figure 1, binary particles composed of phases A and B were modeled, as described below, using the existing procedure [22]. The numerical simulation was conducted with dimensionless parameters, as follows.

• A total of 7463 spherical particles, with diameters ranging from 1 to 2, were generated at random positions in a rectangular prism (30 × 20 × 30), and packed, using the ESyS-Particle open source code (The Centre for Geoscience Computing at the University of Queensland, Brisbane, Australia) [43], with the discrete element method (DEM) [44] (Figure 2a). The particle position data within a prism of (30 × 12 × 30) was digitized into (750 × 300 × 750) voxels, with one voxel representing a volume of 0.04³.
• Three-dimensional structures of phases A and B were generated in a rectangular prism following the rules described in Section 2 below (‘Modeling of Three-Dimensional Phase Structure’) (Figure 2b). The 3D phase data within a (30 × 12 × 30) prism was digitized into (750 × 300 × 750) voxels.
• The internal structures of particles composed of both phase A and B domains were created by comparing the voxel information in steps 1–2 (Figure 2c).

2.1.2. Liberation Analysis

The particle volumes and sectional areas were calculated based on the voxel information described in the previous subsection. Particles whose height at their center of mass was from 2 to 14 were analyzed in 3D. The volume ($V$) and phase A domain ($V_A$) of a given particle were respectively calculated by multiplying the number of voxels constituting the particle and the phase A domain by the single voxel volume ($={0.04}^3$). When the phase A domain was $x (=V_A/V)$, the particle volume was $V\left(x\right)$.

A total of 4316 particle sample sections, with heights of 2, 4, 6, 8, 10, 12, or 14, were analyzed in 2D. The particle section area ($S$) and phase A domain ($S_A$) were respectively calculated by multiplying the number of pixels constituting the particle section and phase A domain by the single pixel area ($={0.04}^2$). The particle sectional area was expressed as $S\left(x\right)$, with the phase A domain corresponding to $x (=S_A/S)$.

The composition distribution was assessed using 12 values of $x$, including $x=0, 0-0.1,$ $0.1-0.2, \dots , 0.9-1 \mathrm{and} 1$, as in previous studies [9,13,15]. The following equations define the volume and areal fractions of phase A ($F_v$, $F_a$), and the 3D and 2D composition distribution indices (${\Lambda }_A^{3D\left(i\right)}, {\Lambda }_A^{2D\left(i\right)}$):

$$F_v=\frac{{{\displaystyle \sum }}_1^M{V}_A}{{{\displaystyle \sum }}_1^{M}V},$$

(1)

$$F_a=\frac{{{\displaystyle \sum }}_1^N{S}_A}{{{\displaystyle \sum }}_1^{N}S},$$

(2)

$${\Lambda }_A^{3D\left(i\right)}=\frac{{\displaystyle \sum }V\left(x\right)}{{{\displaystyle \sum }}_1^{M}V}\mathrm{where}\{\begin{array}{c}i=1 \mathrm{when} x=0.0\\ i=2-10 \mathrm{when} 0.1\left(i-2\right)<x\le 0.1\left(i-1\right)\\ i=11 \mathrm{when} 0.9<x<1.0\\ i=12 \mathrm{when} x=1.0\end{array}$$

(3)

$${\Lambda }_A^{2D\left(i\right)}=\frac{{{\displaystyle \sum }}_1^{N}S\left(x\right)}{{{\displaystyle \sum }}_1^{N}S}\mathrm{where}\{\begin{array}{c}i=1 \mathrm{when} x=0.0\\ i=2-10 \mathrm{when} 0.1\left(i-2\right)<x\le 0.1\left(i-1\right)\\ i=11 \mathrm{when} 0.9<x<1.0\\ i=12 \mathrm{when} x=1.0\end{array}$$

(4)

where $M$ and $N$ denote the number of particles and particle sections, respectively.

The difference between ${\Lambda }_A^{2D\left(i\right)}$ and ${\Lambda }_A^{3D\left(i\right)}$ is due to the stereological bias (${\Lambda }_A^{dif\left(i\right)}$), and is calculated as follows:

$${\Lambda }_A^{dif\left(i\right)}={\Lambda }_A^{2D\left(i\right)}-{\Lambda }_A^{3D\left(i\right)}.$$

(5)

For the sake of convenience, the general term ${\Lambda }_A^{dif\left(i\right)}$ is expressed as ${\Lambda }_A^{dif}$, omitting the index $\left(i\right)$, and similarly for ${\Lambda }_A^{2D}$ and ${\Lambda }_A^{3D}$.

2.2. Modeling of Three-Dimensional Phase Structure

Boolean and Voronoi models were employed to create 3D binary structures. Prior to creating these structures, the rectangular prism (30 × 12 × 30) in which the particles were packed was digitized into (750 × 300 × 750) voxels. The voxels constituting the individual particles were tagged with the particle numbers, and those in the void were tagged as zero. These voxel sets were termed particle voxels for the sake of convenience.

2.2.1. Boolean Model

The following procedure was followed, with steps 2–5 corresponding to Figure 2b, and step 6 corresponding to Figure 2c.

• The rectangular prism was digitized into (750 × 300 × 750) voxels with a resolution of 0.04. All the voxels were tagged with the phase type, initially phase B.
• The following two parameters were given: the volumetric fraction of phase A ($F_v^0$) and the phase A element diameter ($d_A$).
• A spherical phase A element with diameter $d_A$ was randomly located in the rectangular prism, allowing overlap with other elements, and the phase of the voxels in the element was changed from B to A.
• The phase A volume fraction of the rectangular prism was calculated by counting the number of phase A voxels.
• Steps 3–4 were repeated until the phase A volume fraction of the rectangular prism exceeded $F_v^0$. Then, a binary Boolean 3D structure with the given parameters $F_v^0$ and $d_A$ was created. This set of voxels was termed the phase voxels, for the sake of convenience.
• The particle voxels and phase voxels were merged, and binary particles with a Boolean structure were created.

Figure 3a,b show a comparison of the partial sections following steps 5 and 6, when $F_v^0=0.5 \mathrm{and} d_A=1.0$; and Figure 3c–f show the partial sections following step 5, when (c) $F_v^0=0.2, d_A=1.0$, (d) $F_v^0=0.8, d_A=1.0$, (e) $F_v^0=0.5, d_A=3.6$, and (f) $F_v^0=0.5, d_A=0.4$, respectively.

2.2.2. Voronoi Model

The following procedure was followed, with steps 3–7 corresponding to Figure 2b, and step 8 corresponding to Figure 2c.

• The rectangular prism was digitized into (750 × 300 × 750) voxels, with each voxel tagged with two types of information: seed number (initially zero) and phase type (initially B).
• It was assumed that the two parameters $F_v^0$ and $d_A$ were similar to those of the Boolean model.
• $N_A^0$ seeds were randomly located in the rectangular prism, with each seed identified by a unique sequence number. $N_A^0$ was determined by $N_A^0=V_c/V_s$, where $V_c$ denotes the volume of a rectangular prism and $V_s\left(=\pi d_A^3/6\right)$ the volume of a sphere with diameter $d_A$.
• The seed numbers of the voxels were changed to the sequence numbers of the relatively most proximate seeds.
• A seed number was randomly selected, and the voxels with the selected seed number were assigned as phase A.
• The phase A volume fraction of the rectangular prism was calculated by counting the number of phase A voxels.
• Steps 5–6 were repeated until the phase A volume fraction of the rectangular prism exceeded $F_v^0$. Then the binary Voronoi 3D structure was created. This set of voxels was termed the phase voxels.
• The particle voxels and the phase voxels were merged, and binary particles with a Voronoi structure were created.

Figure 4a,b show a comparison of the partial sections after steps 7 and 8, when $F_v^0=0.5, d_A=1.0$; and Figure 3c–f show the partial sections following step 7, when (c) $F_v^0=0.2, d_A=1.0$, (d) $F_v^0=0.8, d_A=1.0$, (e) $F_v^0=0.5, d_A=3.6$, and (f) $F_v^0=0.5, d_A=0.4$, respectively.

2.3. Texture Analysis

As mentioned earlier, Kaneko’s fractal dimension technique (FD) [25], the spatial gray level dependence method (SGLDM) [26,27], and the gray level difference method (GLDM) [27] were employed for texture analysis, as follows.

2.3.1. Fractal Dimension Technique (FD)

Kaneko’s fractal dimension ($\delta$) [25] was employed for particle sectional analysis, based on the procedure described in a previous study [22], which is summarized in Figure 5. The procedure included the following steps.

(i)

Squares of size $d_{max}$ (maximum particle diameter) were superimposed on the particle sections. The squares were subdivided equally into $n^2$ squares of size $r=\left(d_{max}/n\right)$ (Figure 5a).

(ii)

Based on the small squares encompassed by the particle section, 3D structures were posited, with a width and length of $r$, and a height proportional to the image intensity ($g$). The overall surface area of the imaginary 3D structure was estimated by the summation of triangles ABD and BCD (Figure 5b,c). The summation of all such 3D structures in a given particle was defined as $A\left(r\right)$.

(iii)

A total of 50 $A\left(r\right)$ summations, with $n$ ranging from 1 to 50, were plotted with respect to $r$ on a double logarithmic chart. Additionally, $\delta$ was obtained from the least-square fitting line by applying the following equation:

$$\log A\left(r\right)=\left(2-\delta \right)\log r+C,$$

(6)

where $C$ denotes a constant.

It should be noted that FD is a method for assessing textural complexity, not particle or phase shape. In addition, FD has mathematical ambiguity in terms of the dimensional difference between $r$ and $g$; and has limitation in the range of $r$ in which Equation (6) shows linear relationship in the double logarithmic chart. However, it is practically valuable for texture analysis within a certain range of $r$.

2.3.2. Spatial Gray Level Dependence Method (SGLDM)

The co-occurrence matrix $P\left(i,j:\Delta ,\theta \right)$ shows the frequency of gray-level pixel pairs ($i, j$) with a relative distance of $\Delta$ in an angular direction of $\theta$ [26,27]. Therefore, $P\left(i,j:\Delta ,\theta \right)$ is a $g\times g$ matrix where $g$ denotes the number of gray levels. In this study, $g=2$, $\Delta =1$, and $\theta =0°,45°,90°,$ and $135°$. These parameters were determined as follows: $g$ corresponds to the number of phases here, so $g=2$; various values of $\Delta$ were employed in preliminary testing, and it was shown to be insensitive for stereological bias analysis, so $\Delta =1$ was selected; $\Delta$ of adjacent pixels is $\sqrt{2}$ for $\theta =45°$ and $135°$, but it is approximated by one; and four angles of $\theta$ were selected, to assess the omnidirectionality of the texture. Additionally, $P\left(i,j;\Delta ,\theta \right)$ was normalized as $p\left(i, j\right)$, as follows:

$$P\left(i,j\right)={\displaystyle \sum }_{\theta =0, 45, 90,135}P\left(i,j;\Delta ,\theta \right),$$

(7)

$$p\left(i,j\right)= \frac{P\left(i,j\right)}{{{\displaystyle \sum }}_{i=1}^g{{\displaystyle \sum }}_{j=1}^{g}P\left(i,j\right)}.$$

(8)

Figure 6 illustrates the calculation points for $P\left(i,j;\Delta ,\theta \right)$ in this study. As shown in Figure 6a, the particle sections overlapped with squares of size $d_{max}$ and subdivided into $n^2$ small squares, which were generated in step (i) of the fractal dimension procedure in Section 3.1. The points encompassed by the particle section were used for the $P\left(i,j;\Delta ,\theta \right)$ calculation (Figure 6b). Additionally, $n$ was set as 50 in accordance with the fractal dimension technique.

Haralick et al. [26] proposed 14 types of characteristics based on $p\left(i,j\right)$. In the present study, the 14 characteristics were examined in preliminary testing, and 5 candidates applicable to the proposed stereological correction method were selected, based on them having shown a clear correlation with stereological bias in the testing. Equations (9)–(11), (14) and (15) (below) defined the five characteristics as $f_{ASM}, f_{IDM}, f_{SV}, f_{EN},$ and $f_{DV}$, respectively corresponding to $f_1, f_5, f_7, f_9,$ and $f_{10}$ in [26].

Angular Second Moment (ASM)

$$f_{ASM}={\displaystyle \sum }_i{\displaystyle \sum }_j{\left\{p\left(i,j\right)\right\}}^2$$

(9)

Inverse Difference Moment (IDM)

$$f_{IDM}={\displaystyle \sum }_i{\displaystyle \sum }_j\frac{p\left(i,j\right)}{1+{\left(i-j\right)}^2}$$

(10)

Sum Variance (SV)

$$f_{SV}={{\displaystyle \sum }}_{i=2}^{2g}{\left(i-f_{SE}\right)}^2{p}_{x+y}^{}\left(i\right),$$

(11)

where

$$f_{SE}=-{\displaystyle \sum }_{i=2}^{2g}{p}_{x+y}^{}\left(i\right)\ln \left\{p_{x+y}^{}\left(i\right)\right\}$$

(12)

$$p_{x+y}^{}\left(k\right)= {\displaystyle \sum }_i{\displaystyle \sum }_j{p}^{}\left(i,j\right), k=i+j$$

(13)

Entropy (EN)

$$f_{EN}=-{\displaystyle \sum }_i{\displaystyle \sum }_{j}p\left(i,j\right)\ln \left\{p\left(i,j\right)\right\}$$

(14)

Difference Variance (DV)

$$f_{DV}=\frac{{{\displaystyle \sum }}_{k=0}^{g-1}{\left\{p_{\left(x-y\right)}\left(k\right)-{{\displaystyle \sum }}^{\text{}}\frac{p_{x-y}\left(k\right)}{g}\right\}}^2}{g-1},$$

(15)

where

$$p_{x-y}^{}\left(k\right)= {\displaystyle \sum }_i{\displaystyle \sum }_j{p}^{}\left(i,j\right), k=\left|i-j\right|.$$

(16)

2.3.3. Gray-Level Difference Method (GLDM)

A vector $Q\left(i\right)$ shows the frequency of gray-level difference between pixel pairs with a relative distance of $\Delta$ in an angular direction of $\theta$ [27]. In this study, $Q\left(i\right)$ is a $g$-level vector defined by parameters $g=2$, $\Delta =1$, and $\theta =0°,45°,90°,$ and $135°$, for the same reasons as in Section 2.3.2 (‘Spatial gray level dependence method (SGLDM)’). Additionally, $Q\left(i\right)$ was normalized as $q\left(i\right)$, as follows:

$$Q\left(i\right)={\displaystyle \sum }_{\theta =0, 45, 90,135}Q\left(i;\Delta ,\theta \right),$$

(17)

$$q\left(i\right)= \frac{Q\left(i\right)}{{{\displaystyle \sum }}_{i=1}^{g}Q\left(i\right)}.$$

(18)

The calculation points were determined similarly to the SGLDM (Figure 6).

Four $q\left(i\right)$ characteristics (Constant, Angular Second Moment, Entropy, and Mean) were examined. Through preliminary investigation, one characteristic of the constant ($f_{CN}$) was selected (as it alone showed a clear correlation with stereological bias), as follows:

Contrast (CN)

$$f_{CN}^{}={\displaystyle \sum }_{i=0}^{g-1}{i}^{2}q\left(i\right).$$

(19)

It should be noted that $f$ is used as a generic term for the aforementioned texture characteristic indices ($\delta ,f_{ASM}, f_{IDM}, f_{SV}, f_{EN},f_{DV}$, and $f_{CN}$), for the sake of convenience.

2.4. Numerical Simulation

2.4.1. Simulation for Development of the Stereological Correction Model

Based on the procedure in Section 2.1.1, the following spherical particles were created, with a total of 380 Boolean structure patterns: a combination of 19 types of $F_v^0$ ranging from 0.05 to 0.95 with a step size of 0.05, and 20 types of $d_A$ ranging from 0.2 to 4.0 with a step size of 0.2. Mineral liberation assessment, as in Section 2.1.2, was conducted for all 380 cases, and ${\Lambda }_A^{dif}$ and $F_a$ were obtained. The seven $f$ indices were calculated by texture analysis, as in Section 2.3. Finally, isograms of ${\Lambda }_A^{dif}$ were depicted with respect to $f$ and $F_a$, using the thin-plate spline method [45], based on the 380 calculated plots. Twelve types of isograms (${\Lambda }_A^{dif\left(i\right)}$, where $i=1-12$) were created for each $f$.

Figure 7 shows sample ${\Lambda }_A^{dif}$ isograms, with Figure 7a–d showing a comparison of ${\Lambda }_A^{dif\left(1\right)}$, ${\Lambda }_A^{dif\left(5\right)}$, ${\Lambda }_A^{dif\left(8\right)}$, and ${\Lambda }_A^{dif\left(12\right)}$ of $f_{EN}$; and Figure 7e–j respectively showing ${\Lambda }_A^{dif\left(1\right)}$ isograms of $\delta$, $f_{ASM}, f_{IDM}, f_{SV},$ $f_{DV}$, and $f_{CN}$, respectively. These ${\Lambda }_A^{dif}$ isograms, which are the key element in the proposed stereological correction method, are conceptually illustrated in Figure 1(iv).

2.4.2. Simulation for Model Validation

These important ${\Lambda }_A^{dif}$ isograms (illustrated in Figure 7) were based on the combination of spherical particles and Boolean internal structure, and thus the base model was termed the ‘Spherical–Boolean’ model. To investigate the viability of this model for the assessment of other particle systems, the following three systems were created.

• Irregularly shaped particles with a Voronoi internal structure, termed the ‘Irregular–Voronoi’ system.
• Irregularly shaped particles with a Boolean internal structure, termed the ‘Irregular–Boolean’ system.
• Spherical particles with a Voronoi internal structure, termed the ‘Spherical–Voronoi’ system.

The first tested the viability of the Spherical–Boolean model for assessment of particle systems with unfamiliar particle shape and internal structure, while the second and third focused on the relative importance of the specific parameters of particle shape and internal structure, respectively.

Each irregularly shaped particle was modeled using a geodesic grid [32]. Briefly, an ellipsoid was fitted with a triangular mesh, each node of which was randomly relocated in the radical direction to form surface roughness. The resulting particles had an aspect ratio (the quotient of the longest axis length and the shortest orthogonal axis length) of 2.0, and a surface roughness (the quotient of the surface area of the volume-equivalent ellipsoid and the surface area of the particle) of 0.95.

Nine variations of $d_A$ and $F_v^0$ were specified for each particle system.

Table 1. Model validation cases.

Case	Particle Shape	Internal Structure	${\mathit{d}}_{\mathit{A}}$	${\mathit{F}}_{\mathit{v}}^0$
I-A	Irregular	Voronoi	0.4	0.2
I-B	Irregular	Voronoi	0.4	0.5
I-C	Irregular	Voronoi	0.4	0.8
I-D	Irregular	Voronoi	1.0	0.2
I-E	Irregular	Voronoi	1.0	0.5
I-F	Irregular	Voronoi	1.0	0.8
I-G	Irregular	Voronoi	3.0	0.2
I-H	Irregular	Voronoi	3.0	0.5
I-I	Irregular	Voronoi	3.0	0.8
II-A–II-I	Irregular	Boolean	0.4, 1.0, and 3.0	0.2, 0.5, and 0.8
III-A–III-I	Spherical	Voronoi	0.4, 1.0, and 3.0	0.2, 0.5, and 0.8
Proposed correction model	Spherical	Boolean	0.2–4.0, in increments of 0.2	0.05–0.95, in increments of 0.05

tabulates the validation cases. The first letter of the case name (e.g., I) denotes the particle system, and the second (e.g., A) denotes the specification of $d_A$ and $F_v^0$. The second letter (A through I) corresponds to plots A through I in Figure 7. It should be noted that cases II and III each included nine variations, similar to those of case I; however, they are combined in a single line in Table 1 for the sake of simplicity. The specifications in the Spherical-Boolean correction model is shown, for comparison, in the bottom row of Table 1.

The validation of the stereological correction model was performed as follows.

• $F_a$, ${\Lambda }_A^{2D}$, and $f$ of all the particle systems were calculated based on the procedures detailed in Section 2.1 and Section 2.3.
• ${\Lambda }_A^{2D}$ was corrected as ${\Lambda }_A^{3D}$ using the ${\Lambda }_A^{dif}$ isograms (e.g., Figure 7) using the proposed correction method shown in Figure 1. The estimated ${\Lambda }_A^{3D}$ was termed ${\Lambda }_A^{3D}\prime$ for convenience.
• Additionally, the true value of ${\Lambda }_A^{3D}$ was calculated, using the relevant procedure in Section 1, as a correct answer.
• The gaps between ${\Lambda }_A^{2D}$ and ${\Lambda }_A^{3D}$, and ${\Lambda }_A^{3D}\prime$ and ${\Lambda }_A^{3D}$ were assessed in terms of the areal difference between the respective composition distribution curves ($E_1$). Generally, higher $E_1$ values imply a larger magnitude of stereological bias. $E_1$ is detailed in a previous study [21].
• $E_1$ was used to assess the stereological correction efficiency of the seven types of $f$s, and the improvement rate ($I_r$) was evaluated by the following equation:

  $$I_r=100\frac{\left(E_1^o-E_1^c\right)}{E_1^o},$$

  (20)

  where $E_1^o$ denotes $E_1$ without the correction calculated from ${\Lambda }_A^{2D}$ and ${\Lambda }_A^{3D}$, and $E_1^c$ denotes $E_1$ after the correction calculated from ${\Lambda }_A^{3D}\prime$ and ${\Lambda }_A^{3D}$.

3. Results

3.1. Applicability of the Correction Model to Particle Systems with Unfamiliar Particle Shape and Internal Structure

The applicability of the correction model to particle systems with unfamiliar particle shape and internal structure was investigated using Case I (irregularly-shaped particles with a Voronoi internal structure). Figure 8a compares the 2D composition distribution (${\Lambda }_A^{2D}$), true 3D composition distribution (${\Lambda }_A^{3D}$), and seven types of 3D composition distribution (${\Lambda }_A^{3D\prime }$) based on the seven types of $f$s, for Case I-E, along with the corresponding partial sectional image. Note that the particle section image in Figure 8a represents only part of the sample section; as mentioned before, a total of 4316 particle sections were examined. Typically, the difference between the 2D and 3D curves due to stereological bias was considerably reduced by all the correction methods. Since the corrected curves were indistinguishable, a partial close-up of the chart is shown in Figure 8b.

Table 2. Error index ($E_1$) before and after correction using seven types of texture analysis indices (FD, ASM, IDM, SV, EN, DV, and CN), for Case I (irregularly-shaped particles with Voronoi internal structure), along with the corresponding average improvement rate ($I_r$).

Case	${\mathit{d}}_{\mathit{A}}$	${\mathit{F}}_{\mathit{v}}^0$	${\mathit{E}}_1$
			Before Correction	FD ($\mathit{\delta }$)	ASM (${\mathit{f}}_{\mathit{A}\mathit{S}\mathit{M}}$)	IDM (${\mathit{f}}_{\mathit{I}\mathit{D}\mathit{M}}$)	SV (${\mathit{f}}_{\mathit{S}\mathit{V}}$)	EN (${\mathit{f}}_{\mathit{E}\mathit{N}}$)	DV (${\mathit{f}}_{\mathit{D}\mathit{V}}$)	CN (${\mathit{f}}_{\mathit{C}\mathit{N}}$)
I-A	0.4	0.2	0.04660	0.01557	0.00809	0.01663	0.01366	0.00952	0.01631	0.01656
I-B	0.4	0.5	0.06392	0.01582	0.01054	0.01809	0.01897	0.01232	0.01789	0.01801
I-C	0.4	0.8	0.04738	0.01271	0.01178	0.01265	0.01547	0.01176	0.01244	0.01262
I-D	1.0	0.2	0.04472	0.00890	0.00779	0.01646	0.00770	0.00721	0.01580	0.01623
I-E	1.0	0.5	0.06035	0.01167	0.01383	0.01250	0.01621	0.01310	0.01202	0.01238
I-F	1.0	0.8	0.04381	0.01223	0.01177	0.01056	0.01771	0.01170	0.01056	0.01056
I-G	3.0	0.2	0.02004	0.00242	0.00340	0.00522	0.00315	0.00293	0.00472	0.00493
I-H	3.0	0.5	0.02841	0.00422	0.00480	0.00441	0.00578	0.00444	0.00395	0.00432
I-I	3.0	0.8	0.02031	0.00626	0.00514	0.00273	0.00571	0.00505	0.00279	0.00290
Average $I_r$ (%)			-	76.7	79.4	74.7	73.3	79.5	75.6	75.0

FD—fractal dimension; ASM—angular second moment; IDM—inverse difference moment; SV—sum variance; EN—entropy; DV—difference variance; CN—contrast; and $I_r$—Improvement rate (Equation (20)).

tabulates the $E_1$ for Case I, with the seven $f$ indices, together with the respective improvement rates ($I_r$). Generally, the stereological bias was sufficiently overcome by all the corrections. In terms of the applicability of the model to particle systems with unfamiliar particle shape and internal structure, ASM and EN exhibited the highest efficiency (approximately 80% $I_r$), while the others exhibited somewhat lower efficiency (approximately 73%–77% $I_r$).

3.2. Influence of Unfamiliar Particle Shape on the Viability of the Correction Model

The influence of unfamiliar particle shape on the viability of the correction model was investigated using Case II (irregularly-shaped particles with Boolean internal structure). Figure 9a compares ${\Lambda }_A^{2D}$, ${\Lambda }_A^{3D}$, and ${\Lambda }_A^{3D\prime }$, as estimated by the seven types of $f$, for Case II-E, along with the corresponding partial sectional image. The partial close-up of the chart is shown in Figure 9b.

Table 3. Error index ($E_1$) before and after correction using seven types of texture analysis indices (FD, ASM, IDM, SV, EN, DV, and CN), for Case II (irregularly-shaped particles with Boolean internal structure), along with the corresponding average improvement rate ($I_r$).

Case	${\mathit{d}}_{\mathit{A}}$	${\mathit{F}}_{\mathit{v}}^0$	${\mathit{E}}_1$
			Before Correction	FD ($\mathit{\delta }$)	ASM (${\mathit{f}}_{\mathit{A}\mathit{S}\mathit{M}}$)	IDM (${\mathit{f}}_{\mathit{I}\mathit{D}\mathit{M}}$)	SV (${\mathit{f}}_{\mathit{S}\mathit{V}}$)	EN (${\mathit{f}}_{\mathit{E}\mathit{N}}$)	DV (${\mathit{f}}_{\mathit{D}\mathit{V}}$)	CN (${\mathit{f}}_{\mathit{C}\mathit{N}}$)
II-A	0.4	0.2	0.04344	0.00996	0.00615	0.01196	0.01196	0.00680	0.01162	0.01190
II-B	0.4	0.5	0.05409	0.01617	0.00805	0.01665	0.02519	0.01060	0.01633	0.01662
II-C	0.4	0.8	0.03467	0.01245	0.00820	0.01159	0.00766	0.00834	0.01137	0.01156
II-D	1.0	0.2	0.04167	0.00555	0.00364	0.01315	0.00422	0.00299	0.01262	0.01288
II-E	1.0	0.5	0.05268	0.01130	0.00386	0.01704	0.00712	0.00446	0.01634	0.01686
II-F	1.0	0.8	0.03730	0.00828	0.00366	0.01309	0.00760	0.00448	0.01243	0.01297
II-G	3.0	0.2	0.01616	0.00126	0.00152	0.00506	0.00088	0.00096	0.00484	0.00475
II-H	3.0	0.5	0.02887	0.00188	0.00308	0.00357	0.00313	0.00274	0.00311	0.00338
II-I	3.0	0.8	0.02327	0.00110	0.00177	0.00488	0.00189	0.00144	0.00427	0.00471
Average $I_r$ (%)			-	81.7	88.2	71.6	81.7	87.9	73.0	72.2

FD—fractal dimension; ASM—angular second moment; IDM—inverse difference moment; SV—sum variance; EN—entropy; DV—difference variance; CN—contrast; and $I_r$—Improvement rate (Equation (20)).

tabulates the $E_1$ for Case II, with the seven $f$ indices, along with $I_r$. Generally, the stereological bias was considerably overcome by all the corrections. In terms of the proposed model’s applicability to unfamiliar particle shapes, ASM and EN exhibited the highest correction efficiency in this respect (approximately 88% $I_r$), FD and SV exhibited the second highest efficiency (approximately 82% $I_r$), and IDM, DV, and CN exhibited low efficiency (approximately 71%–73% $I_r$).

The stereological bias in the case of the Spherical-Spherical model was perfectly corrected, and the corresponding $I_r$ was 100%. Thus, the influence of unfamiliar particle shape was limited to 12% (=100% − 88%) of $I_r$, for ASM and EN, under the specified simulation conditions.

3.3. Influence of Unfamiliar Internal Structure on the Viability of the Correction Model

The influence of unfamiliar internal structure on the viability of the correction model was investigated using Case III (spherical particles with Voronoi internal structure). Figure 10 and

Table 4. Error index ($E_1$) before and after correction using seven types of texture analysis indices (FD, ASM, IDM, SV, EN, DV, and CN), for Case III (spherical particles with Voronoi internal structure), along with the corresponding average improvement rate ($I_r$).

Case	${\mathit{d}}_{\mathit{A}}$	${\mathit{F}}_{\mathit{v}}^0$	${\mathit{E}}_1$
			Before Correction	FD ($\mathit{\delta }$)	ASM (${\mathit{f}}_{\mathit{A}\mathit{S}\mathit{M}}$)	IDM (${\mathit{f}}_{\mathit{I}\mathit{D}\mathit{M}}$)	SV (${\mathit{f}}_{\mathit{S}\mathit{V}}$)	EN (${\mathit{f}}_{\mathit{E}\mathit{N}}$)	DV (${\mathit{f}}_{\mathit{D}\mathit{V}}$)	CN (${\mathit{f}}_{\mathit{C}\mathit{N}}$)
III-A	0.4	0.2	0.03688	0.00706	0.00437	0.00449	0.00562	0.00431	0.00452	0.00448
III-B	0.4	0.5	0.05149	0.00670	0.00615	0.00605	0.02986	0.00610	0.00603	0.00605
III-C	0.4	0.8	0.03695	0.01094	0.00687	0.00817	0.02677	0.00718	0.00821	0.00817
III-D	1.0	0.2	0.04701	0.00418	0.00356	0.00429	0.00330	0.00359	0.00431	0.00424
III-E	1.0	0.5	0.06185	0.01146	0.00942	0.00990	0.01139	0.00954	0.00990	0.00990
III-F	1.0	0.8	0.04766	0.01734	0.01418	0.01621	0.01771	0.01461	0.01619	0.01622
III-G	3.0	0.2	0.02420	0.00260	0.00290	0.00302	0.00319	0.00276	0.00276	0.00303
III-H	3.0	0.5	0.03466	0.00573	0.00554	0.00567	0.00548	0.00555	0.00548	0.00569
III-I	3.0	0.8	0.02279	0.00466	0.00269	0.00455	0.00303	0.00292	0.00444	0.00459
Average $I_r$ (%)			-	80.7	85.0	82.9	72.2	84.8	83.1	82.9

FD—fractal dimension; ASM—angular second moment; IDM—inverse difference moment; SV—sum variance; EN—entropy; DV—difference variance; CN—contrast; and $I_r$—Improvement rate (Equation (20)).

correspond to Figure 9 and Table 3, respectively. It was noted that the correction enabled the cumulative rate to fall below zero when $x=0$ (Figure 10b); and similarly, it was possible for the correction curves to exceed one when $x=1$. These unrealistic values could occur in the proposed correction algorithm; however, in practice, the unrealistic values could be modified (as zero or one), and their effect was considered to be negligible.

Generally, the stereological bias was considerably overcome by all the corrections (Table 4). In terms of the model’s applicability to unfamiliar particle shapes, ASM and EN again exhibited the highest efficiency (approximately 85% $I_r$), IDM, FD, DV, and CN exhibited the second highest efficiency (approximately 81%–83% $I_r$), and SV exhibited the lowest efficiency (approximately 72% $I_r$). Thus, the influence of unfamiliar internal structure was limited to 15% of $I_r$, for ASM and EN, under the specified simulation conditions.

Comparison of Cases II and III indicated that unfamiliar particle shape and unfamiliar internal structure reduced the viability of the correction model to roughly the same extent, but the reduction was slightly greater in the latter case. This is intuitively understandable because particle shape determined the range of the texture analysis (e.g., number of black dots in Figure 6b), while internal structure (e.g., section pattern in Figure 6b) affected the analysis in a more direct manner.

4. Discussion

It remains to discuss why ASM and EN exhibited the highest correction efficiency. Figure 11 conceptually depicts the $p\left(i,j\right)$ matrix of the SGLDM. Here, the diagonal elements ($p\left(0, 0\right)$ and $p\left(1, 1\right)$, shaded in Figure 11) denote the frequency of pixel pairs in the same phase, while the non-diagonal elements ($p\left(1, 0\right)$ and $p\left(0, 1\right)$, unshaded) denote the frequency of pixel pairs in different phases. Therefore, a texture with relatively larger diagonal than non-diagonal elements indicates a higher degree of phase aggregation. It is intuitively understandable that the degree of phase aggregation has a significant influence on the stereological bias in liberation assessment, since less stereological bias occurs, for example, in completely randomly distributed phase systems such as snow noise.

In the present study, the respective magnitudes of all the $p\left(i,j\right)$ elements were directly calculated for ASM and EN, based on Equations (9) and (14). Hence, ASM and EN efficiently reflected the gap between the diagonal and non-diagonal $p\left(i,j\right)$ elements (i.e., the respective degrees of phase aggregation), and this may have resulted in ASM and EN exhibiting the highest stereological correction efficiency.

An extension of the abovementioned consideration could involve establishing texture characteristic indices ($f$) specifically designed for stereological correction, by assessing the diagonal elements of $p\left(i,j\right)$. Additionally, $f$ can be customized to a certain targeted ore by adjusting $\Delta$ and $\theta$ (Equation (7)), where $\Delta$ and $\theta$ are considered to correlate with the size and anisotropy of the phase grains, respectively.

Boolean and Voronoi models are used for modeling two types of textures, but their validity for modeling real ore texture was not assessed in this study. Some Voronoi variations are specifically designed to model the internal structure of real ore (e.g., the weighted Voronoi model [40]). It may be possible to improve the proposed stereological correction method by using these advanced structural models.