Image Anal Stereol 2010;29:99-110 Original Research Paper STEREOLOGICAL ESTIMATION OF SURFACE AREA FROM DIGITAL IMAGES Johanna Ziegel^ and Markus Kiderlen^ ^Department of Mathematics and Statistics, The University of Melbourne, Victoria 3010, Australia; ^Department of Mathematical Sciences, University of Aarhus, Ny Munkegade Build. 1530, DK-8000 Aarhus C, Denmark e-mail: jziegel@unimelb.edu.au, kiderlen@imfau.dk (Accepted May 12, 2010) ABSTRACT A sampling design of local stereology is combined with a method from digital stereology to yield a novel estimator of surface area based on counts of configurations observed in a digitization of an isotropic 2-dimensional slice with thickness s. As a tool, a resuh of the second author and J. Rataj on infinitesimal increase of volumes of morphological transforms is refined and used. The proposed surface area estimator is asymptotically unbiased in the case of sets contained in the baU centred at the origin with radius s and in the case of baUs centred at the origin with unknown radius. For general shapes bounds for the asymptotic expected relative worst case error are given. A simulation example is discussed for surface area estimation based on 2 X 2 X 2-configurations. Keywords: configurations, digital stereology, local stereology, surface area. INTRODUCTION Algorithms that determine the surface area of a three-dimensional object from a binary image can be grouped into global and local methods (Klette and Rosenfeld, 2004, p. 302). Global methods usually require a digital approximation of the surface or calculation of the surface normals. Local methods estimate the contributions to surface area of small voxel blocks in the binary image independently of the image outside the given block, and add up these contributions to estimate the total surface area. As local methods can be implemented using a linear filter, they are extremely efiicient. Lindblad (2005) suggested a local estimator based on voxel blocks of size 2x2x2 inspired by the common marching cubes algorithm, but without an explicit reconstruction of the surface. There are 2^ = 256 possible combinations of foreground and background voxels in a 2 x 2 x 2 block. They will be called configurations in the following. Two of these configurations contain only background or only foreground voxels and do therefore not contribute to the surface area. The remaining configurations, called boundary configurations, can be grouped into 14 classes, due to symmetry. Thus, local estimators of the surface of a set X in three-dimensional space are of the form 254 i=i (1) where Ay is the contribution and Ni is the total number of occurrences of configuration i. Several natural choices of the weights Ai,...,A254 can be found in the literature (Lindblad and Nyström, 2002; Lindblad, 2005; Schladitz et al, 2006; Ziegel and Kiderlen, 2010). Assuming that the set X was randomly translated before digitizing, Ziegel and Kiderlen (2010) found optimal weights in the sense that the asymptotic average error is minimized among all estimators of the form Eq. 1; see also Gutkowski et al. (2004) for a weaker result. Both papers are based on an asymptotic formula (Kiderlen and Rataj, 2006), where asymptotic refers to increasing resolution of the image. Ziegel and Kiderlen (2010) have shown that only 102 of the 254 configuration types can arise in the case where X is bounded by a planar surface, and it is only these that contribute to the surface area measure asymptotically; see also (Lindblad, 2005). These 102 configurations are therefore called i/j/ormatfve configurations. The estimator Eq. 1 can in principle also be applied to stacks of binary images of horizontal planar sections of X, where configurations are composed of voxels in two neighbouring section planes. A possible application is the estimation of the surface area of particles (e.g. cells) in confocal microscopy, where a stack of focal planes is digitized and analysed. However, due to optical effects, section planes that are hitting a particle close to one of the two horizontal touching planes yield only an extremely blurred image of the section profile. In traditional (continuous) stereology, this problem is solved by considering a random isotropic slice centred at a reference point of the particle under consideration. If the thickness 5 > 0 of the slice is chosen small enough compared to the size of the particle, the blurring effect will not be expressed within the slice. An unbiased surface area estimator of Horvitz-Thompson type is then obtained by weighting the surface area contribution of each infinitesimal surface patch in the slice with its inverse sampling probability. This sampling probability is a Sanction of the distance of the patch from the reference point, and depends on the thickness 2s of the slice (Jensen, 1998). It is the purpose of this paper to show that the two concepts from digital and local stereology described above can be combined. We will obtain a surface area estimator based on a stack of planar parallel digital images in an isotropic local slice, where the planes of the batch are parallel to the slice. It is intuitively clear that an estimator will no longer be a weighted sum of configuration counts like in Eq. 1. Instead, each observed configuration must additionally be weighted according to its inverse sampling probability. However, one cannot expect to obtain an unbiased estimator this way. We will propose a set of weight fianctions such that the surface area estimator is asymptotically unbiased in the cases where X is a ball centred at the reference point or X is contained in such a ball of radius 5; see Corollaries 8 and 9. In all other cases we can determine explicit bounds for the asymptotic worst case error; see Ineq. 11 and Proposition 10. In the next section basic notations and concepts are introduced together with a slight generalization of the asymptotic formula in Kiderlen and Rataj (2006) to weighted volumes of morphological transforms. This is the basis for the main theoretical result. Theorem 3, which describes the asymptotic mean behaviour of an estimator based on weighted configuration counts in an isotropic slice. In the penultimate section. Theorem 3 is used to establish a surface area estimator based on weighted counts of m different configurations in a digitization of an isotropic slice section of X. We determine estimates for the asymptotic relative mean error, and show that these can be improved, if X is known to be contained in a ball centred at the reference point with known radius; see Proposition 10. Up to this point, the results hold for digitizations on general lattices (with possibly different resolutions along the different axes, and with not necessarily orthogonal axes) and for voxel blocks, which may be larger than 2x2x2 configurations. In the final section we specialize these results to the scaled standard lattice L = tl?, define an estimator based on the 102 informative 2 x 2 x 2-configurations, and compare its performance in a simulation example with the theoretical asymptotic results. PRELIMINARIES By we denote the unit sphere in d- dimensional real coordinate space R'^. The standard scalar product on R'^ is (•,•) with associated norm II"II. By a l-subspace we mean a ^-dimensional linear subspace of R'^. Let A,B c. R''.^ The reflection of A at the origin is denoted by Ä = {—x \ x e A}, its complement by A'^ = and its topological boundary by d A. We write A(BB = {a+ b\ae A,b ^ E} for the Minkowski sum of A and B, and AeB= {xeR'' I x+ßc A} for the dilatate of A by B. The positive part of a real valued Sanction / is denoted by f^ = max( f, 0). The support fianction of a convex body K in R is denoted by h{K, •). We use this notion also for compact sets A, A ^ 9, defining h{A,-) = h{com{A),-), where conv(4) is the convex hull of A. The exoskeleton exo(4) of a closed set A is the set of all z G A'^, which do not have a unique nearest point in A. The set exo(4) is measurable and has Lebesgue measure zero (Fremlin, 1997). A closed set A' c R"' is gentle if for -almost all X G dX there are two non-degenerate open balls touching in x such that one of them is contained in ^ and the other in and if also Jf^-^NidXn 5 X 5^-1)) < - for afl bounded Borel sets B c R'^. Here ^^ is the ^-dimensional Hausdorff measure in R'^ and N{A) is the reduced normal bundle of A (Kiderlen and Rataj, 2006). The class of gentle sets is rather large. It contains for instance all convex bodies (compact convex subsets of R'^) with interior points, all topologically regular sets in the convex ring (the family of finite unions of convex bodies), and certain unions of sets of positive reach. At almost all boundary points a of a gentle set X there exists a unique outer unit normal n(a) to X. Let •) be the image measure of ondX under the map a ^ {a,n{a)). The measure •) vanishes outside N{X). Let ^dx R'^XexoidX) dX denote the metric projection. The following theorem is a generalization of (Kiderlen and Rataj, 2006, Theorem 1). Theorem 1. Let X (IW^ be a closed gentle set, f : R'^ ^ R a compactly supported bounded measurable function and B, W and P, Q four non-empty compact subsets ofR'^. Then lim ^ f f{^dx{x))dx 6^0+ e J[{xmP)esB\\i{xesQ)mw\ lN(X) /(a)(A(P© Q,n)-hiB® W,n))+ xQ_,iX;dia,n)). (2) If f is in addition continuous at all points in dX, then ca/3 be replaced by f{x) in Eq. 2. The proof of Theorem 1 can be found in the appendix. Kiderlen and Rataj (2006) show Theorem 1 in the case where / is an indicator function. Their result can be generalized to Theorem 1 essentially by applying the monotone convergence theorem. Let ,..., Xrf be a basis of R'^ and let L = {/31^1 H-----h /JrfXrf I /31,..., G Z} be the lattice generated by this basis. A given lattice L is generated by infinitely many different bases, but the volume of the fiandamental cell Q = [0, xi] © • • • © depends only on L and not on the basis chosen; see for example (Yap, 2000). This number is denoted by det(L). If is a uniform random variable in Co, then the random lattice + L is a stationary random lattice. The distribution of +L does not depend on the choice of Cq. We now define the digitization of a set A' c M'^. The points of + L are interpreted as voxel midpoints of a digital image, where each voxel is a translate of Q. Often, L is the orthogonal standard lattice Z^ and the voxels are small cubes. We work essentially with the Gauss digitization model of X consisting of all voxels having their midpoints in X (Klette and Rosenfeld, 2004, p. 56). As there is a one-to-one correspondence between the set of voxels on the one hand and + L on the other hand, the Gauss digitization is determined by the set Xn + L) of all lattice points in X. In order to vary the resolution of the digitization, the random lattice is often scaled by a variable t>Q and Xf] is called the digitization of X (with resolution 1 /t). In the following we only consider compact gentle sets X. Let the Sanction /be measurable non-negative or integrable. Let + L be a stationary random lattice, and let C L be two non-empty finite subsets of L. The points in B represent 'black' pixels of the digitization, i.e. points which are contained in the set X, whereas W stands for the 'white' points of the background. For t>Q define For / = 1 the random variable Nt counts the number of all translations of the pattern {tB,tW) in the digitization Xf^ + L). Calculation shows that nNt] = t-d det(L) J[xetB]\[x®tW\ f{x)dx. (3) continuous on dX and let B and W be two non-empty ßnite subsets of a lattice L. Then hm ?'^-idet(L)E [Nt] fia)i-hiB®W,n)yQ_,iX;dia,n)). Corollary 2. Let X C R'^ be a compact gentle set. Let f be a locally bounded measurable function, which is m Proof Let C be a compact set such that [Xq tW\ d C iox aW t> Q smaller than some fixed /o > 0 and X (Z C. Replacing /, P, B, Q and W in Theorem 1 by Ic/, {0}, B, {0} and respectively, yields the claim. □ COMBINING LOCAL AND DIGITAL STEREOLOGY In the following we restrict ourselves to R^. The results can be generalized to R'^, c/ > 4, in a straightforward manner. We prefer to present them only in R^ in order to keep the notation concise. Denote the standard basis vectors in R^ by 61,62,63. Let i? be a random proper rotation with distribution given by the normalized Haar measure on the rotation group SO^. Fix the 2-subspace 1q = span(61,62). We define the random 2-subspace L = RIq . It is uniformly distributed in the set J^ of all 2-subspaces of R^. Let ji be the distribution of L. For 7 G ^ and 5 > 0 define Ts = Z{1) = 7© 5(0,5). The set 7; = Z{L) = R{h © -S(0,5)) = I © B{0,5) is called a random 2-slice with thickness 2s. It will be clear from the context whether Tg refers to the deterministic 2-slice Ts{l) or the random 2-slice Ts{L). The following Theorem 3, which is the main theoretical result of the paper, gives a formula for the asymptotic mean of the weighted number of occurrences of two sets B and W (black and white points, respectively), which are specified below. Theorem 3 is an analogous result in a local stereological setup to Kiderlen and Jensen (2003, Theorem 4) in the context of stationary random sets in the plane, and to Gutkowski et al. (2004, Theorem 1) and Kiderlen and Rataj (2006, Theorem 5) in the setting of spatial objects, which are digitized by a stationary random lattice. Theorem 3. Let A' C R^ be a compact gentle set Let R be a random proper rotation and let ^ +1, be a stationary random lattice, which is independent of R. Let B,W d be two non-empty finite subsets of the lattice L and f a continuous non-negative function on R^. The weighted sum Nf of occurrences of{B, IV) in the digitization of X, which he entirely in Tg, is given by xetRi^+V) x+tR{BjW)0. It satisßes : / / f{a)H{(pl)Jf^{daMdl), J J^ J Ts where H{(p) is given by 1 r 4;r|cos^| Js^n{\(e3,-}\=\sm(j>\} (4) Hm ?2(detL)E[Af] =E[F^(i?)], where Fgir) := f g{a,rlo)h{r-'n)C2{Xn Z;d{a,n)) JN{XnTs) forrGSOi withh={-h{B(BW,-)y. (5) The proof of Proposition 4 can be found in the appendix. It uses that given the rotation R, Corollary 2 can be applied to Af. In the following two Lemmas we derive a formula for the expected value of Ej{R). Lemma 5. Let ^: R^ x ^ ^ R öe a non-negative measurable function such that g{-,l) is continuous for all 1 ^ J^. Let Fg{r) be given by Eq. 5. Then the conditional expectation E[Fg{R)\L = 1] can be expressed as idXnz and is the angle between 1 and the outer normal n{a) ofXataedX. The proof of Theorem 3 is based on the following Proposition 4 and Lemmas 5, 6 and can be found in the appendix. The difference of Nf in Theorem 3 and Nf in the following Proposition 4 is that the latter also counts configurations which do not lie entirely in the slice T2. These are of course not observable in practice. Proposition 4. Let A' C R^ be a compact gentle set. Let R be a random proper rotation and let ^ +1, be a stationary random lattice, which is independent of R. Let C L be two non-empty ßnite subsets of L. Let ^: R^ X ^ ^ R öe a non-negative bounded measurable function such that g{-,l) is continuous for all 1 ^ J^. Then the number M = X {x+tRB'z{Xr\Ts)r\tR{^+l.),\-> xetR{^+l.) \x+tRW 0 satisßes J XndTs for ß -almost all 1G where is given by Eq. 4. Proof For jU-almost all 7 G ^ we have that ^^{dXf] dTg) = 0. This implies that there is a unique outer unit normal n{a) for ^^-almost all a G d{Xri Tg). Hence we obtain for r G SOi Fg{r) = f g{a,rh)h{r-'n{a)).^\da). ■J d [Xn Ts) There exists a regular version of the conditional distribution of R given L (Klenke, 2006, Theorem 8.36). Therefore we can use Fubini's theorem to obtain nFg{R)\L=i] = J tdixnTs)ia)g{a,l)E[h{r-'n{a))\L = l]J^-\da). Recall that /x-almost surely la(A'nrs) = ^dXnTs + txndTs- Thß claim now follows from Lemma 6. □ Lemma 6. Forn^ andle^ we have nh{ir'n)\L = l]=H{^), where (j) is the angle between n and 1, and H{(p) is given by Eq. 4. Proof Fix p G SOi, such that pi = k. Then E[h{R-^n)\L = 1]= E[h{R-^n)\pRlo = /0] = E[hiir^pn)\Rlo = lo]. We have (/3,7) = {pn,lo), hence the last conditional expectation in the above equation can be written as the normalized integral over the two small circles C S^ parallel to k at height ± sin (j) with radius cos (j), where (j) is the angle between n and 1. □ AN ESTIMATOR FOR SURFACE AREA In this section we will derive an estimator for surface area using Theorem 3, which is based on a local stereological sampling design. Let A' G R^ be a compact gentle set. Suppose we observe ATl R[{h + ^(0,5)) n + L)] for some random proper rotation R, ^+L, a stationary random lattice, which is independent of R, and s,t > 0. Let {Bi,Wi}, i = \,...,whe boundary configurations ofL, i.e. Bj, are non-empty, disjoint finite subsets of L with Bi\jWi = Q)f^ L, where Co is a fixed fundamental cell of L. We define the following estimator for surface area 5(Z) = ?2(detL)XA?', (6) i where N^ = as defined in Theorem 3 with B= Bj, W= Wj. The continuous Sanctions Ay: [0,0°) ^ [0,oo) have to be suitably chosen according to the choice of {Bj, Wj). We give an example for 2 x 2 x 2-configurations in the concluding section. Proposition 7. We have J oX y (7) where Xf/a ^ [0, is the angle between a and the outer normal n{a) of X at a G dX. The function gi{r, \j/) for \j/ G [0, :7r] and r G [0,0°) is given by , , J/o*/i(arcsin(z))dz, forr s. With the two sided cut-off function x ^ = min{l,max{—l,x}}, the function for y/ G (0, n) and q G (0,1] can be written as = ^ (arccos -arccos , where aY,q{z) = {q - zcosi/A)/(sini/AVl -z^). For \j/ G {0,:7r} we have Go,q{z) = Gn,g{z) = q]{z) . Proof By Theorem 3 we obtain that lim^E[S{X)] =! ! , (8) J ^ J oXr\ Ts j where Hj is given by Eq. 4 for (5, W) = {Bi, Wi). Using Fubini's theorem we can rewrite the right-hand side of Eq. 8 as J oX j J ^ We have = giUl^t), which is a tedious but not difficult calculation. □ Note that for all i/a g [0, :/r] and all r G (5, <>=) = ^ (9) (Jensen, 1998). We assume from now on that none of the functions Z^-, i = 1,..., m, is identically zero. This is fialfilled, if and only if Bj and Wj can be strictly separated by a hyperplane for all i = I,...,m, or, in other words, Wj) is an informative configuration. The following two corollaries to Proposition 7 show that the weight Sanctions Ay can be chosen to yield an asymptotically unbiased estimator S{X) in the case where X is sufficiently small compared to 5 or when A' is a ball centred at the origin with unknown radius. Corollary 8. Let X (Z B{Q, s). SetXi{r) = aigi{Q,Q)-^ with coefficients ai,..., a^ G M that are summing up to one. Then the estimator S{X) as given in Eq. 6 is an asymptotically unbiased estimator of the surface area S{X). Proof By definition Hj>0. The assumption that Hi ^ 0 yields gi{0,Qi) > 0. As gi{r,\if) = gj(0,0) for all (r, y) g [0,5] X [0, n], Eq. 7 implies the claim. □ Corollary 9. Let X be a ball centred at the origin with unknown radius. Suppose that the sets Bi, Wi are such that gi{r,Q) > 0 for all r G [0,oo) and i = I,...,m. Choosing Aj(r) = aigi{r,Q)^^, where Ei=i ay = 1, yields an asymptotically unbiased estimator S{X) of S{X). Proof If A'is a ball centred at the origin we have = 0 for all a G dX. The claim follows from Eq. 7. □ The condition gi{r,0) > 0 in Corollary 9 is equivalent to requiring that the support of Hj contains an interval [0, e) for some e > 0. For general shapes we cannot expect to obtain an unbiased estimator of the form Eq. 7. A suitable choice of Xi for r > 5 will strongly depend on the choice of the pairs Wj). In the sequel we propose a method to choose the weight Sanctions Ay and show how the asymptotic relative worst case error can be determined in this case. Suppose we can determine coefficients ßj > 0 such that for all z G [0,1] we have m ^ ßiHi{arcsm{z)) ^ 1 , i=i then by Eq. 9 we obtain for all Y G [0,:/r] that f{r)-\ where /(r) = max{l,r/5}. The ffinction is the probability that a is contained in the random 2-slice (Jensen, 1998). Setting Aj(r) = ßif{r), we obtain by Proposition 7, that Hm E[5(Z)] ft^ S{X) . We suggest to chose (/Xi,...,/x^) within the set ^ C [0,oo)'" of all {ßi,...,ßm) such that there exists G [0,1]"= with E^iSy = 1 and ßi = This guarantees that the estimator is asymptotically unbiased for sets X (Z B{0,s) by Corollary 8. In the remainder of this section we show how to determine the asymptotic relative worst case error for given coefficients {ßi,...,ßni) G It is immediate from Eq. 8, that if 1 - vf < XMi^(arcsm(z)) < 1 + Vj^ i=i M (10) for some vf, vf^ > 0 and for all z G [0,1], then 1 - vf < „1, This error bound is independent of the size and shape of X If we know that X C B{0, R) for some R>0, then the worst case error is typically smaller and one can determine a bound using the following proposition. Note that the ffinction / occurring here is the inverse probability that a point of distance r from the origin is contained in Tg. Proposition 10. Suppose that X C B{0,R). Let {ßi,...,ßnj) (^y andset Ur)=ßjf{r), where f{r) = max{l,r/5}. Let e > 0, and let L and M{r) be given as in Lemma 11 below. Deßne r^t = (1 + e/{2L)ysforkeN. Letn be minimal such thatr„ > R. For each k=0,.. .,nletO = YkO < ¥ki < • • • < Vkn,, = n be a partition of [0,:7r], such that \xi/k,i+i — Yk,i\ s is Lipschitz continuous with respect to Xj/G [0,n] with Lipschitz constant 1 r2 ^(//(arcsm(z))) where H = Y.ißiHi and ||-||i denotes the L^-norm on [0,1], It is also Lipschitz continuous with respect tor^ [ro,°o), ro > 5, uniformly in y/, with Lipschitz constant (l/ro)/,, where L = ||//(arcsin(- H + 2 — (//(arcsin(z))) dz The proof of Lemma 11 can be found in the appendix. Proof of Proposition 10. Let Y ^ [0, n] ■ Then with 1G {0,... ,/3jt - 1} such that y ^ iVkh Vk,i+i] we obtain m ^Xiirk)giirk, y) i=i = £ Hrk)gj{rk, 1/A) - £ ?ij{rk)gj{rk, xi/m) i=i i=i + ^?ii{rk)gi{rk,\l/M) i=i 0. Let Mg := [{X® eP) 0 e5]\[(^e eQ)®eW\. Then we obtain with h{-) = {h{P® Or)-h{B(BW,-))+that lN(X) fkia)hin)Q_,iX;dia,n)) < lim sup ^ [ f{^dx{x))dx 6^0+ 0. For general f = f^ - r \NQ can treat f^ and f^ separately to obtain Eq. 2. If /is in addition continuous in all points of dX, we obtain uniform continuity in the following sense. For each t] > 0 there exists a 5 > 0 such that for all x G 5A'nsupp(/) and 7G R"' with Hx-yH < 5 it follows that II f{x) - f{y)\\ < 7], Furthermore x G Mg imphes 0 such that y < C. If X G K^n Ts) e Ts) © trW\, then dist(A',5(A'n Ts)) < tC, where C > 0 is a constant depending only on {B,W). Hence we obtain = t -1 dx,L)dx\ l[{XnTs)etRB\\[{XnTs)®tRW\ < (diXn Ts)®BiO,tC)) + Ct-\^^{d{Tsr\B{Q,&imiX))®B{0,tC)). Applying (Kiderlen and Rataj, 2006, Proposition 4), which is derived from a far-reaching generalization of Steiner's formula (Hug et af, 2004), we obtain t-\^\dX®B{Q,tC)) = f JN(dX) Jo 5(dX,a,n) ^ ^^dX®B(o,tO){a + sn) :dsiXd-i{dX-,d{a,n)) 3 . rtC d{a,n)) JN{dX) Jo i=i where ßi{dX; ■) are the support measures of dX, and 0{dX; a, n) = inf{? > 0 | a+ to G exo(5Z)} is the reach function of dX at {a,n). The support measures have locally finite total variation as Xis gentle and hence the compactness of X yields boundedness of the last term in the above inequality for ? < 1. The same argument can also be applied to t'^ {TsnB{0,diamX)) © B{0, tC)) as 7;n 5(0, tC) is compact gentle. □ Proof of Theorem 3. Let 0 < e < s/2. For 1 G define the continuous fianction : R^ ^ [0,1] as a smoothed version of It^ such that Xsei^) = 0, if xG {Ts)'^ and xiei^) = if ^ G T^-g. Substituting g{x, 1) = fix)xl-e e(^) ill Lemma 5 we obtain JdXnz The right hand side converges pointwise in 1 to f{a)H{ 0 such that 5 U C 5(0, q). Then for t < e/q we have that x G implies xl-e,e{^) = 1 and A-© t{Byj W) C Ts. Therefore xetR{^+h) xeTs+e\Ts^2e where It{x) = 1{x+tRB'zxr\tR{^+'V),x+tRW'ztR{^+-h)\X} and we usedthat;^^_g g(x) = 0 forxG {Tss)'^. Using Eq. 3 this yields < t -1 l[Xetre\\[X®trW\ R=r Choose a constant C such that | /| < C on a compact set containing [(A'n Ts) © trB\ for all r G SO^ and all t<\. Then the last expression in the above inequality is bounded by Cr^ I . . ^dTseB(0,3e){^dx{x))dx. J[XetrB]\[X®trW] ^ ' ^ This integral converges for all r G SCh, by Theorem 1 ast^O+ to IdX The function © rlV, •))+ is bounded by a constant C, independent of r. Therefore, using the dominated convergence theorem, the limit of the above integral as e ^ 0 is bounded by Cjf^{dXn dZ) = 0 for almost all r, which yields the claim. □ Proof of Lemma 11. In order to find a Lipschitz constant with respect to y for m Am X V^) = / X Mi^-(arcsin(z)) dz, 1=1 "'O 1=1 we use partial integration to rewrite the fianction for r > 5as H /n\ s — (//(arcsin(z))) / G^ i^{x)dxdz, 2/ r Jo dz Jo ' where //= ßM Then ßigiir, y) is given by Id r d — (//(arcsm(z))) —G^,/,(x)dxdz. Let q = 5/r and y/j such that \aY,q{z)\ < 1. Then we obtain dxj/ arccos(a,,,,g(z)) gcos(i//') — z sin(i/A) a/sm{\j/)^ - z^ - q^ + 2gzcos(i/A) and hence — arccosia^ ^q{z)) dz dxi/ 1 sin(i/A)2 _ _ g2 2gzcos(i/A) . sin(i/A) The above expression is non-negative and bounded by A/l-g2. Therefore {d/dxi/)Yf=ißigi{r,xi/) is bounded by -(//(arcsm(z))) dz < -(//(arcsm(z))) hence M(r) is a Lipschitz constant for (r/s) E™ 1 ßigi{r, y) with respect to y/To find a Lipschitz constant for X™ i h{i)gi{r, Xf/) with respect to r we first differentiate with respect to r and obtain s Jo //(arcsin(z)) G^^^j^z) dz rd ^-sTrlo The first term of the above expression is bounded by l/r||//(arcsin(-))||cx,. In order to find a bound for the second term we use partial integration to rewrite 'ZiLi ßigiii', V) for r > 5 as in Eq. 19. Then {d/dr) E™ 1 ßigi{r, i/a) is given by -"'I 5 r^ /•Id d Let i/A, z such that | a^ ^/^(z) | < 1, then we have (19) ^arccos(a^,,/^(z)) and r2-^sin(i/A)2 - z2 - ^ +2fzcos{\j/) ^arccos(a^^s/^(z)) dz / \ z-fcos(i/A) = ^ arctan sin(i/a)2-z^ - ^+2^zcos(i/a) The term on the right-hand side of the above equation is bounded in absolute value by 5 n fir Therefore, for r > ro, the Sanction (r/ 5) X^i ßigi{i~, V) has Lipschitz constant {l/ro)L in r uniformly in XI/. □ REFERENCES FremlinDH (1997). Skeletons and central sets. Proc London Math Sog 74:701-20. Gutkowski P, Jensen EBV, Kiderlen M (2004). Directional analysis of digitized three-dimensional images by configuration counts. J Microsc 216:175-85. Hug D, Last G, Weil W (2004). A local Steiner-type formula for general closed sets and applications. Math Z 246:237-72. Jensen EBV (1998). Local stereology. London: World Scientific. Kiderlen M, Jensen EBV (2003). Estimation of the directional measure of planar random sets by digitization. Adv Appl Prob SGSA 35:583-602. Kiderlen M, Rataj J (2006). On infinitesimal increase of volumes of morphological transforms. Mathematika 53:103-27. Klenke A (2006). Wahrscheinlichkeitstheorie. Berlin: Springer. Klette R, Rosenfeld A (2004). Digital geometry: Geometric methods for digital picture analysis. San Francisco: Morgan Kaufmann. Lindblad J (2005). Surface area estimation of digitized 3D objects using weighted local configurations. Image Vision Comput 23:111-22. Lindblad J, Nystroom I (2002). Surface area estimation of digitized 3D objects using local computations. In: Discrete Geometry for Computer Imagery: 10th Int Conf, DGCI2002. Berlin: Springer. Schladitz K, Ohser J, Nagel W (2006). Measuring intrinsic volumes in digital 3D images. In: Kuba A, Nyiil LG, Palagyi K, eds., 13 th Int Conf Discrete Geom Computer Imagery. Heidelberg: Springer. Yap CK (2000). Fundamental problems of algorithmic algebra. New York: Oxford University Press. Ziegel J, Kiderlen M (2010). Estimation of surface area and surface area measure of three-dimensional sets from digitizations. Image Vision Comput 28:64-77.