Image Anal Stereol 2006;25:55-61 Original Research Paper PRACTICAL DETERMINATION OF THE COORDINATION NUMBER IN GRANULAR MEDIA Patricia Jouannot-Chesney1, Jean-Paul Jernot1 and Christian Lantuéjoul2 1ESCTM du CRISMAT, UMR 6508, ENSICAEN, 6 Boulevard du Maréchal Juin, 14050 Caen Cedex, France; 2Ecole des Mines, Centre de Géostatistique, 35 rue St Honoré, 77305 Fontainebleau, France e-mail: jouannot@ensicaen.fr; jernot@ensicaen.fr; christian.lantuejoul@ensmp.fr (Accepted March 3 2006.) ABSTRACT In a random stacking of particles, the coordination number is defined as the mean number of particles touching a given one. A classical measurement by image analysis is directly based on this definition. It can be applied in the theoretical case of a global analysis but some problems occur for practical applications. We propose here a new measurement method based on the contribution to the Euler-Poincaré characteristic. A single measurement is performed both on the initial structure (particles linked by contacts) and on the segmented structure (isolated particles). The proposed method is robust w.r.t. the shape of the particles and their contacts. Keywords: coordination number of particles, Euler-Poincaré characteristic, granular media, local measurements. INTRODUCTION The concept of coordination number arises in various contexts, sometimes along with different acceptations: – in a monoatomic structure, the coordination number is nothing but the number of neighbours of each atom (cf. Fig. 1a) – in a random stacking of particles, different particles may have different numbers of neighbours; in this case, it is more appropriate to call «coordination number» the average number of neighbours per particle. The coordination number plays an important role in the description of granular media. It can be related to the microstructural evolution of a powder during sintering (cf. Fig. 1b) and to several physical or mechanical properties of sintered materials (Artz, 1982; Jernot, 1991; Tancret et al., 1997). It has been estimated in the past by using several methods (German, 1989; Guyon and Troadec, 1994). We will focus here on its direct measurement by image analysis without any assumption about the particle shape. a) b) Fig. 1. Contacts in granular media: a) punctual contacts in a periodical arrangement of particles, b) necks between sintered bronze particles. 55 Jouannot-Chesney P et al: Determination of the coordination number GLOBAL MEASUREMENTS The coordination number cannot be obtained using stereological relationships because it is a topological parameter. As such, it must be measured directly in the space of interest (e.g., mostly 3D space in materials science). The coordination number is defined as the mean number of particles touching a given one. Touching particles are called neighbours. In what follows, two global measurement methods are presented. Sequential analysis This method directly derives from the definition of neighbours. It starts with a segmentation by watershed of the population of particles (Lantuéjoul and Beucher, 1981). Then, a sequential analysis is performed on this segmented set. At first, one particle is selected and reconstructed by geodesic dilation (Fig. 2a). Then, after dilation, this particle is intersected with the segmented set and its number of neighbours is counted (Fig. 2b). This procedure is repeated for each particle up to a complete examination of the initial set (note that it is time-consuming). The coordination number is the mean number of neighbours measured for all the particles (Fig. 2c). Parallel analysis This method is based on the measurement of the Euler-Poincaré characteristic (Hadwiger, 1957; Serra, 1982; Stoyan et al., 1995; Ohser and Nagel, 1996) further denoted by EPC. In 2D space, this topological parameter is simply the number of connected components minus the number of holes within them. In 3D space, it is equal to the number of distinct surfaces minus their genus. Consider a finite population of particles {Pi, ieI} satisfying the following two conditions: a) each particle is convex; b) pairs of particles can touch each other but triplets cannot. Assume that one can construct (for instance by erosion or segmentation of the population of particles) another population of particles {Qi, ieI} such that: c) each Qi is non empty, convex and contained within Pi; therefore the Qi’s can be referred to as subparticles; d) distinct subparticles are pairwise disjoint. From both populations, one defines the sets X = U Pi and Y = U Qi . The EPC is denoted by N. i i Let NC(Pi) be the coordination number of particle Pi: NC(Pi)=lN ( PiIPj ) . a) c) Fig. 2. Sequential analysis of particles: a) segmentation and selection of one particle, b) counting of the neighbouring particles, c) coordination number: NC =48/15=3.2. 56 Image Anal Stereol 2006;25:55-61 The coordination number of the population can be written ? NC(Pi) ?? NPinPj ) NC i j ?i ?n(pi) i ?N (Pi ) Because of b), the inclusion-exclusion formula for the EPC yields N ( X ) = ? N ( Pi ) - ?? N ( PinPj ) i j > i ?pj- 12?? N PinPj ) i j ? i On the other hand, because of d), we have ? n(pi ) = ? N(Qi ) = N(Y) from which it follows i i ?? NPinPj ) =2 [N(Y)-N(X)]. i j ? i Finally: NC 1 N( X ) N(Y) (1) This formula is illustrated on Fig. 3 with the two measurements on the initial connected set X and on the disconnected set Y. LOCAL MEASUREMENTS The two previous methods provide exact results only when the structure is globally explored. They must be adapted when the structure has to be analysed locally, which arises when the set of interest cannot be examined at one go. Sequential analysis In order to be applicable, this method requires the neighbourhood of each particle to be totally contained within the measurement field. Equivalently, any particle that has its neighbourhood hitting the edges of the field cannot be properly measured and therefore must be discarded (cf. Figs. 4 and 5). Moreover, the biggest particles have the largest propensity of hitting the edges of the frame. The «bias» incurred can be compensated by the Miles-Lantuéjoul correction (Miles, 1974; Lantuéjoul, 1980). It consists of assigning each particle a weight that is proportional to the chance it has to be contained within the measurement field. N(X) = 1 - 10 = - 9 N(Y) = 15 NC=2 [ 1-(-9)/15 ]=3.2 Fig. 3. Global measurement of the coordination number. Fig. 4. Local measurement by sequential analysis: the global value of the coordination number (4/3) is out of reach with this method. i i 57 Jouannot-Chesney P et al: Determination of the coordination number Fig. 5. Local measurements on the previous set of Fig. 2. The coordination number calculated by summing up the results from each measurement field is NC= 32/11 »2.91. Parallel analysis As will be seen in the following sections, the parallel analysis can be easily adapted to the local case. Local contribution to the Euler-Poincaré characteristic In a recent paper (Jernot et al., 2004), it has been established that the EPC of a bounded set X of Rd can be obtained by superimposing arbitrarily a tessellation on X (i.e., a partition of Rd into convex polyhedral cells), measuring the contribution of each cell to the EPC of X and summing all contributions. The local contribution of a cell Z to the EPC of X is a weighted sum of the EPC of X within each facet F§ of the tessellation contained in Z. The weights are chosen to compensate for the fact that facets belong to several cells. It is explicitly defined by the formula: ( - 1) d-dimF FcZ ordF and ordF stands for its multiplicity order, i.e., the number of cells sharing F as a facet. In the case of a 2D square tessellation of the space the orders of the facets are 4 (vertices), 2 (edges) and 1 (face). The measurement is illustrated in Fig. 6 on a 2D example: in the left field of this figure, the EPC values 1, 2 and 2 come respectively from the non-empty intersections of X with one vertex, two edges and one face. Local measurement of the coordination number for a finite family of particles In terms of local contributions, since N(X)= LC(X;Z), the coordination number can be expressed as: IC(X;Z) NC = 2|1- Z y C(Y;Z) Z (2) or, equivalently, as a ratio of expectations obtained by selecting one cell at random among all the cells hitting X: NC 2 t E { C(Y;Z) } (3) The set presented in Fig. 4 can be analysed using this method: the contributions are equal in both measurement fields, C(X;Z) = 1/2 and C(Y;Z) = 3/2, which gives NC = 4/3 by applying Eq. 2. In Fig. 7, the population of particles presented in Figs. 2 and 3 is analysed through four adjacent fields Z1, ... , Z4. Because this population is bounded, both local and global analyses lead exactly to the same result: NC = 2 in which dimF denotes the dimension of F (0, 1, 2, ..., d) 1- 5/4-11/4-5/4-15/4 11/4 + 17/4 + 11/4 + 21/4 3.2 C(X;Z) = 1/4 - 2/2 + 2 = 5/4 C(X;Z) = 0 - 5/2 + 3 = 1/2 Fig. 6. Examples of local contributions to the Euler-Poincaré characteristic of a set. § In this paper, as in Jernot et al. 2004, a facet of a tessellation is not necessarily a face of a polyhedral cell. A facet is defined as an intersection between cells. 58 Z Image Anal Stereol 2006;25:55-61 C(X;Z1) = - 5/4 C(Y;Z1) = 11/4 C(X;Z3) = - 5/4 C(Y;Z3) = 11/4 as C(X;Z2) = - 11/4 C(Y;Z2) = 17/4 C(X;Z4) = - 15/4 C(Y;Z4) = 21/4 Fig. 7. Local measurement of the coordination number for a finite structure. Local measurement of the coordination number for an infinite family of particles At this stage, only finite populations of deterministic particles have been considered. Now we may wonder whether the definition of the coordination number can be extended to random, stationary and locally finite populations of particles. Let X& (resp. Y& ) be the random set made up of the union of all particles (resp. all subparticles). To define the coordination number, a natural idea is to replace in Eq. 1 the connectivity numbers by their corresponding specific values. Namely we put NC 2 f1- N,(X nJY) (4) where N„(X) and N„(y) are the specific connectivity numbers of X and Y. These are defined as vol(rD) N*(X) lim r—»oo and N,(Y) = lim E{ N(YnrD) } vo l(rD) whatever the convex domain D with nonempty interior (Stoyan et al, 1995). In practice, Eq. 4 is not directly applicable because both specific connectivity numbers are unknown. Various estimators are available for N„(X) and N„(y) depending on the domain where the realizations of X and Y are known: - if the domain consists of two parallel section planes separated by a known distance, the disector technique can be used (Sterio, 1984). if the domain is a cell of a regular tessellation (i.e., all cells are identical up to a shift to a reference cell Z) the contribution technique can also be effective. It yields the intuitive but nontrivial fact that e{ c(x;z) } vol(Z) N*(X) and (5) N*(Y) e{ c(y;z) } vol(Z) the proof of which is left in a forthcoming paper. Starting from (5), unbiased estimators for N„(X) and N*(y) can be obtained, either by averaging the contributions of different cells or by assigning them different weights provided by a Kriging exercise (Chiles and Delfiner, 1999), which makes it possible to integrate the spatial structure of the random set under study. It should be pointed out that the knowledge of one realization X of X in two parallel section planes is insufficient to properly determine the corresponding realization Y of Y . A close look at Fig. 9 will easily convince the reader that the determination of Y requires X to be known in a full-dimensional domain. For the sake of consistency, it is deemed preferable to estimate both connectivity numbers also in a full-dimensional domain, thus resorting to the contribution technique rather than to the disector technique. From Eqs. 4 and 5, it comes NC E(C(X;Z e { c(Y;Z (6) This is basically the same formula as Eq. 3, except that the expectations have different interpre- r^co 2 59 Jouannot-Chesney P et al: Determination of the coordination number tations. They are taken over all cells hitting X in Eq. 3 and over all realizations of X in Eq. 6. As a consequence, suppose that the contributions of the realizations X and Y have been measured on the cells (Z j, j g J) , then we propose to use A NC 5>j CX;Zj ) 1- jeJ (7) as an estimate of the coordination number. In this formula, both sets of coefficients (A,j,jeJ) and (Hj, jeJ) are either constant weights (equal to 1/Card J) or Kriging weights. They add up to one. A PRACTICAL EXAMPLE The population of particles displayed in Fig. 8 is a simulation of a disc packing (Bideau, 1983). It has been chosen because its coordination number is explicitly known to be equal to 4 as a consequence of the underlying physical process (gravity) that governs its construction. On this population of particles both analyses, sequential and parallel, were performed. As can be seen in Fig. 8, both results (4.095 and 4.109) are quite similar even if 28 traces of particles have been discarded during the sequential analysis. It should be pointed out that sequential analysis gives valuable results only when many particles have their whole neighbourhood accessible. For example, examining the population of particles through the fields of Fig. 9 cannot yield any relevant result. In contrast to this, the parallel analysis can cope with this limitation. On each field the contributions for X and Y were determined. The sums of all contributions, 2C(X;Zj)=-58 and Xc(y;Zj)=55, are nothing j j else but the values of Fig. 8. Then, application of Eq. 7 using equal weights (lj = jLij = 1/Card J, Vje J) gives an estimated value of the coordination number equal to 4.109. On the other hand, it is also interesting to estimate the coordination number on each field (see Fig. 9). The average of these estimates is equal to 4.113 which is not very far from the theoretical value 4. This probably comes from the fact that all the particles have the same size and the distribution of the number of neighbours per particle is tight. This 2D example may look somewhat academic but it stresses the limitations of the sequential analysis. It turns out that these limitations are even more severe in 3D space. At the outset, there is a problem in the choice of the resolution that creates a conflict between the size of the particles and that of the measurement field. A high resolution leads to a good definition of contacts but to a small number of particles, whereas a low resolution allows a large number of particles to be analysed but the contacts are poorly defined. Nc = (38 x 4 + 4 x 5)/42 * 4.095 C(X;Z) = 4/4 - 24/2 - 47 = - 58 C(Y;Z) = 4/4 - 32/2 + 70 = 55 Nc=2 ( 1-(-5855) )*4.109 Fig. 8. Local estimation of the coordination number for an infinite structure: comparison between sequential and parallel analyses. 60 2 Image Anal Stereol 2006;25:55-61 C(X;Zj) C(Y;Zj) A NC -5.25 -7.00 -4.25 -7.00 -4.25 -6.50 -5.25 -7.50 -5.25 -5.75 5.75 5.50 5.25 5.50 5.25 6.00 4.75 6.00 6.25 4.75 3.826 4.545 3.619 4.545 3.619 4.167 4.211 4.500 3.680 4.421 Fig. 9. The sequential analysis is impossible when the population of particles of Fig. 8 is examined through ten narrow fields. In this case, only the parallel analysis is operatory. CONCLUSION In order to estimate the coordination number of particles in a powder stacking, a new method, based on the Euler-Poincaré characteristic, has been proposed. This method is independent from the workspace dimension. It is easy to implement and fast. It has been extended to the case of a local analysis using the local contribution to the Euler-Poincaré characteristic. This method is robust w.r.t. the nature of contacts (punctual or not) and the particle shape. ACKNOWLEDGEMENTS The authors want to express their deepest thanks to the referee for his constructive remarks and criticisms. REFERENCES Artz E (1982). The influence of an increasing particle coordination on the densification of spherical powders. Acta Metall 30:1883-90. Bideau D (1983). Relation entre les propriétés de transport de matériaux constitués d'empilements désordonnés isotropes et leurs caractéristiques géométriques. These d'Etat. Université de Rennes I, p. 23. Chiles JP, Delfiner P (1999). Geostatistics: Modeling spatial uncertainty. New-York: Wiley. German RM (1989). Particle packing characteristics. Princeton: MPIF. Guyon E, Troadec JP (1994). Du sac de billes au tas de sable. Paris: Odile Jacob. Hadwiger H (1957). Vorlesungen über Inhalt, Oberfläche und Isoperimetrie. Berlin: Springer Verlag. Jernot JP (1991). Pressureless sintering in solid phase: mechanisms, morphology and modelisation, Physics of granular media. New-York: Nova Science Publishers. Jernot JP, Jouannot-Chesney P, Lantuéjoul C (2004). Local contributions to the Euler-Poincaré characteristic of a set. J Microsc 215:40-9. Lantuéjoul C (1980). On the estimation of mean values in individual analysis of particles. Microscopica Acta 5:266-73. Lantuéjoul C, Beucher S (1981). On the use of the geodesic metric in image analysis. J Microsc 121:39-49. Miles RE (1974). On the elimination of edge effects in planar sampling, Stochastic geometry. In: Harding EF, Kendall DG, eds. New-York: Wiley. Ohser J, Nagel W (1996). The estimation of the Euler-Poincaré characteristic from observations on parallel sections. J Microsc 184:117-26. Serra J (1982). Image analysis and mathematical morphology. London: Academic Press. Sterio DG (1984). The unbiased estimation of numbers and sizes of arbitrary particles using the disector. J Microsc 134:127-36. Stoyan D, Kendall WS, Mecke J (1995). Stochastic geometry and its applications (2nd edition). Chichester: Wiley. Tancret F, Desgardin G, Osterstock F (1997). Influence of porosity on the mechanical properties of cold isostatically pressed and sintered YBa2Cu3O7-x superconductors. Phil Mag A75:505-23. 61