# Stereological estimation of surface area from digital images.

INTRODUCTIONAlgorithms 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 efficient. Lindblad (2005) suggested a local estimator based on voxel blocks of size 2 x 2 x 2 inspired by the common marching cubes algorithm, but without an explicit reconstruction of the surface. There are [2.sup.8] = 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

[??](X) = [254.summation over (i = 1)] [[lambda].sub.i] [N.sub.i], (1)

where [[lambda].sub.i] is the contribution and [N.sub.i] is the total number of occurrences of configuration i. Several natural choices of the weights [[lambda].sub.1], ..., [[lambda].sub.254] can be found in the literature (Lindblad and Nystrom, 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 informative 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 s > 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 function 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 functions 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 s; 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 2 x 2 x 2 configurations. In the final section we specialize these results to the scaled standard lattice L = [tZ.sup.3], 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 [S.sup.d-1] we denote the unit sphere in d-dimensional real coordinate space [R.sup.d]. The standard scalar product on [R.sup.d] is (*,*) with associated norm [parallel]*[parallel]. By a k-subspace we mean a k-dimensional linear subspace of [R.sup.d]. Let A,B [subset] [R.sup.d] The reflection of A at the origin is denoted by [??] = {-x|x [member of] A}, its complement by [A.sup.c] = [R.sup.d]\A, and its topological boundary by [partial derivative]A. We write A [direct sum] B = {a + b|a [member of] A,b [member of] B} for the Minkowski sum of A and B, and A [??] B = {x [member of] [R.sup.d]|x + [??] [subset] A} for the dilatate of A by [??]. The positive part of a real valued function f is denoted by [f.sup.+] = max(f, 0). The support function of a convex body K in R is denoted by h(K, *). We use this notion also for compact sets A, A [not equal to] [empty set], defining h(A, *) = h(conv(A), *), where conv(A) is the convex hull of A. The exoskeleton exo(A) of a closed set A is the set of all z [member of] [A.sup.c], which do not have a unique nearest point in A. The set exo(A) is measurable and has Lebesgue measure zero (Fremlin, 1997).

A closed set X [subset] [R.sup.d] is gentle if for [H.sup.d-1]-almost all x [member of] [partial derivative] X there are two non-degenerate open balls touching in x such that one of them is contained in X and the other in Xc, and if also [H.sup.d-1] (N([partial derivative]X [intersection] B x [S.sup.d-1])) < [infinity] for all bounded Borel sets B [subset] [R.sup.d].

Here [H.sup.k] is the k-dimensional Hausdorff measure in [R.sup.d] 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.sup.d]) 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 [C.sub.d-1] (X, *) be the image measure of [H.sup.d-1] on [partial derivative]X under the map a [right arrow] (a,n(a)). The measure [C.sub.d-1](X, *) vanishes outside N{X). Let [[xi].sub.[partial derivative]X] : [R.sup.d]\exo([partial derivative]X) right arrow] [partial derivative]X denote the metric projection. The following theorem is a generalization of (Kiderlen and Rataj, 2006, Theorem 1).

Theorem 1. Let X [??] [R.sup.d] be a closed gentle set, f: [R.sup.d] [right arrow] R a compactly supported bounded measurable function and B, W and P, Q four non-empty compact subsets of[R.sup.d]. Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. (2)

If f is in addition continuous at all points in [partial derivative]X, then f([[xi].sub.[partial derivative]X](x)) can 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 f is an indicator function. Their result can be generalized to Theorem 1 essentially by applying the monotone convergence theorem.

Let [x.sub.1], ..., [x.sub.d] be a basis of [R.sup.d] and let

L = {[n.sub.1][x.sub.1] + *** + [n.sub.d] [x.sub.d]|[n.sub.1], ..., [n.sub.d] [member of] Z}

