Image Anal Stereol 2007;26:129-136 Original Research Paper UNBIASED ESTIMATORS OF SPECIFIC CONNECTIVITY J ean -P aul J ernot , P atricia J ouannot and C hristian L antuejoul CRISMAT, Ensicaen/CNRS, 6 boulevard du mare´chal Juin, 14050 Caen, France; 2Ecole des Mines, 35 rue Saint-Honore´, 77305 Fontainebleau, France e-mail: jean-paul.jernot@ensicaen.fr, patricia.jouannot@ensicaen.fr, christian.lantuejoul@ensmp.fr (AcceptedNovember 7, 2007) ABSTRACT This paper deals with the estimation of the specific connectivity of a stationary random set in Rd. It turns out that the “natural” estimator is only asymptotically unbiased. The example of a boolean model of hypercubes illustrates the amplitude of the bias produced when the measurement field is relatively small with respect to the range of the random set. For that reason unbiased estimators are desired. Such an estimator can be found in the literature in the case where the measurement field is a right parallelotope. In this paper, this estimator is extended to apply to measurement fields of various shapes, and to possess a smaller variance. Finally an example from quantitative metallography (specific connectivity of a population of sintered bronze particles) is given. Keywords: Boolean model, edge effects, Euler-Poincare´ characteristic, specific connectivity, unbiased estimation. INTRODUCTION A UNIDIMENSIONAL EXAMPLE A few years ago, one of us was reviewing a paper in which an asymptotically unbiased estimator was used exactly as an unbiased one. As we made reservations about the way of handling this type of estimator, the authors reply “This estimator is UFAPP (unbiased for all practical purposes)”. Although humorous and witty, this reply should not hide the difficulty of using asymptotically unbiased estimators and interpreting the results they produce. This paper is mainly devoted to the estimation of the specific connectivity of a stationary random set. It turns out that the “natural” estimator is only asymptotically unbiased. The example of the boolean model illustrates the importance of the bias produced when the measurement field is relatively small with respect to the range of the random set. In the case where this measurement field is a right parallelotope, literature shows that this natural estimator can be made unbiased by adding a correcting term, e.g.,, Mecke et al. (1990) or Prasad et al. (1990). The unbiased estimator thus obtained can actually be generalized to apply to measurement fields of various shapes, and to possess a smaller variance. This generalisation is presented here only in two dimensions, but there is no doubt that it extends to spaces of any dimension. It should be pointed out that the bias studied in this paper originates from a support problem (limited size of the measurement field), not a discretisation one. The reader interested in the latter problem can consult Serra (1982) or Ohser et al. (2003). 129 Let X be a stationary, unidimensional boolean model of segments. It is shown in Appendix 1 that the mean connectivity ofX on [0,a] can be written as E{N(Xn[0,a])} = eae-em+1-e-em, (1) where 0 is the Poisson intensity and m is the mean segment length. Starting from Eq. 1, the specific connectivity of X is obtained first by dividing by a and then letting a tend to infinity NL = lim E{N(Xn[0, a ])} = Be -*™ . (2) Eqs. 1 and 2 show that N(Xn [0,a])/a is a biased estimator of Ni, although an asymptotically unbiased one. In order to quantify the importance of the bias, it is natural to consider the difference E{N(Xn[0,a])\ 1-e-Qm ------------------------Nl =-----------. Up to the factor a, the right hand side member is nothing but the probability that a point belongs to the boolean model (Matheron, 1975), or equivalently, the expected connectivity of X at one point. Accordingly NdÒ) = - [N ( Xn [0,a]) -N ( Xn {0})] , a N(a) = - [N(Xn [0,a]) -N(Xn {a})] , a a Jernot JP et al: Unbiased estimators of specific connectivity are unbiased estimators of Nl. It is actually better to consider their average, Nl = - N Xn [o,a]) - JNfXn {0}) - JN Xn{a})l , a 2 2 (3) which is also unbiased and has a smaller variance than each of them. Indeed, for x = 0 or x = a, we have Var{NL(x)} — Var{NL} 1 4a2 Var {loeX — laeX} > 0 SPACE OF ARBITRARY DIMENSION Let X be a stationary, ergodic random set in IRd. Suppose that a realization of X is available in a domain Z. Exactly as in one dimension, the natural estimator that is obtained by taking the connectivity of X n Z divided by the volume |Z| of Z is only asymptotically unbiased. In order to illustrate the importance of the bias incurred when the domain Z is too small, we have considered a d-dimensional stationary boolean model with intensity q and unit hypercubic objects. In the case where Z is a hypercube with edges parallel to those of the objects, the expected value of this natural estimator, or reduced mean connectivity (RMC), is equal to E{N(XnZ)} \Z\ ad -e-e d( J ai-dJlS(i,j)(-ey (4) In this expression, a is the edge length of the hypercubic field, and the S(i, j)’s are the Stirling numbers of the second kind (cf. Appendix 2). This gives a 1 a2 1 a3 -0 -0 - + 1 +02 a -0 O O O a a in one, two and three dimensions respectively. Fig. 1 shows this expected value as a function of the proportion of space (p(q) = 1 - e-q) occupied by the objects, for various field sizes ranging from 1 to ¥. The two-dimensional case is considered here in order to understand the shape of the curves produced. For small p, the objects are isolated and the RMC is nothing but the mean number of objects per unit area. As p increases, objects start to overlap and the RMC keeps on growing, but not so fast. Letting p growing again, the objects become numerous enough to create holes and the RMC starts decreasing. These holes contribute effectively to the RMC insofar as they are totally included within the field, which explains why their influence is more notable for large fields than for small ones. When p becomes close to 1, the holes become gradually filled by the arrival of new objects and the RMC starts again growing. Finally, when p = 1, the whole field is covered with objects, in which case the RMC is equal to 1/ad. a=1 a=2 a=4 a=16 a=64 — a=infinity 1 dimension ^—------------ —~^ — ^^ ^^^ ^^*^ 0.0 0.2 0.4 0.6 0.8 space proportion occupied by objects 1.0 2 dimensions 0.0 0.2 0.4 0.6 0.8 space proportion occupied by objects 1.0 a=1 3 dimensions a=2 — a=4 ^ — a=16 y^ — a=64 / — a=infinjty ^^ 0.0 0.2 0.4 0.6 0.8 space proportion occupied by objects 1.0 Fig. 1. Plotting the expected value of the “natural estimator” versus the proportion of space occupied by the objects in one, two and three dimensions, and for various f ield sizes. The specif ic connectivity curves are obtained for an inf inite size f ield. They are represented in black. 1 -- I 1 a 1 130 Image Anal Stereol 2007;26:129-136 Let us come back to our general problem. Unbiased estimators of the specific connectivity of X can be found in the literature when a realization of X is available in a right parallelotope Z, e.g.,, Mecke et al. (1990). They derive from the expression NV 1 Z\ [E{N(XnZ)}-E{N(XnuF3vF)}} , (5) where the union is taken over all faces of Z (vertices, edges. . .) containing the vertex v. In this equation, the mean connectivity of X in Z is amended by a term referred to as the “shell correction” in Prasad et al. (1990). Of course, Eq. 5 is valid whatever the vertex v. By averaging Eq. 5 over all vertices of Z, we obtain a new expression that can be written as NV 1 5j (-2) im E{N(XnF)} , (6) lZl FeF where F is the family of all faces of Z, including Z itself which is conventionally considered as a face of dimension d (cf. Appendix 3). From the estimation standpoint, Eqs. 5 and 6 yield the unbiased estimators Ny(v) NV [N(XnZ)-N(Xn(uF3vF))} 1 — ^ (-2)im ~ N(XnF) (7) (8) By construction, Ny and the N V/ (v) ’s are related by the expression NV 1 2d åNV(v and if the Ny (v) ’s are identically distributed, then the distribution of Ny is less dispersed than that of the NV(v)’s (cf. Appendix 4). More generally, Eq. 6 can be written as NV 1 Z\ X (-1)d-dimFaFE{N(XnF)} (9) FeF with aF = 2dimF-d. It turns out that this type of equation is valid for a wide range of polytopes. In what follows, we are going to make precise the role played by the coefficients aF in two dimensions, in order to facilitate their assessment. BIDIMENSIONAL SPACE Right triangle Suppose at first that Z = Tabc is a right triangle with vertices a, b and c, as in Fig. 2. Let Tbcd the triangle deduced from Tabc by symmetry with respect to the midpoint of [b,c\. The union of both triangles R = Tabc U Tabd is a rectangle on which Eq. 6 applies. d a b Fig. 2. Case of a right triangle. In order to get more concise expressions, we put cf = E{N(X H F)} from now onwards. Because cab = ccd, cac = cbd as well as ca = cb = cc = cd, we can write NA 1 \R\ \cR-cab-cac + ca] But we also have cr = cTabc + cTbcd - cbc. Moreover cTabc = cTbcd if the distribution of X is symmetric w.r.t. the origin. Under this latter assumption, we obtain NA 1 Tabc cT abc 11 (cab +cac +cbc)+ ca 22 (10) Scalene triangle More generally, let us take for Z a scalene triangle Tabc. Let d be the projection of a onto the line that supports b and c. Eq. 10 applies on both right triangles Tabd and Tacd: 2TabdNA = 2cabd - cab - cad- cbd + ca , (11) 2TacdNA = 2cacd - cac - cad- ccd + ca ¦ (12) At this stage, two cases must be distinguished, depending on whether d lies inside or outside [b,c] (cf. Fig. 3): ^_ Fig. 3. Case of a scalene triangle. - if d € [b,c], then Tabc is the union of Tabd and Tacd. By adding Eqs. 11 and 12, and owing to c abd + c acd = cabc + cad as well as cbd + ccd = bc + a, we arrive at Eq. 10. _ c _ _ _ _ _ _ v a a _ b d b d c c 131 Jernot JP et al: Unbiased estimators of specific connectivity – if d ? [b,c], then Tabc is the set difference between Tabd and Tacd. By subtracting Eqs. 11 and 12, and owing to cabd -cacd = cabc -cac as well as cbd -ccd =cbc -ca, we end up again with Eq. 10. It should be pointed out that formula (Eq. 10) is not fully satisfactory because it does not assign the same role to the vertices a, b and c. Yet, we can resort to the fact that the angles of a triangle add up to p and write 1 ca =aaca +abcb +accc , 2 where av denotes the internal angle of vertex v in the triangle divided by 2p. This interpretation is interesting because it is applicable to all faces of the triangle and to the triangle itself. The factor 1/2 associated with an edge can be written as p/(2p), and pis precisely the internal angle at any inner point of that edge (cf. Fig. 4). Similarly, the coefficient 1 associated with the triangle can be written as (2p)/(2p), and 2p is precisely the internal angle at any inner point of the triangle. Fig. 4. Internal angles in a triangle. The function a can be therefore defined on all faces of Tabc including Tabc itself. aF is called the internal angular proportion of face F. Using this function, formula (Eq. 9) is valid for triangles. Polygon Let Z be a simple polygon (polygon limited by a simple Jordan curve). How to estimate NA starting from measurements in Z? Let (Ti,i ? I) be a triangulation of Z. This means that the Ti’s have disjoint interiors and that their union is exactly Z (cf. Fig. 5). Fig. 5. Triangulation of a polygon. Resorting to a triangulation is advantageous because formula (Eq. 9) can be applied to each element of the triangulation. The results obtained can then be combined via the inclusion-exclusion formula. This finally gives NA 1 å (-1)dimFaFcF |Z| F?F (13) which once again involves the internal angular proportion of all faces of Z (cf. Fig. 6). Details of the proof are reported in Appendix 5. Fig. 6. Internal angles in a polygon. It should be pointed out that this approach is applicable whenever the field Z can be triangulated. Accordingly, formula (Eq. 13) is valid not only for simple polygons, but also for polygonal domains that will be examined below. There is another way of writing formula (Eq. 13). Let e1,...,en be the edges of Z, and let v be an arbitrary point (vertex or not). As each edge has an internal angular proportion of 1/2, and because the internal angular proportions of each vertex add up to (n 2)p/2p= (n-2)/2, this formula becomes NA 1 Z\ 1 n (n-2) cZ - åcei +cv 22 i=1 Note also that the inclusion-exclusion formula relates the sum of all edge contributions to the mean connectivity associated with the boundary ¶Z of Z n åcei i=1 c¶Z +ncv Replacing the edge contributions by their expression, a second formula is obtained for NA XZ- 2XdZ-Xv 1 / 1 NA= ZAXZ- 2XdZ - Xv (14) Finite union of pairwise disjoint polygons Here Z is the union of a finite family (Zi,i ? I) of pairwise disjoint polygons. Formula (Eq. 14) applies to each polygon \Zi \Na = XZi -XdZi - Xv i (ž I But åicZi =cZ and åic¶Zi =c¶Z. By summation over all polygons and division by |Z|, one obtains NA 1 Z\ cZ 1xdZ-\I\Xv (15) _ a b c _ _ 132 Image Anal Stereol 2007;26:129-136 Polygonal domain Here Z is a domain limited by a union of simple polygonal curves. Suppose that Z has |I| connected components and |J| holes. Note that cZ and c¶Z are not affected by a shift of the connected components of Z as long as they remain mutually pairwise disjoint. Accordingly, it can be assumed without loss of generality that Z has no connected component located within one of its holes, as in Fig. 7. Let H be the closure of the union of the holes. In what follows, H is supposed to be contained within the interior of Z ?H. Fig. 7. Example of polygonal domain. Let us apply formula (Eq. 15) to Z U H and to H \ZUH\Na = XZUH----Xd(ZUH) ~ \I\Xv > \H\Na = Xh — 2XdH ~ \J\Xv . By difference, this gives 1 \z\na = Xzuh-Xh - z {Xd(ZuH) - XdH ~ \I\ ~ \JXv Now we have (16) Xzuh = XZ + Xh - XdH because ZC\H = dH. On the other hand, the boundary of Z can be split into its “inner” boundary (produced by the holes of Z) and its “outer” boundary (obtained by filling the holes of Z). Accordingly, we have d Z = dH Ud (Z UH), from which we immediately derive Xd (ZuH) = XdZ-XdH . Replacing Xzuh and XdlZuH) by their expression in Eq. 16, we obtain \Z\Na = XZ 2XdZ (\I\ - \J)xv and it remains to observe that |I| - |J| is nothing but the Euler-Poincare´ N(Z) of Z to conclude NA \Z [XZ c¶Z -N(Z)cv 2 (17) AN EXAMPLE Classical stereological relationships are often used to characterise microstructures in many scientific fields (biology, geology, materials sciences etc...). The practical application of these relationships rests on the measurement of the specific connectivity. In order to show the difference between its natural estimates and its unbiased ones, as provided by Eqs. 6 and 13, measurements have been carried out on polished sections through a sintered material (Jernot, 1982). The microstructures presented in Fig. 8 correspond to a 60 mm bronze powder, sintered under argon at 730 ?C, 750 ?C and 770 ?C during 2 hrs. The volume fraction of the solid phase is respectively equal to 0.657, 0.716 and 0.790. Fig. 8. Three polished sections of bronze particles sintered during 2 hrs at three dif ferent temperatures. 1 _ 133 Jernot JP et al: Unbiased estimators of specific connectivity The beginning of the sintering is associated to the growth of the necks between the particles and followed by a progressive closing of the porous network due to the densification of the material. This coalescence of particles is reflected in the microstructure by a simultaneous modification of the area and curvature of the solid/pore interface and, consequently, of its specific connectivity. A quick examination of the results obtained from these images shows that the discrepancies between the natural and the unbiased estimates are far from negligible at every stage of sintering. This is particularly manifest on the two-dimensional results. Matheron G (1975). Random sets and integral geometry. New York: Wiley. Mecke J, Schneider RG, Stoyan D, Weil WRR (1990). Stochastiche Geometrie. Basel: Birkha¨user. Ohser J, Nagel W, Schladitz K (2003). The Euler number of discretised sets – surprising results in three dimensions. Image Anal Stereol 22:11–9. Prasad BP, Lantue´joul C, Jernot JP (1990). Use of the shell correction for quantification of three-dimensional images. Trans Roy Microsc Soc 161-3:387–403. Serra J (1982). Image analysis and mathematical morphology. London: Academic Press. 730 °C 750 °C 770 °C 1Dnat. (mm 1) 15.17 11.92 9.35 1Dunb. (mm-1) 13.88 10.37 7.97 2D nat. (mm-2) 74.31 -83.60 -46.44 2D unb. (mm-2) 11.61 -125.40 -84.76 CONCLUSION The connectivity of a stationary random set X in a field of measurement, divided by the d-volume of the field, is only an asymptotically unbiased estimator of the specific connectivity of X. Explicit calculations on the boolean model show that the bias can be important when the field is relatively small with respect to the size of the objects. In two dimensions, two unbiased estimators have been proposed when the distribution of X is symmetric w.r.t. the origin, and when the field of measurement is a polygonal domain. The first one (Eq. 13) requires the connectivity of X in the whole field, along each edge and at each vertex. It involves the internal angular proportion of each face. The second one (Eq. 17) requires the connectivity of X in the whole field, along its boundary and at a point. It involves the Euler-Poincare´ characteristic of the field. ACKNOWLEDGEMENTS The authors want to thank Dr W. Nagel for his careful reading of the paper and his judicious remarks. REFERENCES Aigner M (1997). Combinatorial theory. Berlin: Springer. Jernot JP (1982). Analyse morphologique et mode´lisation du frittage et des mate´riaux fritte´s. PhD Thesis, Caen University. Lantue´joul C (2002). Geostatistical simulation. Models and algorithms. Berlin: Springer. APPENDICES APPENDIX 1: PROOF OF EQ. 1 Let X be a unidimensional, stationary boolean model of segments with Poisson intensity 0 and length cumulative distribution function F. Let also a > 0. The aim of this appendix is to calculate the mean connectivity E {N(X n [0, a])} of X over [0, a]. Let us assume here that F has a probability density function f as well as a finite mean m. Following Lantue´joul (2002), X C\ [0,a] is a boolean model, but its parameters are not those of X. More precisely, Xn [0,a] is distributed as the union in [0, a] of a Poisson number K of independent “typical” segments. The mean of K is û = Q(a + m). By typical segments, we mean segments [Xk ,Xk+Lk ] with respective probability density functions Xk ~ h(x) 1-F(-xA0) 1x 0. To do so, observe first that p„=E{N(n1<„[X1,X1+L1])}. Now, r\i<„[Xi,Xj + Lj] = 0/ if and only if Xm = \liof=i a + m The calculation of I\ is straightforward. We obtain with Pnï=lA, ^ 0 or equivalently \ à id d P{n^A^0} = (lj (l ij , by putting l = m/a. Plugging this expression into Eq. 22, we obtain after some calculations E \N(X n Z)} = 2_, (_l)»-i 1 n\ d / d\ nil-i . (23) Il m n a+m (19) To go further, it is convenient to introduce the Stirling numbers of the second kind. They can be defined by their generating function (Aigner, 1997) To calculate I2, we write r n I=n\ i TT 1 \ 1 —F(xn —xi) .=1 a + m dx\ ¦ ¦ ¦ dxn ni = å S(i, j)(n)j j=0 (24) the integration domain being xi < x2 < • • • < x„ and 0 < xn < a. Using the change of variables y i = xn — x, for i = 1,...,n— 1 and y„ = x„, the n integrals separate and we get i m \ a a + m with the usual convention (n)^¦ = n (n — 1 ) • • • (n — j + 1 ) if j > 0 and (n)o = 1. The following table gives the first Stirling numbers I2 = n a+m a+m Combining Eqs. 19 and 20, we get mn-1 (20) j = 0 j = 1 j = 2 j = 3 i = 0 1 i = 1 0 1 i = 2 0 1 1 i = 3 0 1 3 1 pn (a+m)n (m+na) (21) It remains to plug the expression of pn into Eq. 18. This finally gives We now replace nl by its generating function in Eq. 23, and because 00 Jn mn-^ E{N(Xn[o,a])} = y(-i)"- —------- (m+na) n{aTmy that is E{N(Xn[0,a})} = qae-dm + l-e-dm , after summation. APPENDIX 2: PROOF OF EQ. 4 Actually, we are going to prove a little more by considering a boolean model whose objects are right parallelotopes with independent and identically when j > 0, whereas n=\ n! (qmd)"(n)o = l-e-emd it follows E{N(XnZ)\ = l-e-emd Y (d)l-iS(i,j)(-qmdy , which gives Eq. 4 after dividing by |Z| = ad. Note that one gets the specific connectivity by letting a tend to 00: Nv = -e-e»>dLy d ¦ 135 Jernot JP et al: Unbiased estimators of specific connectivity APPENDIX 3: PROOF OF EQ. 6 Here we put M(v)=N(Xn(uF3vXnF)) . Let {Fi,i e I(v)} the set of all faces of Z that contain v. Because UF3vF = UieItv\Fi, and owing to the distributivity of the intersection over the union, we have M{v)=N{\JieI{v)XC\Fi . By the inclusion-exclusion formula, this becomes M(v)= £ (-1 J -NÇXnFJ) , 0^JCI(v) where FJ = djeJFj and |J| is the number of elements of J. Now it can be noticed that for each face F of Z that contains v there exists one J and only one such that F = FJ. Moreover |J| =d - dimF. Accordingly M(v) = Jj(-1)d-1-d™FN{XnF). F3v We now take the mean over all vertices of Z and use the fact that each face of Z has 2dimF vertices. Then we arrive at 1 å 1 å M(v) = d 2dimF( 2d 1 d- 1 - dimF N(XnF) -£(-2)dimF-dN Xf]F Various inequalities can be obtained depending on the choice of j. For instance, j(t) = t gives E{NV} < E{Ny(v)}. Similarly, j(t) = -t gives - E{Ny} < -E{Ny(v)}. Therefore Ny and Ny(v) have the same -—-2 mean. If we also take j(t) = t , we obtain E{NV } < -—~ 2 E{Ny(v) }. In combination with the equality of the means, this gives Var{Ny} < Var{Ny(v)}. APPENDIX 5: PROOF OF EQ. 13 It is always possible to choose the triangulation (Ti, i G I) in such a way that the vertices of its triangles are vertices of Z. The standard inclusion-exclusion formula gives c Z =E{N(XriUieITi)} = X (- 1)\j\-1cTj , with TJ = C\j(zJTj. The TJ are either triangles (for all J with |J| = 1), or internal edges (edges of the triangles that are not edges of Z, for certain J’s with |J| = 2), or even vertices of Z (for all others J such that |J| = 2 and all J satisfying |J| > 3). It turns out that the contribution of each vertex v of Z vanishes. Indeed, assume that v belongs to n triangles. Then v is incident to n - 1 internal edges. Its contribution is therefore It remains to take the expectation and plug this equality into Eq. 5 that has been averaged over all vertices of Z to get Eq. 6 as desired. APPENDIX 4: STATISTICAL COMPARISON OF VARIOUS ESTIMATORS Suppose that all Ny (v) ’s are identically distributed. Let v be a vertex selected at random. Then Ny(v) shares the same distribution. Moreover E{NV(v) | NV} 1 2d ^E{NHv) | Ny} v = E{Ny | Ny} — N V/ Using Jensen inequality, this implies that we have j(Ny) = j(E{NyXv) | Ny})