3D XMT Data Preprocessing and Segmentation

In this section, we give a short overview of the image preprocessing steps and the particle-wise segmentation of the 3D XMT data considered in the present paper; for more details see Furat et al. (Reference Furat, Leißner, Ditscherlein, Sedivy, Weber, Bachmann, Gutzmer, Peuker and Schmidt2018). The first step consists of the reduction of noise in the XMT data with a non-local means denoising algorithm, see Buades et al. (Reference Buades, Coll and Morel2005). In contrast to the commonly used Gaussian filter, this nonlinear image filter has the advantage of smoothing homogeneous regions while preserving edges. In order to separate the particles from the background we use the local adaptive Sauvola thresholding algorithm, which determines binarization thresholds for each voxel based on its local neighborhood, see Shafait et al. (Reference Shafait, Keysers and Breuel2008). Such local thresholding techniques can produce better results than global thresholds, due to globally inconsistent grayscale values in XMT images. One crucial and nontrivial task in image processing is the particle-wise segmentation, which allows the extraction of single particles from the image data for quantitative analysis. In a first step, we use a marker-based watershed algorithm, see Spettl et al. (Reference Spettl, Wimmer, Werz, Heinze, Odenbach, Krill and Schmidt2015), to obtain an initial segmentation. This results in some particles being wrongly separated into multiple fragments. To remove such oversegmentations, we train a neural network (Hastie et al., Reference Hastie, Tibshirani and Friedman2009) to decide whether adjacent segments should be merged or not. To make this decision the neural network is supplied with information regarding the local morphology and grayscale values around two adjacent fragments. By applying the neural network as a post-processing step on the initial segmentation, we receive our final segmentation, see Figure 1b.

Fig. 1. a: Volumetric XMT data (gray) with two registered 2D SEM-EDS planes (blue and red). b: Cut-out of the particle-wise segmentation of the XMT data.

2D SEM-EDS Data and Registration

In addition to the 3D XMT image data, SEM-EDS data were considered by Furat et al. (Reference Furat, Leißner, Ditscherlein, Sedivy, Weber, Bachmann, Gutzmer, Peuker and Schmidt2018). The latter were obtained with a mineral liberation analyzer from the same sample at two spatially different planar sections. An illustration of the SEM-EDS data can be seen in Figure 2a. It provides 2D information about the morphology of the particles in planar sections, but in contrast to the XMT data, see Figure 2b, it also provides information about the mineralogical composition of particles. For example, in the false color image of Figure 2a the blue phase indicates zinnwaldite.

Fig. 2. a: 2D SEM-EDS image, the blue phase indicates zinnwaldite and pink indicates quartz. b: The corresponding section registered in the 3D XMT data.

Thus, by localizing the 2D SEM-EDS data in the segmented 3D image (registration), see Figure 1a, we have knowledge about the mineralogical composition of each 3D particle that intersects with the planar SEM-EDS data. Due to the particle-wise segmentation described in the “3D XMT Data Preprocessing and Segmentation” section we also know the 3D morphology of each of these particles. Furthermore, we can link the particle-wise segmentation to the grayscale values of the 3D XMT image, since the grayscale values provide valuable information about the composition of particles, as shown in Furat et al. (Reference Furat, Leißner, Ditscherlein, Sedivy, Weber, Bachmann, Gutzmer, Peuker and Schmidt2018). Therefore, both the grayscale values, and the morphology of particles in the XMT image that hit the SEM-EDS plane, allow us to adjust a classifier that can predict the mineralogical composition of an arbitrary particle solely based on information gained from the XMT image.

To be precise, let I : W ⊂ ℤ³ → [0, 1] be the grayscale XMT image, where W is a cuboid observation window. The 2D SEM-EDS data can be regarded as a map $L:H\cap W\subset { {\opf Z}}^3\to \lcub {0,1,2, \ldots} \rcub $, where H is a set of voxels representing a plane and the values of L indicate different minerals or the background, e.g., we put L (x) = 0 if background was observed at $x\in H\cap W$, or L (x) = 1 if zinnwaldite was observed at $x\in H\cap W$. Furthermore, let P ₁, …, P _n ⊂ W ⊂ ℤ³ be the sets of voxels corresponding to the particles that are obtained by the particle-wise segmentation of the 3D XMT image and that hit the plane H, i.e., $P_k\cap H = \lcub {x_1^{\lpar k \rpar }, \ldots, x_{\ell_k}^{\lpar k \rpar }} \rcub \ne \emptyset $ for each k = 1, …, n. Each particle P _k gets the label

(1)$$L(P_k) = {\rm mode}\,\left( {L\left( {x_1^{( k ) }} \right) , \ldots, L\left( {x_{\ell_k}^{\lpar k \rpar }} \right) } \right) ,$$