be the lattice generated by this basis. A given lattice L is generated by infinitely many different bases, but the volume of the fundamental cell [C.sub.0] = [0, [x.sub.1]] [direct sum] *** [direct sum] [0, [x.sub.d] depends only on L and not on the basis chosen; see for example (Yap, 2000). This number is denoted by det(L). If [xi], is a uniform random variable in [C.sub.0], then the random lattice [xi], + L is a stationary random lattice. The distribution of [xi], + L does not depend on the choice of Co. We now define the digitization of a set X [subset] [R.sup.d]. The points of [xi] + L are interpreted as voxel midpoints of a digital image, where each voxel is a translate of [C.sub.0]. Often, L is the orthogonal standard lattice [Z.sup.d] 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 [xi] + L on the other hand, the Gauss digitization is determined by the set X [intersection] ([xi] + 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 > 0 and X [intersection] t ([xi] + L) is called the digitization of X (with resolution 1/t).

In the following we only consider compact gentle sets X. Let the function f be measurable non-negative or integrable. Let [xi] + L be a stationary random lattice, and let B, W [??] 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 > 0 define

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

For f [??] 1 the random variable [N.sub.i] counts the number of all translations of the pattern (tB,tW) in the digitization X [intersection] t([xi] + L). Calculation shows that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (3)

Corollary 2. Let X [??] [R.sup.d] be a compact gentle set. Let f be a locally bounded measurable function, which is continuous on [partial derivative]X and let B and W be two non-empty finite subsets of a lattice L. Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Proof. Let C be a compact set such that [X[??] t[??]]\[X [direct sum] t[??]] [subset] C for all t > 0 smaller than some fixed [t.sub.0] > 0 and X [subset] C. Replacing f, P, B, Q and W in Theorem 1 by [1.sub.c] f, {0}, [??], {0} and [??], respectively, yields the claim.

COMBINING LOCAL AND DIGITAL STEREOLOGY

In the following we restrict ourselves to [R.sup.3]. The results can be generalized to [R.sup.d], d [greater than or equal to] 4, in a straightforward manner. We prefer to present them only in [R.sup.3] in order to keep the notation concise.

Denote the standard basis vectors in [R.sup.3] by [e.sub.1], [e.sub.2], [e.sub.3]. Let R be a random proper rotation with distribution given by the normalized Haar measure on the rotation group S[O.sub.3]. Fix the 2-subspace [l.sub.0] = span([e.sub.1], [e.sub.2]). We define the random 2-subspace L = [Rl.sub.0]. It is uniformly distributed in the set L of all 2-subspaces of [R.sup.3]. Let [mu] be the distribution of L. For l [member of] L and s > 0 define [T.sub.s] = [T.sub.s](l) = l [direct sum] 5(0, s). The set [T.sub.s] = [T.sub.s](L) = R([l.sub.0] [direct sum] B(0,s)) = L [direct sum] B(0,s) is called a random 2-slice with thickness 2s. It will be clear from the context whether [T.sub.s] refers to the deterministic 2-slice [T.sub.s](l) or the random 2-slice [T.sub.s](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 X [??] [R.sup.3] be a compact gentle set. Let R be a random proper rotation and let [xi] + L be a stationary random lattice, which is independent of R. Let B, W [subset] L be two non-empty finite subsets of the lattice L and f a continuous non-negative function on [R.sup.3]. The weighted sum [N.sup.f.sub.t] of occurrences of (B, W) in the digitization of X, which lie entirely in [T.sub.s], is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

for t > 0. It satisfies

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

where H([PHI]) is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (4)

and [[PHI].sup.X.sub.a,l] is the angle between I and the outer normal n(a) of X at a [member of] [partial derivative]X.

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 [N.sup.f.sub.t] in Theorem 3 and [[??].sup.g.sub.t] in the following Proposition 4 is that the latter also counts configurations which do not lie entirely in the slice [T.sub.2]. These are of course not observable in practice.

Proposition 4. Let X [??] [R.sup.3] be a compact gentle set. Let R be a random proper rotation and let [xi] + L be a stationary random lattice, which is independent of R. Let B, W [subset] L be two non-empty finite subsets of L. Let g: [R.sup.3] x L right arrow] R be a non-negative bounded measurable function such that g(*, l) is continuous for all l [member of] L. Then the number

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

for t > 0 satisfies

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (5)

for r [member of] S[O.sub.3] with h = [(-h(B [direct sum] [??], *)).sup.+].

The proof of Proposition 4 can be found in the appendix. It uses that given the rotation R, Corollary 2 can be applied to [[??].sup.g.sub.t]. In the following two Lemmas we derive a formula for the expected value of [F.sub.g](R).

Lemma 5. Let g : [R.sup.3] x L right arrow] R be a non-negative measurable function such that g(*, l) is continuous for all l [member of] L. Let [F.sub.g](r) be given by Eq. 5. Then the conditional expectation E[[F.sub.g](R)|L = l] can be expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

for [mu]-almost all l [member of] L, where H([PHI]) is given by Eq. 4.

Proof. For [mu]-almost all l [member of] L we have that [H.sub.2]([partial derivative]X[intersection] [partial derivatives][T.sub.s]) = 0. This implies that there is a unique outer unit normal n(a) for [H.sup.2]-almost all a [member of] [partial derivatives](X [intersection] [T.sub.s]). Hence we obtain for r [member of] S[O.sub.3]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

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

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Recall that [mu]-almost surely [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. The claim now follows from Lemma 6.

Lemma 6. For n [member of] [S.sup.2] and l [member of] L we have

E[h([R.sup.-1]n)|L = l] = H([PHI]),

where [PHI] is the angle between n and l, and H([PHI]) is given by Eq. 4.

Proof. Fix [rho] [member of] S[O.sub.3] such that [rho]l = [l.sub.0]. Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

We have <n,l> = <[rho]n,[l.sub.0]>, hence the last conditional expectation in the above equation can be written as the normalized integral over the two small circles [subset] [S.sup.2] parallel to [l.sub.0] at height [+ or -] sin [PHI] with radius cos [PHI], where [PHI] is the angle between n and l.

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 X [member of] [R.sup.3] be a compact gentle set. Suppose we observe X [intersection] R[([l.sub.0]+B(0,s)) [intersection] t([xi] + L)] for some random proper rotation R, [xi] + L a stationary random lattice, which is independent of R, and s,t > 0. Let ([B.sub.i], [W.sub.i]),i = 1, ..., m be boundary configurations of L, i.e. [B.sub.i], [W.sub.i] are non-empty, disjoint finite subsets of L with [B.sub.i] [union] [W.sub.i] = [C.sub.0] [intersection] L, where [C.sub.0] is a fixed fundamental cell of L. We define the following estimator for surface area

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (6)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] as defined in Theorem 3 with B = [B.sub.i],W = [W.sub.i]. The continuous functions [[lambda].sub.i]: [0, [infinity]) [right arrow] [0,[infinity]) have to be suitably chosen according to the choice of ([B.sub.i], [W.sub.i]). We give an example for 2 x 2 x 2-configurations in the concluding section.

Proposition 7. We have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (7)

where [[PSI].sup.X.sub.a] [member of] [0, [pi]] is the angle between a and the outer normal n(a) of X at a [member of] [partial derivative]X. The function [g.sub.i](r, [psi]) for [psi] [member of] [0, [pi]] and r [member of] [0,[infinity]) is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

With the two sided cut-off function x [right arrow] [x.sup.*] = min{1, max(-1, x}}, the function [G.sub.[psi], q](z) for [psi] [member of] (0, [pi]) and q [member of] (0,1] can be written as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. For [psi] [member of] {0, [pi]} we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Proof. By Theorem 3 we obtain that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (8)

where [H.sub.i] is given by Eq. 4 for (B,W) = ([B.sub.i], [W.sub.i]). Using Fubini's theorem we can rewrite the right-hand side of Eq. 8 as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

We have [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], which is a tedious but not difficult calculation.

Note that for all [psi] [member of] [0, [pi]] and all r [member of] (s, [infinity])

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (9)

(Jensen, 1998). We assume from now on that none of the functions [H.sub.i],i = 1,..., m, is identically zero. This is fulfilled, if and only if [B.sub.i] and [W.sub.i] can be strictly separated by a hyperplane for all i = 1, ..., m, or, in other words, ([B.sub.i],[W.sub.i]) is an informative configuration.

The following two corollaries to Proposition 7 show that the weight functions [[lambda].sub.i] can be chosen to yield an asymptotically unbiased estimator [??](X) in the case where X is sufficiently small compared to s or when X is a ball centred at the origin with unknown radius.

Corollary 8. Let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] with coefficients [a.sub.1], ..., [a.sub.m] [member of] R 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 [H.sub.i] [greater than or equal to] 0. The assumption that [H.sub.i] [??] 0 yields [g.sub.i](0,0) > 0. As [g.sub.i](r, [psi]) = [g.sub.i](0,0) for all (r, [psi]) [member of] [0, s] x [0, n], Eq. 7 implies the claim.

