Image Anal Stereol 2001;20:101-104 Original Research Paper PREDICTION OF THE EXTREMAL SHAPE FACTOR OF SPHEROIDAL PARTICLES Daniel Hlubinka and Viktor Benes¡ Department of Probability and Statistics, Charles University, Sokolovska´ 83, 186 75 Praha 8, Czech Republic E-mail: hlubinka@karlin.mff.cuni.cz (Accepted June 1, 2001) ABSTRACT In the stereological unfolding problem for spheroidal particles the extremal shape factor is predicted. The theory of extreme values has been used to show that extremes of the planar shape factor of particle sections tend to the same limit distribution as extremes of the original shape factor for both the conditional and marginal distribution. Attention is then paid to the extreme shape factor conditioned by the particle size. Normalizing constants are evaluated for a parametric model and the numerical procedure is tested on real data from metallography. Keywords: extremal shape factor, spheroidal particles, unfolding. INTRODUCTION Consider a random collection of spheroidal particles in a reference volume. We shall assume that the particles are all oblate or all prolate. The size is the length of a particle’s major semiaxis X, while W is the length of the minor semiaxis. The shape factor in our setting is given by T = X2/W2 — 1. It is clear that 0 < W < X < T], and 0 < T < (O, where the equalities X = W and T = 0 hold when the particles are balls. Both X and W are nonnegative real random variables. Values r\ and (O are fixed non-negative real numbers, possibly infinite, the upper end-points of the supports of distribution of X and T respectively. We study random planar sections of the particles. These sections form ellipses. An ellipse is characterized by its size (major semiaxis length) Y, and its shape factor Z = Y2/V2 — 1, where V is the length of the minor semiaxis. It holds that 0 < Y < X, and 0 < Z < T. Let g(x, t) be the joint probability density function of size and shape factor (X, T) of a particle. The orientation of the particle is assumed to be isotropic, i.e. a uniform random variable on the hemisphere independent of size and shape. On the other hand we do not assume independence between size and shape factor. Following (Cruz-Orive, 1976) the distribution of the section size and shape factor (Y,Z) for oblate particles has the joint probability density f(y,z) 2M r\ pa> yz g(x,t)dtdx y/ty1+tyt zy/x2' -y2 (1) where 2M is the population mean caliper diameter of particles. We will need the joint distribution of the original size X and planar shape factor Z as well. Using similar arguments we get the probability density h(x,z) 1+ /1+z I0 Mx Jz g(x, t) dt where Mx = yty1+tyftz dz dt (2) arctan \ft and gx(t) = g(x,t)/z®g(x,t)dt is the conditional density of T given X = x. Then the complementary equation which together with (2) yields (1) is f(y,z) y n Mxh(x,z) 2M \fx2 -y dx. (3) We denote fy(z) the conditional density of shape factor Z given size Y = y. The conditional probability density of Z given X = x is equal to r gx(t)dt z t 1+tt Vty1+ty/t-z' (4) Further let h(z) denote the marginal density of the transformed shape factor Z. We have h{z) = 1 + z 1 yfty1+tytz 0 Mx Finally define fz(y), hz{x) as the conditional density of size Y given shape Z = z, and size X given shape Z = z, respectively. 101 Hlubinka D, Benes¡ V: Extremal shape factor For the sake of completeness we provide the formula for the density of the section characteristics for prolate particles which is expressed in terms of shorter semiaxes fv z 2M(1+z) {1+t)3/2g(w,t)dtdw yftyjt - zy/w2 - v2 In the following, however, we restrict consideration to oblate particles. For the prolate particle case analogous results hold. The theorem says that the domain of attraction for the shape factor is the same for a particle and its section. We will not study further part c) since it is less important for practical applications. Indeed, the total extremal shape factor is useless if it is not related to particle size. To use the results of theorem 1, part a), we need to estimate normalizing constants an,bn,an,b'n such that conditioned by the known particle’s size X = x T(n) bn 9 A, EXTREMAL PROPERTIES Z(n) bn 9 A, A univariate distribution function K is said to belong to the domain of attraction of a distribution function L if there exist normalizing constants {an}, {bn} such that the n-th power Kn{anx + bn) converges weakly to L, where L is one of the following distributions Li,r(x exp(- -x~ ~r), x>0, i exp(- -(- -x)" -r), x<0, i exp(- -e~ "x), xel, i 1, Fr echet 2, Weibull 3, Gumbel (5) where y > 0. We shall write Kg S>{L) if K is in the domain of attraction of L. If the distribution K has a density k, there are sufficient conditions (de Haan, 1975) for K to be in @(L). These conditions are denoted (C1 ), (C2 ), (C3). In Drees and Reiss (1992) the extremal properties in the Wicksell corpuscule problem were studied. This was developed so as to apply to the size-shape unfolding problem of spheroidal particles. Let Hx, H, Gx, Fy, Fz, Hz be the distribution functions corresponding to densities hx, h, gx, fy, fz, hz respectively. Then the following theorem holds (Hlubinka, 2000): Theorem 1 a) Suppose that for any fixed size x the density gx(t) satisfies condition Ci . Then the distribution Hx G &(Li „), i = 1,2,3, where ß = J for i = 1, and ß = y+ 1/2 for i = 2. i,y b) Assume that gx(t) satisfies condition Ci uniformly in x. Then Fy G 3>(Li „) for all y, i = 1,2,3 where ß = J for i = 1, and ß = y+ 1/2 for i = 2. c) Assume that gx{t) satisfies condition Ci uniformly in x. Then H G ^(Li„), i = 1,2,3, where ß = jfor i = 1, and ß = y+ 1/2 for i = 2. where T , is the maximum of n observations, and (n) A is a random variable with the extremal type of distribution. A parametric model of gamma distribution will be considered with the density M—1 ..a gx(t;a,ß) = -^-e-^,t>0. (6) where a and \i are positive constants which can possibly depend on the given size x. Recall that the gamma distribution is in the domain of attraction of the Gumbel distribution. Using the techniques of Takahashi (1987) the following theorem was proved in Hlubinka (2000): Theorem 2 Assume that for any fixed size X = x the shape factor T follows a gamma distribution. Then both T and Z conditioned on X = x belong to the Gumbel domain of attraction and their normalizing constants are 1 1 M 1 M 1_ r(«) logn + (a— 1)loglogn + log logn + | a - - ) loglogn + log 2) MxT(a) Using Theorem 2 we can approximate the distribution Gnx(t) ofT, by A((t — bn)/an). Quantiles of T,s are estimated via qp = bn + an[— log(—logp)], and ET,, = bn + anC, where C = 0.5772 is the Euler constant. Concerning the estimation of normalizing constants an,bn using estimation of an,bn, several methods concerning the corpuscule problem are 102 Image Anal Stereol 2001;20:101-104 discussed in Takahashi and Sibuya (1998). We can use for ecample a maximum likelihood method based on the joint distribution of {Tn_k+1-,,..., T,,) which is derived in Weissman (1968) to estimate an and b'n. Note that since an = an = 1j\i does not depend on n the main problem is to find a reasonable estimation of bn based on 1//x and b'n. This means the use of numerical methods since a is unknown. For the situation b) from theorem 1, unfortunately, such a simple result as Theorem 2 is not available. PRACTICAL APPLICATION The prediction of the extremal particle shape factor is important in metallography. Engineers claim that damage to construction materials depends on extremal rather than mean characteristics of the microstructure. Therefore corresponding methods for prediction of spatial extremes based on data from sections (metallographic samples) have to be developed. In Benes¡ et al. (1997) the problem of unfolding the size-shape-orientation distribution has been theoretically solved. We use data from that study, i.e. measurement of an aluminium alloy specimens with oblate Si particles. The assumptions of spheroidal shape and isotropy are approximately fulfilled. Metallographic samples were measured by means of an image analyser to obtain samples from the joint distribution with density f(y,z). The aim is to use Theorems 1 and 2 for the prediction of the extremal shape factor of particles for given size classes. The gamma distribution parametric model (6) for the shape factor given the size is used. Our approach is to use formula (3) for unfolding the size distributionof particles to obtain an empirical distribution corresponding to the density h(x,z). This is then fitted to the theoretical distribution obtained by plugging (6) into (4). This serves two purposes: first we can evaluate the degree of fit and conclude whether the model is appropriate, secondly the parameters a and \i are estimated and normalizing constants an, bn obtained from Theorem 2. In the first step we use discretization and the EM-algorithm for the unfolding as described generally e.g. in Ohser and Mu¨cklich (2000). Equation (3) is transformed to: 1 n 1-Fz(y) = 2 My p(x,y)hz(x)dx, (7) where p(x,y) = Mx^/x2-y2. It holds that where NV is the mean number of paticles per unit volume and NA is the mean number of particle sections per unit area of the section plane. Equation (7) is thus transformed to NA{1-Fz(y)) =NV yVp(x,y)hz(x)dx, (8) where NA is estimated from the data. The y variable is subdivided into classes with endpoints yk = bk, b a given constant. Then for t,k = NA (Fz (yk) - Fz (yk_ 1)) we have rn Ck = NV [pxyk-1) -p(xyk)] hz(x)dx. Jy Assume that Hz(x) can be discretized to have values ai for xi1 < x < xi, where xi = bi, then k = N V ^p ik(a i - a i_1) = 5paL-, where pik = p{xiyk_l) - p{xlyk)J i > k, pik = 0 otherwise. It holds that 1 i = NV i ai — ai_1), from which the estimator of NV is NV = X, %i Now the iteration of the EM-algorithm is given by iv + 1)—st step: &+1=q-i pck (9) 'ii k=1 r k where ^v is the v-th iterated ^, qi = Y^kpiki rk = ^i>kpik i i'•• and the initial value i,i 0 can be taken to be Li. In this way probabilities B,JNV of the discrete conditional distribution Hz{x) given z are estimated. Using the empirical distribution function H(z) = F(z), also the contrary conditional distribution function Hx(z) can be directly estimated from Hz{x). Then for the shape factor subintervals {zi_1zi denote He m = Hx(zi - Hx{zi_1) for fixed x. The corresponding theoretical probabilities Hjh for the gamma model are obtained from (4) and (6) as rzj 1 ,,a where z"00 ta1 e^ßt I2 = zj t 1 +t)Bt{zj-1zj)dt' 103 Hlubinka D, Benes¡ V: Extremal shape factor and the integral Bt(s,u) = su y/( 1+z)/{t — z)dz is given in the expression forMx (following (2)). Since in fact the unknown parameters play a role already in the unfolding step (through the Mx term in kernel function p), it is useful to suggest an analytical function modelling the dependence of the parameters \i, a on x. For a given size class we estimate the parameters ll, a by minimizing either the ^(Hi th — Hi em)2 or max, lY.Hi th-Y.<,Hi em!. "i