where the mode of the sample of labels $\lpar {L\lpar {x_1^{\lpar k \rpar }} \rpar , \ldots, L\lpar {x_{\ell_k}^{\lpar k \rpar }} \rpar } \rpar $ is the label that appears most often. Thus, we assign to each particle P _k the mineral that is observed most frequently in the intersecting plane H. This means, for example, that a composite particle P _k which consists of zinnwaldite, quartz, and other minerals will be referred to as a zinnwaldite-composite particle if its main component in the intersection $P_k\cap H$ is zinnwaldite. Since we make this labeling only based on the mineralogical observation in the plane H we assume stationarity of the mineralogical composition of the particles P _k, i.e., that the composition of a particle P _k outside of the plane H is adequately represented by $P_k\cap H$. In our data, the set of particles observed in the plane H can be decomposed in three disjoint sets:

(2)$$\lcub {P_1, \ldots, P_n} \rcub = Z\cup Q\cup O$$

where Z = {P _k : L(P _k) = 1} are the zinnwaldite-composite particles and Q = {P _k : L(P _k) = 2} are the quartz-composite particles. The set O = {P _k : L(P _k) ≥ 3} contains the remaining particles that were observed in H. In the considered plane, H contains 861 particles, from which 342 belong to the set of zinnwaldite-composite particles Z and 462 to the set of quartz-composite particles Q. Since the set O contains only 57 particles, we disregard these observations from now on. Note that due to the labeling performed according to equation (1), a particle P ∈ Z does not necessarily consist solely of zinnwaldite. Figure 3a shows the distribution of the zinnwaldite volume fraction for particles labeled as zinnwaldite-composite particles in the cross-section $P\cap H$. Analogously, the quartz particles P ∈ Q are consolidated, see Figure 3b. The goal of the present paper is to develop a method which allows us to predict the main mineralogical component L(P) of particles when the additional planar SEM-EDS data are not available.

Fig. 3. a: Histogram of the zinnwaldite volume fractions for particles for which the majority mineral phase is zinnwaldite. The mean zinnwaldite volume fraction is 0.925 with a variance of 0.018. b: Histogram of the quartz volume fractions for particles for which the majority mineral phase is quartz. The mean quartz volume fraction is 0.937 with a variance of 0.013.

Size Characteristics

Relevant size descriptors of a particle P′ ⊂ ℝ³, which we only observe on the grid as P ⊂ ℤ³, are the volume ν ₃(P′) and the surface area a(P′). It is clear that the volume can be estimated by counting voxels belonging to the set P and the surface area can be estimated by considering suitably defined voxel configurations, see Schladitz et al. (Reference Schladitz, Ohser, Nagel, Kuba, Nyúl and Palágyi2006). We denote the estimated particle volume and surface area by ν ₃(P) and a(P), respectively.

For our prediction models we use the following estimator of the volume equivalent radius

(3)$$r(P) = \root \scriptstyle{\vskip3pt \hskip-1.5pt 3} \of {\displaystyle{3 \over {4\pi}} \nu _3(P)}, $$

which is in a one-to-one relationship with the volume. The surface area a(P) will not be directly regarded as a particle characteristic, but it is required for the sphericity factor which is a shape characteristic considered in the next section.

Shape Characteristics

In the previous section, we only considered the volume equivalent radius r as a size characteristic, but this is not enough to reliably distinguish between zinnwaldite- and quartz-composite particles, since particles of similar sizes can have very different appearance. Therefore, we examine some further characteristics for describing the shape of particles. One of the shape characteristics we consider is the sphericity factor

(4)$$s(P) = \displaystyle{{\root \scriptstyle{\vskip1pt \hskip-1pt 3} \of {36\pi \nu _3^2 \,(P)}} \over {a(P)}},$$

which takes values in [0, 1], where s(P′) = 1 holds if P′ ⊂ ℝ³ is a sphere.

Another shape characteristic of interest is the convexity factor

(5)$$c(P) = \displaystyle{{\nu _3(P)} \over {\nu _3(q(P))}},$$

where q(P) ⊂ ℤ³ is the convex hull of P on the lattice ℤ³. Like the sphericity factor, the convexity factor takes values in [0, 1], where c(P) = 1 holds if and only if P is convex on ℤ³. There is a causality between the sphericity factor s and the convexity factor c since more spherically shaped particles also have higher convexity factors. Yet, a convexity factor of 1 does not necessarily imply a large sphericity factor, since, for example, ellipsoids have always a convexity factor of 1, yet their sphericity factor can be arbitrarily close to 0. Since we observe a large number of flat particles in our image data, we also consider characteristics that quantify elongation of particles. For each particle P ⊂ ℤ³ with the barycenter $\bar{x} = \,(\bar{x}_1,\bar{x}_3,\bar{x}_3)\in {\opf R}^3$, whose coordinates are given by

(6)$$\bar{x}_i = \displaystyle{1 \over {\nu _3(P)}}\mathop \sum \limits_{x\in P} x_i\quad {\rm for}\;i = 1,2,3,$$

where x = (x ₁, x ₂, x ₃), we consider the (positive-semidefinite) covariance matrix:

(7)$$C = \left( {\displaystyle{1 \over {\nu_3(P)}}\mathop \sum \limits_{x\in P} (x_i-{\bar{x}}_i)(x_j-{\bar{x}}_j)} \right)_{i,j = 1,2,3}.$$

The eigenvalues 0 ≤ λ ₁ ≤ λ ₂ ≤ λ ₃ of C, where we assume that λ ₃ >0, are closely related to the axis lengths of the best fitting ellipsoid ε(P) corresponding to the particle P, i.e., the axis lengths a ₁, a ₂, a ₃ of ε(P) are given by $a_i = \gamma \sqrt {\lambda _i} $ for i = 1, 2, 3, where γ is some scaling factor. The elongation factor e(P) of the particle P is then given by

(8)$$e(P) = \displaystyle{{a_1} \over {a_3}} = \sqrt {\displaystyle{{\lambda _1} \over {\lambda _3}}}. $$

Note that, like the previously considered characteristics, the elongation factor e(P) is normalized and takes values in [0, 1], where values close to 0 indicate elongated particles such as rod- or plate-shaped particles. By additionally analyzing the relationship between λ ₂ and λ ₃ it is possible to distinguish between rods and plates, but since we observed no rod-like particles in our data set, this was not necessary in the present paper.

Grayscale Characteristics

The XMT data do not only provide information about the 3D morphology of particles. In Furat et al. (Reference Furat, Leißner, Ditscherlein, Sedivy, Weber, Bachmann, Gutzmer, Peuker and Schmidt2018), we have seen that the grayscale values of XMT images contain substantial information about the mineralogical composition of particles. Recall that the grayscale values of a particle P = {x ₁, …, x _n} ⊂ W are given by y ₁ = I(x ₁), …, y _n = I(x _n) ∈ [0, 1], where I is the XMT image. Furthermore, we assume without loss of generality that the grayscale values are ordered, i.e., y ₁ ≤ y ₂ ≤ … ≤ y _n. This allows us to define the following grayscale characteristics.

A natural choice could be the mean

(9)$$\bar{y} = \displaystyle{1 \over n}\mathop \sum \limits_{k = 1}^n y_k.$$

But, since particles often have some imperfections, such as small regions of high density, see Figure 2, the mean is not well suited for representing the dominant grayscale value. Therefore, we used the more robust median

(10)$$m\,(P) = y_{0.5}$$

as grayscale characteristic, where the median is defined by the empirical quantiles of the sample