Corollary 9. LetX be a ball centred at the origin with unknown radius. Suppose that the sets Bh Wt are such that [g.sub.i](r,0) > 0 for all r [member of] [0, [infinity]) and i = 1, ..., m. Choosing [[lambda].sub.i](r) = [a.sub.i][g.sub.i][(r, 0).sup.-1], where [[SIGAMA].sup.m.sub.i = 1] [a.sub.i] = 1 yields an asymptotically unbiased estimator [??](X) of S(X).

Proof. If X is a ball centred at the origin we have [[psi].sup.X.sub.a] = 0 for all a [member of] [partial derivative]X. The claim follows from Eq. 7.

The condition [g.sub.i](r,0) > 0 in Corollary 9 is equivalent to requiring that the support of [H.sub.i] contains an interval [0, [epsilon]) for some [epsilon] > 0.

For general shapes we cannot expect to obtain an unbiased estimator of the form Eq. 7. A suitable choice of [[lambda].sub.i] for r > s will strongly depend on the choice of the pairs ([B.sub.i],[W.sub.i]). In the sequel we propose a method to choose the weight functions [[lambda].sub.i] and show how the asymptotic relative worst case error can be determined in this case. Suppose we can determine coefficients [[mu].sub.i] [greater than or equal to] 0 such that for all z [member of] [0,1] we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

