Image Anal Stereol 2002;21:175-181 Original Research Paper CHARACTERIZATION OF BIVARIATE SIZE-ORIENTATION DISTRIBUTION OF CIRCULAR PLATE PARTICLES Karel Bodläk1, Arun Balasundaram2, Arun M Gokhale2 and Viktor Beneš1 'Department of Probability and Mathematical Statistics, Charles University, Sokolovska 83, 18675 Praha 8, Czech Republic, 2School of Materials Science and Engineering, Georgia Institute of Technology, Atlanta, GA 30332-0245, USA e-mail: benesv@karlin.mff.cuni.cz (Accepted 26 August, 2002) ABSTRACT The paper is devoted to the stereological unfolding problem of bivariate size-orientation distribution of platelike particles in metallography. Gokhale (1996) derived an integral equation which relates this bivariate distribution in three-dimensional (3D) space to the corresponding size-orientation distribution of planar sections of the specimen. The present paper yields a numerical algorithm which enables to transform a bivariate histogram of observed quantities to the histogram of 3D characteristics. The use of the method is demonstrated in examples with simulated data, where an easy analytical solution is available and can be compared with the results of estimation. The spectrum of unfolding problems solved numerically (cf. Ohser and Muecklich, 2000) is thus extended. Keywords: circular plate particles, EM algorithm, unfolding VUR sampling. INTRODUCTION The stereological estimation of the joint distribution of particle parameters from planar sections enables to quantify the 3D geometrical properties of the materials microstructure which is related to mechanical properties. The bivariate distribution of the size and shape factor of spheroidal particles was studied by Cruz-Orive (1976) and later for particles of some other shapes in Ohser and Muecklich (2000). The trivariate size-shape-orientation unfolding problem for spheroidal particles was solved in Beneš and Krejèif, (1997). Two main sampling schemes were used, either the IUR (isotropic uniform random) or VUR (vertical uniform random) sections. Particles or cracks which can be modelled by circular plates have zero shape factor. In this case the investigation of size-orientation distribution is sufficient, which was investigated first by Gokhale (1996). In that paper the theoretical solution is presented for the VUR sampling design, i.e., the integral equation connecting the probability densities is derived as well as its inversion. This inversion theoretically enables to evaluate the bivariate distribution of the spatial characteristics (the size and orientation in 3D) from the planar bivariate distribution. However, the formula is based on two-fold integration followed by a two-fold differentiation and an analytical form of planar bivariate joint probability density is required. Given the bivariate histogram of section data measured by means of an image analyser this would need first a parametrical bivariate density model to be fit. Altogether it is not recommendable to proceed using inversion formula directly, therefore a numerical histogram-based method has to be developed which is known, e.g., for size-shape problems (Ohser and Muecklich, 2000). To do this the first step is to transform the integral equation in terms of probability densities to an equivalent integral equation in terms of distribution functions. This corresponds to the discrete data of the transformation of relative frequencies to cummulative frequencies. The integral equation for distribution function is then discretized which leads to a system of linear equations. Since a direct solution of this system may lead to negative frequencies, special numerical approaches have been adapted. Among them the EM-algorithm (Ohser and Muecklich, 2000) is the most popular one recently, since it always leads to a non-negative solution. First attempt to the solution of the formulated problem was published by Beneš et al. (1996). The integral Eq. 1 was split in two subsequent univariate unfolding problems which were discretized in its original form. Thus the use of elliptic integrals was unavoidable. In the present paper an improvement consists in a direct discretization of bivariate equation which decreases estimation errors of EM-algorithm. 175 BodläkKet AL: Size-orientation distribution of circular plate particles Further an ingenious transformation (Eq. 4) of data leads to a simpler evaluation of kernel coefficients in a closed form. In the next Section the derivation of the integral equation between distribution functions is presented which arises to be difficult for the size-orientation problem. A special solution which needs a transformation of functions was derived, so that finally one does not deal with a distribution function exactly, but with an alternative form which is convenient for the numerical solution. A simple transformation of data then enables a standard numerical solution using EM-algorithm (Ohser and Muecklich, 2000) in Section 3. In Section 4 a synthetic example is present which demonstrates the use of the developed algorithm. Model bivariate densities in Gokhale (1996) are tested such that the inverse solution is known theoretically. Then a sample is simulated from the planar bivariate distribution and transformed to the desired bivariate histogram estimating the spatial size-orientation distribution. True and estimated marginal distributions of spatial characteristics are plotted to demonstrate the results of the method. ESTIMATION OF BIVARIATE SIZE-ORIENTATION DISTRIBUTION BY EM ALGORITHM Fix an arbitrary direction in the space which is called a vertical axis. Consider a coordinate system in the 3D space, such that the vertical axis coincides with z-axis, and a system of circular plates dispersed in a reference volume. Each plate is described by its radius R and a normal orientation (0,. That means we will 2k 0j consider the distribution of size and colatitude only which arises in many applications, e.g. if a material specimen is under tension or compression in the direction of the vertical axis. Let f(r,oc) be the bivariate probability density of half chord length r and chord orientation a of plates if the planar section is VUR. The relationship between f(r,oc) and g(R,9) is given by a double integral equation in Gokhale (1996) \*'t cos20sin0g(i?,0)d0dR ' K) r n/\-a (R2 - S j'" (sÜl2 ÖT- COS2 fff* where Nv is the mean number of particles per unit volume, Na is the mean number of particle sections per unit area of a vertical plane, rm is the largest particle's radius. 176 Image Anal Stereol 2002;21:175-181 From practical reasons it is better to use y = n/2-a in f(r,oc). Then we get after this change of variables *(r,r) 4rN, #NAcos y Y't cos20sm0g(R,0)d0dR I I (R2-r2f (sin*y-cos20f (2) In order to present a numerical solution we transform the Eq. 2 from the relationship between the density functions to the relationship between the distribution functions. The procedure is described as follows. After the differentiation the right hand side of Eq. 2 transforms to relatively to the amount of data because otherwise the method can lead to large disturbances in an output. Assume that plate sections were observed using the VUR sampling design and their parameters were measured. Suppose now that a number M of size classes and a number N of orientation classes is used for discretization. Thus we can define class intervals with limits as follows: ri= R; = i8, 8 = rm/M, i = 1,...,M; Yj = 0j =ja>, a) = 7t/(2N),j = 1,...,N. Typically rm is unknown and it is estimated by means of a half of the largest observed plate trace length.Under this choice of the discretatization the platelike particle's trace with size and orientation r, y, respectively, belongs to an (ij)-th class if r;_i < r < r; andYj-i^Yr/2 \ [ \(R2 -r2)m(cos2y-cos20)mcos20sm0g(R,0)d0dR r y kj>! 71 (s-i)8 0 we obtain 177 BodläkK ET AL: Size-orientation distribution of circular plate particles js(R,r) dR = -\b(b2-r2)U2-a(a2-r2f - r2 log and for setting r = 0 we get b + (b2-r2f a + (a2-r2)1 1 fs(i?,r) dR = —\b2 -a2\. Integrating the right side J 9 of Eq. 8 in 6 we can write Jt(0,y) dG: cos 9 (cos2 y-cos2 6) cos y cos4y (cos2 y-cos2 9) arctan cos G (cos2 y-cos2 G) In order to solve the system of linear Eq. 6 the iterative EM algorithm in Ohser and Muecklich (2000) is used. Suppose that y is said to be the matrix of incomplete data obtained by planar sampling and x represents the matrix of parameters to be estimated. Each iteration consists of two steps, an E-step and an M-step, where E stands for expectation and M for maximization. Let z^ be the number of events occuring in the class (i j) which contributes to the counts in the class (k,l). The E-step yields the expected value of Zijki given input data y under the current estimate x( ^ of the parameter matrix x as follows ^=xfyuPim'ru')- (9) The M-step yields a maximum likelihood estimate of the parameter x using the estimated data zp from the M N M N E-step.Define ttj =YLPw md r^ = E5>#hx, k=l 7=1 (X) ij i=l j=\ Then the M-step is given by ,w i M N fzz zijkl. Combining these two steps we get an EM iteration given by the updating formula yaPijU r(X+\) fc=l /=1 ru (X) (10) i = l,...,M,j = l, ,N The Eq. 10 produces a sequence x„' which W converges to a solution, i.e., the matrix of parameters (°) — i> This choice yy x. As an initial iteration we setx~ of nonnegative initial values ensures that each of x is nonnegative. There is in fact not a unique solution. Convergence of the EM algorithm is guaranteed in theory but the solution can depend on the choice of the initial values used in the iteration scheme. See Ohser and Muecklich (2000) for more experience with the use of the EM algorithm in metallographical practice when solving unfolding problems. SIMULATION STUDY The following examples are given to illustrate the use of the EM algorithm for estimation of the true bivariate size and orientation distribution from simulated planar data according to given joint density functions. The number of iterations of the EM algorithm depends on how fine is the discretization as well as on the accuracy of measured planar data. For a sample size about 10000 and a total of 10 classes the number of iterations need not exceed 32 as recommended in Ohser and Muecklich (2000). We have chosen 5 classes of the radii and 5 classes of the orientation. The number of the simulated platelike particle traces was 2000 and the number of iterations of the EM algorithm was suggested equal to 16. Suppose first that the platelike particle's size and orientation distribution is fi(r,a) = U-lO^/Ti-r^-r2)2, the radius of the largest platelike particle is given by rm = 0.001 cm. The size of a class is equal 0.0002 cm. The derived true platelike particle's size and orientation distibution (Gokhale, 1996) has a form gi(R,9) = 1016/7iR(Rm2-R2)3/2 where Rm = rm, i.e. orientations are uniform random in this example. If we set the number of platelike particles per cm2 of the vertical plane area at value Na = 104, then the total number of platelike particles Nv is equal to 64/ttIO6. Using the EM algorithm to evaluate the data simulated from the density f we get the estimated bivariate size and orientation distribution est(g) in 3D, Fig. 3. The estimated total number of microcraks per a unit volume is est(Nv) = 2.04-107 which is very close to the true value Nv. The true and estimated marginal distribution of the radius and the orientation is plotted in Fig. 4. In order to compare the true and estimated marginal distributions we have to transfer a scale of the true distribution to classes from 1 to 5 on the horizontal axis and use a relative scale for the vertical axis. 178 Image Anal Stereol 2002;21:175-181 Fig. 3. Platelike particle's distribution with the uniform distribution of the angle, (a) The bivariate histogram given by simulated planar data from f. (b) Estimated bivariate distribution in 3D. Fig. 4. The solid lines (T) show the true marginal distributions of the spatial size and orientation distribution corresponding to Fig. 3. We indicate smoothed estimated marginal distribution corresponding to the size and orientation by dotted line (E). (a) The curve E is fitted by the spline of the third order, (b) For smoothing the curve E we use the least squares method. In the next example suppose that platelike particle's trace r and orientation a are still uncorrected but the orientations are not uniform random. The planar probability density has a form f2(r,oc) = 24-1018/7i-r(rm2-r2)2 sin2oc and true platelike particle's size and orientation distribution is g2(R,6) = 2-1016/7tR(Rm2-R2)3/2 cos26. For the same setting of NA the total number of platelike particles Nv is equal to 256/(37i)106. From the simulation we receive the estimated value est(Nv) = 2.8-107. The estimated bivariate distribution in 3D is depicted in Fig. 5, analogously to Fig. 4 the true and estimated marginal distribution of 3D characteristics is plotted in Fig. 6. The last simulation example covers the situation of correlated platelike particle's trace r and its orientation angle a. The planar probability density is given by f3(r,a) = 4-1012/7i-r(rm2-r2)+12-1018/7tr(rm2-r2)2 sin2a. The associated distribution in 3D has a form g3(R,6) = 3- 109/7t-R(Rm2-R2)1/2 + 1016/rc-R(Rm2-R2)3/2 cos20 and the associated total number of platelike particles is Nv = 16-107/(6rc) + 128-106/(3ti). The estimated total number of platelike particles is est(Nv) = 2.MO7 which is again very close to the true value Ny. The simulated planar and estimated spatial bivariate distributions are showed in Fig. 7, the true and estimated spatial marginal distribution is in Fig. 8. Fig. 5. Anisotropic uncorrelatedplatelike particle's distribution, (a) The bivariate histogram given by simulated planar data from f2. (b) Estimated bivariate distribution in 3D. 179 Bodläk K ET AL: Size-orientation distribution of circular plate particles Fig. 6. The solid lines (T) show the true marginal distributions of the spatial size and orientation distribution corresponding to Fig. 5. We indicate the smoothed estimated marginal distribution corresponding to the size and orientation by dotted line (E). The both curves E in (a) and (b) are fitted by the spline of the third order. Fig. 7. Anisotropic and correlated platelike particle's distribution, (a) The bivariate histogram given by simulated data from fa. (b) Estimated bivariate distribution in 3D. Fig. 8. The solid lines (T) show the true marginal distributions of the spatial size and orientation distribution corresponding to Fig. 7. We indicate smoothed estimated marginal distribution corresponding to the size and orientation by dotted line (E). (a) The curve E is fitted by the method of least squares, (b) For smoothing the curve E we use by the spline of the third order. DISCUSSION A method for the estimation of the spatial bivariate size-orientation distribution of circular plates using stereological methods was developed. The unfolding of true plate parameters is based on the observation in a vertical uniform random section plane. In metallographical practice the circular plates present an approximate model for cracks or platelike particles in particle reinforcement. The vertical axis then coincides with the direction of tension or compression of the material. If a close analytical form of the true spatial bivariate platelike particle's distribution is available then (sometimes after tedious mathematical calculations) one could express the associated planar bivariate distribution in the vertical planes and vice versa. But for a real dataset it is difficult to fit its planar bivariate distribution by an analytical formula. In such a case (especially when the data are presented in a bivariate histogram form) it is easier to proceed in discretization and apply an efficient numerical procedure like EM algorithm to calculate the bivariate Image Anal Stereol 2002;21:175-181 distribution of 3D characteristics from the measured data. Such a method was developed in the present paper and its use was demonstrated in synthetic examples where the data were simulated from a known distribution. The results are promising for the application of the procedure in metallographical practice, they extend the supply of unfolding methods discussed in Ohser and Muecklich (2000). ACKNOWLEDGEMENTS Karel Bodläk acknowledges U.S. NSF sponsored programme (DMR-9816618) and Czech Ministrery of Education project KONTAKT ME335 for international collaboration between Charles University and Georgia Institute of Technology for financial support. He also thanks Georgia Institute of Technology for allowing his stay in Atlanta. A.B. and A.M.G acknoledge U.S. National Science of Foundations Research Grant No. DMR-9816618 (for which Dr MacDonald was the Program Monitor) for financial support. V.B. was supported by project MSM11300008 of the Czech Ministery of Education. REFERENCES Beneš V, Krejèif P (1997). Decomposition in stereological unfolding problems. Kybernetika 33(3):245-58. Beneš V, Gokhale AM, Slämovä M (1996). Unfolding the bivariate size-orientation distribution. Acta Stereol 15(1):9-14. Cruz-Orive LM (1976). Particle size-shape Distributions: The General Spheroid Problem. J Microsc 107(3):235-53. Gokhale AM (1996). Estimation of bivariate size and orientation distribution of platelike particles. Acta Metall and Mater 44(2):475-85. Ohser J, Muecklich F (2000). Statistical Analysis of Microstructures in Material Science. New York: J Wiley & Sons. 181