(11)$$y_p = \left\{ {\matrix{ {y_{[ {np} ] + 1},} \hfill & {\;{\rm if}\;np\in {\rm {\opf N}},} \hfill \cr {\displaystyle{{(y_{[ {np} ] } + y_{[{np} ] + 1})} \over 2},} \hfill & {\;{\rm if}\;np\notin {\rm {\opf N}},} \hfill \cr}} \right.$$

for p = 0.5. A natural choice for measuring the variability of grayscale values of a particle could be the sample standard deviation

(12)$$\sigma = \sqrt {\displaystyle{1 \over {n-1}}\mathop \sum \limits_{k = 1}^n \,{(y_k-\bar{y})}^2}. $$

Similarly to the mean, the standard deviation is not a robust particle descriptor. Therefore, we use the interquartile range

(13)$${\rm iqr}\,(P) = y_{0.75}-y_{0.25}$$

to characterize the variability of the grayscale values of the dominant grayscale region of the particle. Further characteristics that describe the texture of particles (Shapiro & Stockman, Reference Shapiro and Stockman2001) will be investigated in a forthcoming paper.

Modeling of Single Particle Characteristics

We now fit parametric probability distributions to the 1D particle characteristics described in the “Particle Characteristics” section. To be precise, we first fit parametric distributions to the volume equivalent radius, sphericity factor, convexity factor, elongation factor, median grayscale value, and the iqr for the 3D zinnwaldite-composite particles that intersect with the SEM-EDS plane. The chosen families of parametric distributions and the corresponding parameters that were determined to model the distributions of these characteristics are listed in Table 3. Note that we used gamma- and beta-distributions, whose probability density functions are given by

(17)$$f_{{\rm gamma}}(x) = \displaystyle{1 \over {\theta ^k\Gamma (k)}}x^{k-1}\exp \left( {-\displaystyle{x \over \theta}} \right){\rm \;} {\tf="PSX238E" 0}_{x \gt 0}$$

and

(18)$$f_{{\rm beta}}\,(x) = \displaystyle{{\Gamma (\alpha )\Gamma (\beta )} \over {\Gamma (\alpha + \beta )}}x^{\alpha -1}(1-x)^{\beta -1}{\tf="PSX238E" 0}_{x\in [{0,1}] },$$

respectively, where k, θ, α, β > 0 are some model parameters and the gamma function Γ:[0, ∞) → [0, ∞) is defined by the integral

(19)$${\rm \Gamma (}\alpha {\rm )} = \mathop \int \limits_0^\infty x^{\alpha -1}{\rm e}^{-x}{\rm d}x\quad {\rm for}\;{\rm each}\;\alpha \gt 0.$$

For modeling the volume equivalent radius we used the gamma distribution. Note that there are also other families of distributions which result in a good fit, like the log-normal distribution. The remaining characteristics were modeled with beta distributions, since these particle characteristics have only values in the interval [0, 1] which coincides with the support of the beta distribution. A visual comparison with the histograms of particle characteristics is given in Figure 6. In order to determine the model parameters (k, θ for the gamma and α, β for the beta distribution) for these 1D fits we used maximum-likelihood estimators (Held & Bové, Reference Held and Bové2014). Analogously, the 1D distributions of single particle characteristics were fitted for quartz-composite particles, see Figure 7 and Table 3.

Fig. 6. Histograms (blue) and fitted probability densities of zinnwaldite-composite particle characteristics for (a) volume equivalent radius, (b) sphericity factor, (c) convexity factor, (d) elongation factor, (e) median grayscale value, and (f) iqr, see also Table 3.

Fig. 7. Histograms (blue) and fitted probability densities of quartz-composite particle characteristics for (a) volume equivalent radius, (b) sphericity factor, (c) convexity factor, (d) elongation factor, (e) median grayscale value, and (f) iqr.

Table 3. Parameters of Fitted Distributions.

Modeling of Multidimensional Particle Characteristics

When modeling the joint distribution of independent random particle characteristics X,  Y with probability density functions f _X and f _Y, respectively, the joint distribution is immediately given by f _(X, Y) (x, y) = f _X(x)f _Y(y). However, this approach is not available for correlated characteristics, as it is the case for our data, where, for example, we obtained a correlation coefficient of 0.36 for the sphericity factor and median grayscale value of zinnwaldite-composite particles.

Thus, in order to describe the joint distribution of vectors of particle characteristics, whose marginal distributions are given by the 1D distributions we fitted in the “Modeling of Single Particle Characteristics” section, see Table 3, we consider a more general approach, using the so-called Archimedean copulas, see Nelsen (Reference Nelsen2006).

To make the paper more self-contained, we first explain some basic ideas of the copula approach, which will be used to construct multivariate probability distributions with non-normal marginals. Therefore, for some d ≥ 2, let f ₁, …, f _d : ℝ → [0, ∞) denote the 1D probability densities of particle characteristics for which we want to determine the d-dimensional joint probability density f : ℝ^d → [0, ∞) whose marginal densities are given by f ₁, …, f _d. In order to construct the multivariate density f we need to consider, for technical reasons, the 1D (cumulative) distribution functions F ₁, …, F _d : ℝ → [0, 1] which are given by

(20)$$F_i\,(x) = \mathop \int \limits_{-\infty} ^x f_i\,(y)\,{\rm d}y\quad {\rm for}\;i = 1, \ldots, d.$$

Note that F ₁(x) denotes the probability that the particle characteristic described by the density f ₁ does not exceed the value x ∈ ℝ and that the density can be obtained from the distribution function F _i by

(21)$$f_i(x) = \displaystyle{{\rm d} \over {{\rm d}x}}F_i(x)\quad {\rm for}\;i = 1, \ldots, d.$$

Analogously, the d-dimensional density f has a cumulative distribution function F : ℝ^d → [0, 1] which is given by

(22)$$\eqalign{&F(x) = \mathop \int \limits_{-\infty} ^{x_1} \ldots \mathop \int \limits_{-\infty} ^{x_d} f\,(y_1, \ldots, y_d)\,{\rm d}y_d \ldots {\rm d}y_1 \cr & {\rm for}\;x = (x_1, \ldots, x_d)\in {\rm {\opf R}}^d.$$

The joint density f is then obtained by the d-fold partial derivative

(23)$$f\,(x) = \displaystyle{{\partial ^d} \over {\partial x_1\cdots \partial x_d}}F(x)\quad {\rm for}\;x = (x_1, \ldots, x_d)\in {\rm {\opf R}}^d.$$

In order to show the connection between the multivariate distribution function f and the notion of copulas we begin with the definition of the latter. A d-dimensional copula is a multivariate cumulative distribution function K : ℝ^d → [0, 1], whose 1D marginal distributions are uniform distributions in the interval [0, 1]. For example, the marginal cumulative distribution function K ₁ : ℝ → [0, 1] of the first component is given by

(24)$$\eqalign{&K_1(x_1) & = \mathop {\lim} \limits_{x_2, \ldots, x_d\to \infty} K\lpar {x_1,x_2, \ldots, x_d} \rpar \cr & = \left\{ {\matrix{ {0,} \hfill & {\;{\rm if}\;x_1 \lt 0,} \hfill \cr {x_1,} \hfill & {\;{\rm if}\;x_1\in [ {0,1} ],} \hfill \cr {1,} \hfill & {\;{\rm if}\;x_1 \gt 1.} \hfill \cr}} \right.}$$

Copulas are of special interest because of Sklar's theorem (Nelsen, Reference Nelsen2006). It states that F is a d-dimensional cumulative distribution function with marginal distribution functions F ₁, …, F _d if and only if there is a copula function K such that

(25)$$F(x) = K(F_1(x_1), \ldots, F_d(x_d))\quad {\rm for}\;x = (x_1, \ldots, x_d)\in {\rm {\opf R}}^d.$$

This means that every multivariate distribution function F can be represented by its marginals F ₁, …, F _d and a copula function K. Thus, when modeling the joint distribution function F of random vectors whose marginal distributions are given by F ₁, …, F _d, it suffices to model the copula K. Note that there are numerous classes of copula functions, for which many of them do not have an analytical representation. In the present paper we limit the search for an appropriate copula function to the so-called Archimedean copulas.

Archimedean Copulas and Differential Variant of Sklar's Theorem

For that purpose, let φ : [0, 1] → [0, ∞] be a continuous, strictly decreasing and convex function with φ(1) = 0 and φ(0) = ∞. Note that φ is called an Archimedean generator. Then, it can be shown that the function K : [0, 1]^d → [0, 1] given by

(26)$$\eqalign{&K(u) = \varphi ^{-1}(\varphi (u_1) + \cdots + \varphi (u_d)) \cr & {\rm for}\;u = (u_1, \ldots, u_d)\in [0,1]^d$$

can be uniquely extended to a copula on ℝ^d, i.e., a function possessing the properties mentioned in the “Modeling of Multidimensional Particle Characteristics” section. It is called an Archimedean copula (Nelsen, Reference Nelsen2006). With the copula function K given in equation (26) which still depends on some abstract function φ we can construct a d-dimensional distribution function F by means of equation (25) which has the marginal distribution functions F ₁, …, F _d. Since, finally, we are interested in the joint density f of particle characteristics, applying equations (23)–(25) we get that

(27)$$\eqalign{&f\,(x) = f_1(x_1) \ldots f_d(x_d)\,k(F_1(x_1), \ldots, F_d(x_d)) \cr & {\rm for}\;x = (x_1, \ldots, x_d)\in {\rm {\opf R}}^d,$$

which can be seen as a differential variant of Sklar's theorem given in equation (25), where the function k : ℝ^d → [0, ∞) is the d-fold derivative

(28)$$k(u) = \displaystyle{{\partial ^d} \over {\partial u_1 \ldots \partial u_d}}K(u)\quad {\rm for}\;u = (u_1, \ldots, u_d)\in {\rm {\opf R}}^d.$$

Recall that the d-dimensional probability density f given by equation (27) has the marginal densities f ₁, …, f _d and, besides this, depends on the choice of a suitable copula function K or, equivalently, the choice of its derivative k. Thus, the task in modeling multidimensional probability distributions for given (1D) marginal densities is the determination of a suitable copula function K, which describes the correlation structure of the 1D marginals. Since in the present paper we only consider Archimedean copulas, which are defined by their generator φ, this task is equivalent to finding a suitable Archimedean generator.

Parametric Families of Archimedean Generators

Instead of analyzing single generator functions φ one by one, we can instead consider parametric families {φ _θ : θ ∈ Θ} of generators. For such families the function φ _θ is an Archimedean generator for each parameter θ ∈ Θ, where Θ is some set of admissible parameters. This has the advantage that we can consider a whole range of generators for which it is rather easy to determine an optimal choice of a generator among the family {φ _θ : θ ∈ Θ} for modeling the joint density f.

An example of a one-parametric family of generators is the Ali-Mikhail-Haq generator, see Nelsen (Reference Nelsen2006), which is given by

(29)$$\varphi _\theta (u) = \log \left( {\displaystyle{{1-\theta \,(1-u)} \over u}} \right),$$

for each u ∈ [0, 1] and some θ ∈ Θ = ( − 1, 1). Note that the parameter θ regulates the correlation between the uniformly distributed marginals of the corresponding copula K _θ. The relationship between the copula parameter θ and the common correlation coefficients, namely, the Pearson, Kendall, and Spearman correlation coefficients (Nelsen, Reference Nelsen2006), is shown in Figure 8 for 2D Ali-Mikhail-Haq copulas.

Fig. 8. Relationship between the copula parameter θ and the Pearson, Kendall, and Spearman correlation coefficients of the uniformly distributed marginals of 2D Ali-Mikhail-Haq copulas.

Note, that in the present paper, other one-parametric families of copulas are considered as well, e.g., the Frank copula whose generator is given by

(30)$$\varphi _\theta (u) = -\log \left( {\displaystyle{{{\rm e}^{-\theta u}-1} \over {{\rm e}^{-\theta} -1}}} \right),$$

for θ ≠ 0. An extensive list on even further copulas, such as the Clayton, Gumbel, Joe, and Plackett copula, can be found in Joe (Reference Joe1997). Beyond that, it is possible to obtain further copulas by rotating a given 2D copula by a multiple of 90°. Thus, there are four rotated versions of the aforementioned parametric families of copulas, which increases the list of considered parametric copula families even further.

Fitting a model from a parametric family of generators {φ _θ : θ ∈ Θ} is equivalent to determining an optimal parameter $\hskip2.5pt \hat{\hskip-2.8pt\theta} $. This problem can be described as an optimization problem using maximum-likelihood methods in the following way.

The family of generators {φ _θ : θ ∈ Θ} induces a family of copulas {K _θ : θ ∈ Θ} given by

(31)$$\eqalign{&K_\theta (u) = \varphi _\theta ^{-1} (\varphi _\theta (u_1) + \cdots + \varphi _\theta (u_d)) \cr & {\rm for}\;u = (u_1, \ldots, u_d)\in [0,1]^d,$$

for which the optimal parameter $\hskip2.5pt \hat{\hskip-2.8pt\theta} $ has to be determined. Assuming that K _θ is differentiable we can consider the d-fold derivative

(32)$$k_\theta (u) = \displaystyle{{\partial ^d} \over {\partial u_1 \ldots \partial u_d}}K_\theta (u),$$

and the corresponding d-dimensional distribution function F _θ, which is given by equation (25), has the probability density function

(33)$$\eqalign{& f_\theta (x) = f_1(x_1) \ldots f_d(x_d)k_\theta (F_1(x_1), \ldots, F_d(x_d)) \cr & {\rm for}\;x = (x_1, \ldots, x_d)\in {\rm {\opf R}}^d,$$

where f ₁, …, f _d and F ₁, …, F _d are the 1D densities and distribution functions fitted in the “Modeling of Single Particle Characteristics” section to the image data considered in the present paper.

Note that, for a sample of ℓ ∈ ℕ observations x ⁽¹⁾, …, x ^(ℓ) ∈ ℝ^d of the considered d-dimensional vector of particle characteristics, the so-called log-likelihood function is given by

(34)$$\eqalign{&\log {\rm {\cal L}}(\theta \vert x^{(1)}, \ldots, x^{(\ell )}) = \mathop \sum \limits_{k = 1}^\ell \log \,(\,f_\theta (x^{(k)})) \cr &\, {\rm for}\;\theta \in \Theta.$$

This leads to the maximum-likelihood estimator

(35)$$\hskip2.5pt \hat{\hskip-2.8pt\theta} = \mathop {{\rm argmax}}\limits_{\theta \in {\rm \Theta}} \;\log \,{\rm {\cal L}}(\theta \vert x^{(1)}, \ldots, x^{(\ell )}),$$

which can be considered as the optimal choice for the parameter θ when the sample x ⁽¹⁾, …, x ^(ℓ) was observed. The fitted joint density f is then given by $f = f_{\hskip1.5pt \hat{\hskip-2pt\theta}}. $

Until now, the procedure discussed above solely explains how to estimate the copula parameter θ for a given parametric family of copulas. In order to select an adequate model type among several parametric copula families, one can use the so-called Akaike information criterion (AIC) which utilizes the log-likelihood function considered in equation (34) to compare parametric families of copulas and select the best performing model. Note that, besides utilizing the log-likelihood function considered in equation (34), the AIC additionally takes into consideration the number of model parameters to avoid overfitting models. For more details regarding the AIC and model selection, see e.g., Held & Bové (Reference Held and Bové2014).

Application to Particle Characteristics

In the context considered in the present paper, the estimation procedure described above can be applied in the following way.

In order to obtain the (2D) joint density of the sphericity and median grayscale value of zinnwaldite-composite particles, see Figure 5, we utilize the 2D vectors of particle characteristics x ^(j) = (s(P _j), m(P _j)) which have been computed in the “Particle Characteristics” section for each zinnwaldite-composite particle P ₁, …, P _ℓ ∈ Z that hits the SEM-EDS plane. For both characteristics the parameters of the corresponding 1D probability densities $f_s^{\rm z}, f_m^{\rm z} $ for zinnwaldite-composite particles are given in Table 3. Now we consider, for example, the family of Frank copulas, assuming that it is able to describe the 2D density of the sphericity and median grayscale value. Then, we utilize the formula given in equation (33) and the previously determined 1D probability densities $f_s^{\rm z}, f_m^{\rm z} $ to obtain a parametric family of possible candidates for the 2D joint density of the considered particle characteristics. By inserting the computed vectors of particle characteristics x ^(j) into equation (34), we obtain the log-likelihood function $\log {\rm \;} {\rm {\cal L}}$, see Figure 9, which after maximization provides us the copula parameter $\hskip2.5pt \hat{\hskip-2.8pt\theta} = 2.45$. Thus, by means of equation (33) we obtain the 2D probability density of zinnwaldite-composite particle characteristics:

(36)$$f_{s,m}^{\rm z} (x_1,x_2) = f_s^{\rm z} (x_1)f_m^{\rm z} (x_2)k_{\hskip0.8pt \hat{\hskip-1.8pt\theta}} (F_s^{\rm z} (x_1),F_m^{\rm z} (x_2)),$$

where $F_s^{\rm z} $ and $F_m^{\rm z} $ are the distribution functions corresponding to the densities $f_s^{\rm z} $ and $f_m^{\rm z} $, respectively. Note that we found that the family of Frank copulas indeed provides a good fit, in comparison with other parametric families, by comparing the AIC of various families of copulas mentioned in the “Parametric Families of Archimedean Generators” section. Similarly, we showed that the Frank copula rotated by 90° provides a good fit for the 2D probability density $f_{s,m}^{\rm q} $ of the sphericity and median grayscale value of quartz-composite particles with the copula parameter $\hskip2.5pt \hat{\hskip-2.8pt\theta} = -1.14$, see Figure 5.

Fig. 9. Log-likelihood functions for fitting 2D copulas to the sphericity and median grayscale value of zinnwaldite-composite particles using a Frank copula (left) and of quartz-composite particles using a Frank copula rotated by 90° (right).

Recall that the quality of the decision, whether a particle mainly comprises zinnwaldite or quartz, improved significantly once we considered the 2D probability densities $f_{s \comma m}^{\rm z},\, $$f_{\vskip0.01pt s,m}^{\rm q} $ over 1D densities, compare Tables 1 and 2. Furthermore, recall that the estimation procedure described in the “Parametric Families of Archimedean Generators” section is not limited to the simultaneous consideration of two particle characteristics. But, in the same way, it can be applied for general d-dimensional vectors of particle characteristics. For example, we can investigate the six particle characteristics introduced in the “Particle Characteristics” section, thus considering the case d = 6.

For simplicity, we first fit one-parametric copulas to the six-dimensional (6D) data, which leads to better prediction results than just considering 2D particle characteristics, see Table 4. However, we are aware of the fact that one-parametric copulas do not entirely capture the correlation structure of the 6D data, because the empirical (pairwise) correlation coefficients can vary distinctly for each individual pair of particle characteristics. It turned out that, among the one-parametric copula families mentioned in the “Parametric Families of Archimedean Generators” section, the Ali-Mikhail-Haq copula yielded the best fit. In particular, using both the 6D vectors of characteristics x ^(j) = (r(P _j),  s(P _j),  c(P _j),  e(P _j),  m(P _j),  iqr(P _j)) of zinnwaldite particles P ₁, …, P _ℓ ∈ Z and the corresponding 1D probability densities $f_r^{\rm z}, f_s^{\rm z}, f_c^{\rm z}, f_e^{\rm z}, f_m^{\rm z}, f_{{\rm iqr}}^{\rm z} $, see the “Modeling of Single Particle Characteristics” section, by means of equations (31)–(34) we obtain the copula parameter $\hskip2.5pt \hat{\hskip-2.8pt\theta} = 0.69$, and thus a parametric fit for the 6D density f _z of zinnwaldite-composite particle characteristics.

Table 4. Confusion Matrices of the Copula-Based Decision Rule Given in Equation (41) for Particles Observed in the SEM-EDS Section Which Has Been Used for Adjusting the Prediction Model and for the Particles Observed in the SEM-EDS Section Which Was Withheld for Validation.

Analogously, we obtained the Ali-Mikhail-Haq copula with parameter $\hskip2.5pt \hat{\hskip-2.8pt\theta} = 0.47$ and thus the multidimensional density function f _q for the characteristics of quartz particles.

One way to achieve even better prediction results would be to utilize multi-parametric copula families to fully capture the correlation structure of the 6D data, such as, e.g., vine copulas, see Bedford & Cooke (Reference Bedford and Cooke2002). However, in the following, we consider yet another approach to achieve this goal, which uses the 6D data, but projects it to a suitably chosen 2D subspace. In this way, we still exploit information of the 6D vector of particle characteristics while being able to consider 2D copulas, or even 2D kernel density estimators.

Kernel Density Estimation and Dimension Reduction

Recall that for describing multivariate probability densities of particle characteristics, which are necessary for improving the decision rule of the 1D case considered in equation (14), we used a parametric copula approach in the “Application to Particle Characteristics” section. Alternatively, such densities could be estimated in a nonparametric way, using the so-called kernel density estimators (Scott, Reference Scott2015). To be precise, from a sample x ⁽¹⁾, …, x ^(ℓ) ∈ ℝ of a 1D characteristic of zinnwaldite-composite particles a non-parametric density $\hskip1pt\hat{\hskip-1pt f}_{\rm z}$, can be estimated by

(37)$$\hskip1pt\hat{\hskip-1pt f}_{\rm z}(x) = \displaystyle{1 \over \ell} \mathop \sum \limits_{\,j = 1}^\ell \kappa \left( {\displaystyle{{x-x^{(\,j)}} \over h}} \right)\quad {\rm for}\;x\in {\rm {\opf R}},$$

where h >0 is called the bandwidth and κ : ℝ → [0, ∞) is a kernel function, e.g.,

(38)$$\kappa (x) = \displaystyle{1 \over {\sqrt {2\pi}}} \exp {\rm \;} \left( {-\displaystyle{{x^2} \over 2}} \right)\quad {\rm for}\;x\in {\rm {\opf R}}.$$

Similarly, it would be possible to estimate non-parametric multivariate densities for vectors of particle characteristics. However, note that due to the “curse of dimensionality” the required sample size for adequate kernel density estimation increases exponentially with the considered dimension, see Scott (Reference Scott2015). Thus, it is difficult to utilize kernel density estimators for computing accurate 6D densities representing our 6D data. Nevertheless, the issue concerning the sample size of the 6D vectors of particle characteristics, caused by the “curse of dimensionality,” can be remedied by dimension reduction. Common techniques such as the principal component analysis (PCA) project the data to a lower dimensional subspace which maximizes the variance of the data (Hastie et al., Reference Hastie, Tibshirani and Friedman2009). However, PCA does not ensure a good separability of the two classes, namely zinnwaldite- and quartz-composite particles, after projection of the data. Other methods, such as linear discriminant analysis (LDA) can reduce the dimension of the data such that classification can still be performed reliably. Note that LDA always projects data on a c − 1 dimensional subspace, where c is the number of classes. In the setting of the present paper, LDA would project the 6D vectors of zinnwaldite- and quartz-composite particle characteristics on a 1D subspace such that the normalized discrepancy measure

(39)$$\Delta _1 = \displaystyle{{{(\mu _{\rm z}-\mu _{\rm q})}^2} \over {\sigma _{\rm z}^{\,2} + \sigma _{\rm q}^{\,2}}} $$

is maximized, where μ _z, μ _q ∈ ℝ are the means of the projected characteristics vectors of zinnwaldite and quartz, respectively. Analogously, the values $\sigma _{\rm z}^{\,2}, \sigma _{\rm q}^{\,2} \in {\rm {\open R}}$ are their variances. Heuristically speaking, the LDA maximizes the distance between clusters, while minimizing their variances. In order to fit multidimensional probability densities or to perform multidimensional kernel density estimation on the data after dimension reduction we modified the expression considered in equation (39) such that we reduce the data to two dimensions. More precisely, we project the 6D data on a 2D subspace such that the normalized discrepancy measure

(40)$$\Delta _2 = \displaystyle{{\parallel {\tilde{\mu}} _{\rm z}-{\tilde{\mu}} _{\rm q}\parallel ^2} \over {\sigma _{{\rm z},1}^{\,2} + \sigma _{{\rm z},2}^{\,2} + \sigma _{{\rm q},1}^{\,2} + \sigma _{{\rm q},2}^{\,2}}} $$

is maximized, where $\parallel \cdot \parallel $ denotes the Euclidean norm in ℝ² and $\tilde{\mu} _{\rm z},\tilde{\mu} _{\rm q}\in {\opf R}^2$ are the means of the clusters corresponding to zinnwaldite- and quartz-composite particles after projection. The values $\sigma _{{\rm z},1}^{\,2}, \sigma _{{\rm z},2}^{\,2} $ are the respective variances of the first and second components for vectors of characteristics of zinnwaldite after projection. Analogously, $\sigma _{{\rm q},1}^{\,2}, \sigma _{{\rm q},2}^{\,2} $ denote the variances of the components for the quartz case. By maximizing the expression given in equation (40) we obtain the dimension reduction by projecting for each particle P the corresponding vectors of characteristics x = (r(P), s(P), c(P), e(P), m(P), iqr(P)) onto the plane with span vectors v ₁ = (0.0014, 0.3189, 0.0189, 0.0562, − 0.9455, − 0.0285) and v ₂ = (0.0001, − 0.0944, 0.1410, − 0.0193, − 0.0598, 0.9835). Figure 10 indicates that the vectors of characteristics for zinnwaldite- and quartz-composite particles after projections are relatively distinct. Furthermore, this reduction of dimension makes kernel density estimation more viable again. Note that the “ragged” probability densities in Figures 11a and 11b still indicate too small sample sizes for adequate kernel density estimation. However, the methods discussed in the “Modeling of Single Particle Characteristics” and “Parametric Families of Archimedean Generators” sections can be used to fit parametric copulas to the 2D data after dimension reduction, see Figures 11c and 11d.

Fig. 10. 6D data after dimension reduction obtained by maximizing cluster distance, while minimizing cluster variances.

Fig. 11. Kernel density estimation of vectors of characteristics after dimension reduction for (a) zinnwaldite- and (b) quartz-composite particles. The corresponding probability densities obtained by fitting parametric copulas for (c) zinnwaldite- and (d) quartz-composite particles.