then by Eq. 9 we obtain for all [psi] [member of] [0, [pi]] that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], where f(r) = max(1,r/s}. The function f[([parallel] a [parallel]).sup.-1] is the probability that a is contained in the random 2-slice [T.sub.s] (Jensen, 1998). Setting [[lambda].sub.i](r) = [[mu].sub.i] f (r), we obtain by Proposition 7, that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

We suggest to chose ([[mu].sub.1], ..., [[mu].sub.m] within the set L [??] [[0, [infinity]).sup.m] of all ([[mu].sub.1], ..., [[mu].sub.m]) such that there exists [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] and [[mu].sub.i] = [a.sub.i] [g.sub.i] [(0, 0).sup.-1]. This guarantees that the estimator is asymptotically unbiased for sets X [??] 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 ([[mu].sub.i], ..., [[mu].sub.m]) [member of] L. It is immediate from Eq. 8, that if

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (10)

for some [v.sup.m.sub.1], [v.sup.M.sub.1] > 0 and for all z [member of] [0,1], then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (11)

This error bound is independent of the size and shape of X. If we know that X [??] 5(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 function f occurring here is the inverse probability that a point of distance r from the origin is contained in [T.sub.s].

Proposition 10. Suppose that X [??] B(0, R). Let ([[mu].sub.i], ..., [[mu].sub.m]) [member of] L and set

[[lambda].sub.i] (r) = [[mu].sub.i] f(r),

where f(r) = max{1, r/s}. Let [epsilon] > 0, and let L and M(r) be given as in Lemma 11 below. Define [r.sub.k] = [(1 + [epsilon]/ (2L)).sup.k] s for k [member of] N. Letn be minimal such that [r.sub.n] [greater than or equal to] R. For each k = 0, ..., n let 0 = [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] be a partition of [0, [pi]], such that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] for all l = 0, ..., [n.sub.k]. Set

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. (12)

For a proof of Proposition 10 we make use of a Lipschitz result for the integrand in Eq. 7.

Lemma 11. The function (r/s)[[SIGMA].sup.m.sub.i = 1] [[mu].sub.i] [g.sub.i] (r, [psi]) for r > s is Lipschitz continuous with respect to [psi] [member of] [0, [pi]] with Lipschitz constant

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

where H = [[SIGMA].sub.i] [[mu].sub.i][H.sub.i] and [parallel] * [[parallel].sub.1] denotes the [L.sup.1]-norm on [0,1]. It is also Lipschitz continuous with respect to r [member of] [[r.sub.0], [infinity]), [r.sub.0] [greater than or equal to] s, uniformly in [psi], with Lipschitz constant (1/[r.sub.0])L, where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

The proof of Lemma 11 can be found in the appendix.

Proof of Proposition 10. Let [psi] [member of] [0, [pi]]. Then with l [member of] (0, ..., [n.sub.k] - 1} such that [psi] [member of] [[[psi].sub.kl], [[psi].sub.k, l + 1] we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Hence for r [member of] [s, R] with k such that r [member of] [[r.sub.k], [r.sub.k+1]) we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

By the same arguments one obtains the analogous lower bound for [[SIGMA].sup.m.sub.i = 1] [[lambda].sub.i] (r) [g.sub.i] (r, [psi]). Using Eq. 7 this yields Ineq. 12.

COEFFICIENTS FOR 2 x 2 x 2 CONFIGURATIONS

Recall that a configuration of size 2 x 2 x 2 is a pair (B,W) of non-empty disjoint subsets of L = [Z.sup.3], such that B [union] W = [Z.sup.3] [intersection] [[0, 1].sup.3]. It is informative, if there is a hyperplane, which strictly separates B and W. In this section we want to investigate the surface area estimator Eq. 6 in the case, where ([B.sub.i], [W.sub.i]) runs through the family of all 102 informative configurations. These configurations are thoroughly investigated in Ziegel and Kiderlen, (2010). In particularthe functions [(-h([B.sub.i] + [[??].sub.i])).sup.+] are explicitly given, hence we can numerically determine the functions [H.sub.i]. As in Ziegel and Kiderlen (2010) we classify the informative configurations into five types, depending on the number and position of black points B or white points W. A configuration of type one has exactly one black point or exactly one white point, a configuration of type two has exactly two black points or exactly two white points, and a configuration of type three has exactly three black points or exactly three white points. Configurations of type four and five have exactly four white and four black points, which are afnnely dependent in the case of type four, and afhnely independent in the case of type five.

For configurations of type one, all functions [H.sub.i] are identical and we denote them by [H.sup.1]. The function [H.sup.1](arcsin(z)) for z [member of] [0,1] is shown in Fig. 1. All functions [H.sub.i](arcsin(z)), z [member of] [-1,1], are symmetric with respect to the origin, which is why we only display them for values z [member of] [0,1]. For configurations of type two, three and four there are two different functions [H.sub.i] occurring per type. We denote them by [H.sup.2,1], [H.sup.2,2], [H.sup.3,1], [H.sup.3,2] and [H.sup.4,1], [H.sup.4,2] respectively; see Figs. 2, 3 and 4. For configurations of type five all functions Hi coincide and we denote them by [H.sup.5], which is displayed in Fig. 5. Fig. 6 shows all functions [H.sub.i] scaled by the number of their occurrence amongst all functions [H.sub.i] induced by informative configurations. The numbers of occurrence of the functions [H.sub.1] [H.sup.2,1] [H.sup.2,2], [H.sup.3,1] [H.sup.3,2], [H.sup.4,1], [H.sup.4,2] and [H.sup.5] are given in the last column of Table 2.

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

We define [[eta].sub.i] = [g.sub.i][(0,0).sup.-1]. The values [[eta].sub.i], for informative configurations are given in Table 1, where [[eta].sub.i] = [[eta].sup.1] whenever ([B.sub.i], [W.sub.i]) is a configuration of type one, and analogously in the other cases.

We have seen that the estimator in Eq. 6 with

[[lambda].sub.i](r) = [[mu].sub.i] f (r) (13)

is unbiased for S(X) if X is a subset of B(0,s), whenever ([[mu].sub.1], ..., [[mu].sub.102]) [member of] L. In Ziegel and Kiderlen, 2010) a one-parameter family of coefficients ([[mu]'.sub.1](u)), ..., [[mu]'.sub.102](u)), u [member] [0,1] was derived that minimizes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

for each u [member of] [0,1]. Adapting the coefficients ([[mu]'.sub.1](u), ..., [[mu]'.sub.102] (u)) to yield an asymptotically unbiased estimator for spherical shapes, we obtain the family of coefficients ([[mu].sub.1] (u), ..., [[mu].sub.102](u)), u [member of] [0, 1] given in Table 2; for details see Ziegel and Kiderlen (2010). It turns out that ([[mu].sub.1](0), ..., [[mu].sub.102](0)) [member of] L or in other words

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

We therefore suggest to set [[mu].sub.i] = [[mu].sub.i](0) and [[lambda].sub.i] as in Eq. 13. Hence we obtain the estimator

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], (15)

which has the same formal structure as the estimator in equation Eq. 1. However, the ordinary configuration counts [N.sub.i] in the classical estimator are replaced by the weighted configuration counts [W.sup.i.sub.t] given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

where f(r) = max(l,r/s} is the inverse sampling probability. The fusion of digital and local stereology becomes apparent in the structure of this estimator: While the weights [[mu].sub.i](0) stem from the ordinary digital surface area estimation, the weighting of the individual configurations in [W.sup.i.sub.t] with f reflects the local design. With this choice we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], (16)

using the definition Eq. 4 of the [H.sub.i], hence [v.sup.m.sub.1] = 0.0461, [v.sup.M.sub.1] = 0.0225 fulfill Eq. 10. Therefore, by Ineq. 11, the upper and lower bound in Ineq. 16 are bounds for the asymptotic relative worst case error of [??](X), as given in Eq. 15, independent of the shape and size of X. In the next paragraph we illustrate the application of Proposition 10 with a simulation example. This yields better error bounds but requires information about the size and preferably also the shape of X.

We consider a cylinder with radius 1 and height 2 centred at 0 which is contained in B(0, [square root of 2]). We observe an isotropic slice of thickness 5 = 1. For a lattice distance of t = 0.055 we obtain a mean estimated surface area of 18.115 with varianc[e.sub.1].146 in 1000 Monte Carlo simulations. This corresponds to a mean relative error of 3.8%. For t = 0.020 the mean estimated surface area in 1000 simulations is 18.553 with varianc[e.sub.1].205 and mean relative error 1.6%. The asymptotic relative mean error for sets X with X [??] B(0, [square root of (2)]) and [[psi].sup.X.sub.a] [member of] [0, [pi]/4] for all a [member of] X is less than 1.2%. We determined this value numerically using Proposition 10. We obtain [parallel](d/dz))H(arcsin(z)))[parallel].sub.1] = 5.7002, where H = [[SIGMA].sub.i][[mu].sub.i] [H.sub.i] and L = 13.3768 as defined in Lemma 11. Note that the information [[psi].sup.X.sub.a] [member of] [0, [pi]/4] reduces the asymptotic relative mean error as it is then sufficient to work with partitions [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] in Proposition 10. With [epsilon] = = 0.0030 we obtain [v.sup.M.sub.2] = 0.0042 and [v.sup.m.sub.2] = 0.0082.

In summary this simulation example indicates that the bias of the simple linear surface area estimator Eq. 6 with weight functions Eq. 13 is of reasonable magnitude. It appears that the theoretical asymptotic error bounds are often too optimistic, and should only be used in the case where good to very good resolution images are available.

ACKNOWLEDGEMENTS

The first author was generously supported by RiskLab Switzerland. The second author was supported by the Danish Council of Strategic Research.

APPENDIX

In this appendix we present the proofs of Theorems 1 and 3, Proposition 4 and Lemma 11.

Proof of Theorem 1. Let C [subset][R.sup.d] be a bounded Borel set. Then Eq. 2 holds for f = [1.sub.C] by Kiderlen and Rataj (2006, Theorem 1). It is immediate that Eq. 2 also holds for compactly supported measurable step functions. For a non-negative compactly supported bounded measurable function f, let [([f.sub.k]).sub.k[member of]N], [([g.sub.k]).sub.k[member of]N] be sequences of step functions such that [f.sub.k] [up arrow] f and [g.sub.k] [down arrow] f and [f.sub.k] [greater than or equal to] 0. Let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. Then we obtain with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (17)

Using the monotone convergence theorem we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Note that for applying the monotone convergence theorem to the right-hand side of Ineq. 17 we can assume that [U.sub.k [member of] N] supp ([g.sub.k]) is compact and that the sequence [([g.sub.k]).sub.k [member of] N] is uniformly bounded. The same argument holds if we take lim [inf.sub.[epsilon][right arrow] 0+] in Ineq. 17 and hence the claim is shown for f [greater than or equal to] 0. For general f = [f.sup.+] - [f.sup.-] we can treat [f.sup.+] and [f.sup.-] separately to obtain Eq. 2. If f is in addition continuous in all points of [partial derivative]X, we obtain uniform continuity in the following sense.

For each [eta] > 0 there exists a [delta] > 0 such that for all x [member of] [partial derivative]X [intersection] supp(f) and y [member of] [R.sup.d] with [parallel]x-j[parallel] [less than or equal to] [delta] it follows that [parallel]f(x) - f(y)[parallel] < [eta] Furthermore x [member of] [M.sub.[epsilon]] implies [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. This implies the second claim.

Proof of Proposition 4. For [mu]-almost all l [member of] L we have that [H.sup.2] ([partial derivative]X [intersection] [partial derivative][T.sub.s]) = 0, which can be seen using Fubini's theorem. This implies that the set X [intersection] [T.sub.s] is compact gentle for [mu]-almost all l [member of] L. Applying Corollary 2 we obtain that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] for t [right arrow] 0+ pointwise for almost all r [member of] S[O.sub.3] as [(-h(rB [direct sum] r[??}, n)).sup.+] = h([r.sup.-1] n). We claim that the conditional expectation [t.sup.2]E[[??].sup.g.sub.t]|R] is uniformly bounded for t [less than or equal to] 1, hence Lebesgue's dominated convergence theorem yields the assertion. In fact Eq. 3 implies

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

By assumption there is a constant C > 0 such that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], then dist(x, [delta](X [intersection] [T.sub.s])) [less than or equal to] tC', where C' > 0 is a constant depending only on (B,W). Hence we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Applying (Kiderlen and Rataj, 2006, Proposition 4), which is derived from a far-reaching generalization of Steiner's formula (Hug et al., 2004), we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

where [[mu].sub.i]([partial derivative]X;*) are the support measures of [partial derivative]X, and [delta]([partial derivative]X; a, n) = inf(t [greater than or equal to] 0 | a + tn [member of] exo([partial derivative] X)} is the reach function of [patial derivative] X at (a, n). The support measures have locally finite total variation as X is gentle and hence the compactness of X yields boundedness of the last term in the above inequality for t [less than or equal to] 1. The same argument can also be applied to [t.sup.-1] [H.sup.3] ([partial derivative]([T.sub.s] [intersection]B(0,diamX)) [direct sum] B(0, tC')) as [T.sub.s] [intersection] B(0, tC') is compact gentle.

Proof of Theorem 3. Let 0 < [epsilon] < s/2. For f [member of] L, define the continuous function [x.sup.l.sub.s,[epsilon]]: [R.sup.3] [right arrow] [0,1] as a smoothed version of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] such that [X.sup.l.sub.s,[epsilon]] (x) = 0> if x [member of] [([T.sub.s]).sup.c] and [X.sup.l.sub.s, [epsilon]] (x) = 1, if x [member of] [T.sub.s-[epsilon]]. Substituting g(x,l) = f(x) [x.sup.l.sub.s-[epsilon], [epsilon]] (x) in Lemma 5 we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

The right hand side converges pointwise in l to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

as [epsilon] [right arrow] 0. It is also bounded independently of [epsilon] and l, hence by the dominated convergence theorem [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] converges to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. (18)

Proposition 4 yields that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

is also given by Eq. 18. It remains to show that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Choose q > 0 such that B [union] W [??] B(0, q). Then for t < [epsilon]/q we have that x [member of] [T.sub.s-2[epsilon]] implies [x.sup.L.sub.s-[epsilon], [epsilon]](x) = 1 and x + t(B [union] W) [??] [T.sub.s]. Therefore

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] and we used that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. Using Eq. 3 this yields

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Choose a constant C such that [absolute value of f] [less than or equal to] C on a compact set containing [(X [intersection] [T.sub.s]) [??] tr[??]] for all r [member of] S[O.sub.3] and all t [less than or equal to] 1. Then the last expression in the above inequality is bounded by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

This integral converges for all r [member of] S[O.sub.3] by Theorem 1 as t [right arrow] 0+ to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

The function [(-h(rB [direct sum] r[??],*)).sup.+] is bounded by a constant C', independent of r. Therefore, using the dominated convergence theorem, the limit of the above integral as [epsilon] [right arrow] 0 is bounded by C' [H.sup.2]([partial derivative]X [intersection] [partial derivative][T.sub.s]) = 0 for almost all r, which yields the claim.

Proof of Lemma 11. In order to find a Lipschitz constant with respect to yr for

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

we use partial integration to rewrite the function for r > s as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (19)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.].

Let q = s/r and [psi], q, z such that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. Then we obtain

[partial derivative/[partial derivative][psi] arccos ([[alpha].sub.[psi], q](z))

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

and hence

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

The above expression is non-negative and bounded by [square root of (1-[q.sup.2])]. Therefore [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] is bounded by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

hence M(r) is a Lipschitz constant for (r/s)[[SIGMA].sup.m.sub.i = 1] [[mu].sub.i] [g.sub.i](r, [psi]) witn respect to [psi].

To find a Lipschitz constant for [[SIGMA].sup.m.sub.i = 1] [[lambda].sub.i](r)[g.sub.i](r, [psi]) with respect to r we first differentiate with respect to r and obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

The first term of the above expression is bounded by 1/r[parallel]H(arcsin(*))[[parallel].sub.[infinity]. In order to find a bound for the second term we use partial integration to rewrite [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] for r > 5 as in Eq. 19. Then [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Let [psi], r, z such that [absolute value of[[alpha].sub.[psi], s/r](z)] [less than or equal to] 1, then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

The term on the right-hand side of the above equation is bounded in absolute value by

s/[r.sup.2] [pi]/2.

Therefore, for r > [r.sub.0], the function [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] has Lipschitz constant (l/[r.sub.0])L in r uniformly in yr.

(Accepted May 12, 2010)

REFERENCES

FremlinDH (1997). Skeletons and central sets. Proc London MathSoc 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, Nystrom I (2002). Surface area estimation of digitized 3D objects using local computations. In: Discrete Geometry for Computer Imagery: 10th Int Conf, DGCI 2002. Berlin: Springer.

Schladitz K, Ohser J, Nagel W (2006). Measuring intrinsic volumes in digital 3D images. In: Kuba A, Nyul LG, Palagyi K, eds., 13th 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.

JOHANNA ZIEGEL (1) AND MARKUS KIDERLEN (2)

(1) Department of Mathematics and Statistics, The University of Melbourne, Victoria 3010, Australia; (2) Department of Mathematical Sciences, University of Aarhus, Ny Munkegade Build. 1530, DK-8000 Aarhus C, Denmark e-mail: jziegel@unimelb.edu.au, kiderlen@imf.au.dk

Table 1. Values of [[eta].sub.i] for informative configurations. [H.sub.1] [[eta].sub.i] = [([integral].sup.1.sub.0] [H.sub.i](arcsin)(z))dz).sup.-1] [H.sup.1] [[eta].sup.1] = 38.119 [H.sup.2,1] [[eta].sup.2,1] = 48.179 [H.sup.2,2] [[eta].sup.2,2] = 48.179 [H.sup.3,1] [[eta].sup.3,1] = 171.187 [H.sup.3,2] [[eta].sup.3,2] = 171.187 [H.sup.4,1] [[eta].sup.4,1] = 25.856 [H.sup.4,2] [[eta].sup.4,2] = 25.856 [H.sup.5] [[eta].sup.5] = 114.825 Table 2. One-parameter family of coefficients [[mu].sub.i] for u [member of] [0,1]. The number in the last column indicates the number of occurrences of the function [H.sup.1], [H.sup.2,1], [H.sup.2,2], [H.sup.3,1], [H.sup.3,2], [H.sup.4,1], [H.sup.4,2], [H.sup.5] respectively amongst all Junctions Hi for informative 2 x 2 x 2-configurations. u [member of] [H.sup.1] [[mu].sub.i](u), [0, 1] [H.sup.1] [[mu].sup.1] (u) = (1.652/2)u 16 [H.sup.2,1] [[mu].sup.2,1] (u) = 0.675 8 [H.sup.2,2] [[mu].sup.2,2] (u) = 0.675 16 [H.sup.3,1] [[mu].sup.3,1] (u) = 1.168 - (1.652/4)u 32 [H.sup.3,2] [[mu].sup.3,2] (u) = 1.168 - (1.652/2)u 16 [H.sup.4,1] [[mu].sup.4,1] (u) = 0.954 4 [H.sup.4,2] [[mu].sup.4,2] (u) = 0.954 2 [H.sup.5] [[mu].sup.5] (u) = 1.652 (1-u) 8

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Original Research Paper |
---|---|

Author: | Ziegel, Johanna; Kiderlen, Markus |

Publication: | Image Analysis and Stereology |

Article Type: | Report |

Geographic Code: | 4EUDE |

Date: | Jul 1, 2010 |

Words: | 7882 |

Previous Article: | Image analysis of pearlite spheroidization based on the morphological characterization of cementite particles. |

Next Article: | On the specific area of inhomogeneous Boolean models. Existence results and applications. |

Topics: |