Stereological estimation of surface area from digital images

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. The method is based on a result in (Kiderlen M, Rataj J (2006). On infinitesimal increase of volumes of morphological transforms. Mathematica 53: 103–127) which is generalized in this paper. The proposed surface area estimator is asymptotically unbiased in the case of sets contained in the ball centred at the origin with radius s and in the case of balls 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× 2× 2-configurations.


Introduction
Methods of local stereology are in world-wide use in microscopic studies of biological tissue.The sampling designs of local stereology are based on sections through a fixed reference point of the structure under consideration.It is often encountered that the observations within the sections are obtained by methods of digital stereology, such as point counting for area estimation.In this paper we are taking the idea of combining local and digital stereology a step further.We derive an estimator for the surface area of three-dimensional objects based on discrete binary images obtained by a local sampling design under very weak assumptions on the shape.
The estimator is based on a generalization of an asymptotic result in (Kiderlen and Rataj, 2006).The original result was successfully used for boundary length estimation in (Kiderlen and Jensen, 2003;Jensen and Kiderlen, 2003) and for surface area estimation in (Gutkowski et al., 2004;Ziegel and Kiderlen, 2009).Let X ⊆ R 3 be a compact gentle set; see the next section for details.Suppose we observe a digitization of X ∩ T s , where T s is an isotropic thickened 2-dimensional linear subspace of thickness 2s.Then the weighted number N f t of occurrences of a boundary configuration (B, W ) of black B and white W points in the digitized picture of X ∩ T s behaves in mean like ct −2 as t → 0, where t > 0 is the lattice distance.The estimator we present relies on the fact that the normalization constant c can be given explicitly as the double integral over the set L of all 2-dimensional linear subspaces of R 3 and over the observed part of the boundary of X.The value f (a) determines the weight given to an observed configuration located at a and H is a positive function on [0, π], which depends only on (B, W ). The function φ X a,l gives the angle between l and the outer normal of ∂X at a.The constant det L depends on the lattice L chosen for the digitization.
For surface area estimation we choose a set of boundary configurations (B i , W i ), i = 1, . . ., m and weight functions λ i on the positive halfline to estimate S(X) by where N λ i t is the λ i -weighted number of occurrences of (B i , W i ) in the digitized picture of X ∩ T s .For a set X that is contained in the ball B(0, s), centred at the origin 0 with radius s, the weight functions can be chosen to yield an asymptotically unbiased estimator.This is also possible if the set under consideration is a ball centred at the origin with unknown radius.For general shapes one cannot expect to obtain an unbiased estimator.We propose a method for determining the weight functions, which yields bounds for the expected asymptotic relative worst case error.We illustrate this method in the case of 2 × 2 × 2-configurations.
In the next section basic notations and concepts are introduced together with an asymptotic result (Theorem 2.1) on weighted volumes of morphological transforms.The subsequent section is devoted to the proof of the main theoretical result, formula (1), which is stated with all necessary assumptions in Theorem 3.1.In the penultimate section Theorem 3.1 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, which can be improved whenever X ⊂ B(0, R) for some known radius R > 0; see (15).In the final section we specialize these results to the scaled standard lattice L = tZ 3 , define an estimator based on the m = 102 informative 2 × 2 × 2-configurations, and compare its performance in a simulation example with the theoretical asymptotic results.

By
The positive part of a real valued function f is denoted by f + = max(f, 0).The support function of a convex body K in R d is denoted by h(K, •).We use this notion also for compact sets A, A = ∅, 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 ∈ A c , which do not have a unique nearest point in A. The set exo(A) is measurable and has Lebesgue measure zero, see (Fremlin, 1997).
A closed set X ⊂ R d is gentle if for H d−1 -almost all x ∈ ∂X there are two nondegenerate open balls touching in x such that one of them is contained in X and the other in X c , and if also Here H k is the k-dimensional Hausdorff measure in R d and N(∂X) is the reduced normal bundle of ∂X; for further details see (Kiderlen and Rataj, 2006).The class of gentle sets is rather large.It contains for instance all convex bodies (compact convex subsets of R 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 a unique outer unit normal n(a) to X exists.Let C d−1 (X, •) be the image measure of H d−1 on ∂X under the map a → (a, n(a)).The measure C d−1 (X, •) vanishes outside N(X).Let ξ ∂X : R d \ exo(∂X) → ∂X denote the metric projection.The following theorem is a generalization of (Kiderlen and Rataj, 2006, Theorem 1).
Theorem 2.1.Let X ⊆ R d be a closed gentle set, f : R d → R a compactly supported bounded measurable function and B, W and P, Q four non-empty compact subsets of R d .Then If f is in addition continuous in all points of ∂X, then f (ξ ∂X (x)) can be replaced by f (x) in (2).
Proof.Let C ⊆ R d be a bounded Borel set.Then (2) holds for f = 1 C by (Kiderlen and Rataj, 2006, Theorem 1).It is immediate that (2) also holds for compactly supported measurable step functions.For a non-negative compactly supported bounded measurable function f , let (f k ) k∈N , (g k ) k∈N be sequences of step functions such that (3) Using the monotone convergence theorem we obtain lim sup Note that for applying the monotone convergence theorem to the right-hand side of (3) we can assume that ∪ k∈N supp(g k ) is compact and that the sequence (g k ) k∈N is uniformly bounded.The same argument holds if we take lim inf ε→0+ in (3) and hence the claim is shown for f ≥ 0. For general f = f + − f − we can treat f + and f − separately to obtain (2).If f is in addition continuous in all points of ∂X, we obtain uniform continuity in the following sense.For each η > 0 there exists a δ > 0 such that for all x ∈ ∂X ∩ supp(f ) and This implies the second claim.
Let x 1 , . . ., x d be a basis of R d and let 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 depends only on L and not on the basis chosen.This number is denoted by det(L).
If ξ is a uniform random variable in C 0 , then the random lattice ξ + L is a stationary random lattice.Let X ⊆ R d be a compact gentle set and let the function f be measurable nonnegative or integrable.Let ξ + L be a stationary random lattice, and let B, W ⊆ L be two non-empty finite subsets of L. Define Corollary 2.2.Let X ⊆ R d be a compact gentle set.Let f be a locally bounded measurable function, which is continuous on ∂X and let B and W be two non-empty finite subsets of a lattice L. Then ⊂ C for all t smaller than some fixed t 0 > 0 and X ⊂ C. Replacing f , P , B, Q and W in Theorem 2.1 by 1 C f , {0}, B, {0} and W , respectively, yields the claim.

Combining local and digital stereology
In the following we restrict ourselves to R 3 .The results can be generalized to R d , d ≥ 4, in a straightforward manner.We prefer to present them only in R 3 in order to keep the notation concise.Denote the standard basis vectors in R 3 by e 1 , e 2 , e 3 .Let R be a random proper rotation with distribution given by the normalized Haar measure on SO 3 .Fix the 2-subspace l 0 = span(e 1 , e 2 ).We define the random 2-subspace L = Rl 0 .It is uniformly distributed in the set L of all 2-subspaces of R 3 .Let µ be the distribution of L.
) is called a random 2-slice with thickness 2s.It will be clear from the context whether T s refers to the deterministic 2-slice T s (l) or the random 2-slice T s (L).
Theorem 3.1.Let X ⊆ R 3 be a compact gentle set.Let R be a random proper rotation and let ξ + L be a stationary random lattice, which is independent of R. Let B, W ⊂ L be two non-empty finite subsets of the lattice L and f a continuous non-negative function on R 3 .The weighted sum N f t of occurrences of (B, W ) in the digitization of X, which lie entirely in T s , is given by where H(φ) is given by and φ X a,l is the angle between l and the outer normal of X at a ∈ ∂X.
The proof of Theorem 3.1 is based on the following Proposition 3.2 and Lemmas 3.3, 3.4.The difference of N f t in Theorem 3.1 and Ñg t in Proposition 3.2 is that the latter also counts configurations which do not lie entirely in the slice T 2 .These are of course not observable in practice.
Proposition 3.2.Let X ⊆ R 3 be a compact gentle set.Let R be a random proper rotation and let ξ + L be a stationary random lattice, which is independent of R. Let B, W ⊂ L be two non-empty finite subsets of L. Let g : R 3 × L → R be a nonnegative bounded measurable function such that g(•, l) is continuous for all l ∈ L. Then the number where Proof.For µ-almost all l ∈ L we have that H 2 (∂X ∩ ∂T s ) = 0, which can be seen using Fubini's theorem.This implies that the set X ∩ T s is compact gentle for µalmost all l ∈ L. Applying Corollary 2.2 we obtain that is uniformly bounded for t ≤ 1, hence Lebesgue's dominated convergence theorem yields the assertion.In fact ( 5) implies By assumption there is a constant Hence we obtain Applying (Kiderlen and Rataj, 2006, Proposition 4), which is derived from a farreaching generalization of Steiner's formula (see Hug et al. (2004)), we obtain where µ i (∂X; •) are the support measures of ∂X, and δ(∂X; a, n) = inf{t ≥ 0 | a + tn ∈ exo(∂X)} is the reach function of ∂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 ≤ 1.The same argument can also be applied to Lemma 3.3.Let g : R 3 × L → R be a non-negative measurable function such that g(•, l) is continuous for all l ∈ L. Let F g (r) be given by (7).Then the conditional expectation E[F g (R)|L = l] can be expressed as for µ-almost all l ∈ L, where H(φ) is given by (6).
Proof.For µ-almost all l ∈ L we have that H 2 (∂X ∩ ∂T s ) = 0.This implies that there is a unique normal n(a) for H 2 -almost all a ∈ ∂(X ∩ T s ).Hence we obtain for r ∈ SO 3 There exists a regular version of the conditional distribution of R given L (Klenke, 2006, Satz 8.36).Therefore we can use Fubini's theorem to obtain Recall that µ-almost surely ½ ∂(X∩Ts) = ½ ∂X∩Ts + ½ X∩∂Ts .The claim now follows from Lemma 3.4.
Lemma 3.4.For n ∈ S 2 and l ∈ L we have where φ is the angle between n and l, and H(φ) is given by (6).
Proof.Fix ρ ∈ SO 3 such that ρl = l 0 .Then We have n, l = ρn, l 0 , hence the last conditional expectation in the above equation can be written as the normalized integral over the two small circles ⊂ S 2 parallel to l 0 at height ± sin φ with radius cos φ, where φ is the angle between n and l.
Proof of Theorem 3.1.Let 0 < ε < s/2.For l ∈ L, define the continuous function The right hand side converges pointwise in l to ∂X∩Ts f (a)H(φ X a,l ) H 2 (da), as ε → 0. It is also bounded independently of ε and l, hence by dominated convergence is also given by ( 8).It remains to show that lim Choose q > 0 such that B ∪ W ⊆ B(0, q).Then for t < ε/q we have that where I t (x) = ½ {x+tRB⊆X∩tR(ξ+L),x+tRW ⊆tR(ξ+L)\X} and we used that χ L s−ε,ε (x) = 0 for x ∈ (T s−ε ) C .Using (5) this yields for all r ∈ SO 3 and all t ≤ 1.Then the last expression in the above inequality is bounded by This integral converges for all r ∈ SO 3 by Theorem 2.1 as t → 0+ to where n(a) is the unique normal of X at a.It exists for H 2 -almost all a ∈ ∂X.
The function (−h(rB ⊕ r W , •)) + is bounded by a constant C ′ , independent of r.Therefore, using dominated convergence, the limit of the above integral as ε → 0 is bounded by C ′ H 2 (∂X ∩ ∂T s ) = 0 for almost all r, which yields the claim.

An estimator for surface area
In the following we will derive an estimator for surface area using Theorem 3.1, which is based on a local stereological sampling design.Let X ∈ R 3 be a compact gentle set.Suppose we observe X ∩ R[(l 0 + B(0, s)) ∩ t(ξ + L)] for some random proper rotation R, ξ + L a stationary random lattice, which is independent of R, and t > 0. Let (B i , W i ), i = 1, . . ., m be boundary configurations of L, i.e.B i , W i are non-empty, disjoint finite subsets of L with B i ∪ W i = C 0 ∩ L, where C 0 is the fundamental cell of L. We define the following estimator for surface area where N i t = N λ i t as defined in Theorem 3.1 with B = B i , W = W i .The continuous functions λ i : [0, ∞) → [0, ∞) are weight functions, which have to be suitably chosen according to the choice of (B i , W i ).We give an example for 2 × 2 × 2-configurations in the following section.By Theorem 3.1 we obtain that where H i is given by ( 6) for (B i , W i ).Using Fubini's theorem we can rewrite the right-hand side of (10) as We have L ½ Ts (a)H i (φ X a,l )µ(dl) = g i ( a , ψ X a ), where ψ X a ∈ [0, π] is the angle between a and n, where n is the outer normal of X at a ∈ ∂X.The function g i (r, ψ) for ψ ∈ [0, π] and r ∈ [0, ∞) is given by Let x * denote the two sided cut-off function x → x * = min{1, max{−1, x}}.The function G ψ,q (z) for ψ ∈ (0, π) and q ∈ (0, 1] is given by where α ψ,q (z) = (q − z cos ψ)/(sin ψ √ 1 − z 2 ).For ψ ∈ {0, π} we have Note that for all ψ ∈ [0, π] and all r ∈ (s, ∞) see (Jensen, 1998).We assume that none of the functions H i , i = 1, . . ., m, is identical zero.This is fulfilled, when B i and W i can be strictly separated by a hyperplane for all i = 1, . . ., m.As g i (0, 0) = g i (r, ψ) for all (r, ψ) and as (10) can be rewritten as the choice λ i (r) = a i η i with η i = g i (0, 0) −1 implies that Ŝ(X) is asymptotically unbiased for all X ⊂ B(0, s) whenever the coefficients a 1 , . . ., a m ∈ R are summing up to one.
If the set X under consideration is a ball centred at the origin with unknown radius, we have ψ X a = 0 for all a ∈ ∂X.In this case it is sensible to assume that the sets B i , W i are such that we have g i (r, 0) −1 > 0 for all r ∈ [0, ∞) and i = 1, . . ., m.This is equivalent to requiring that the support of each H i contains an interval [0, ε) for some ε > 0. Choosing λ i (r) = a i g i (r, 0) −1 , where m i=1 a i = 1, yields an asymptotically unbiased estimator Ŝ(X) of S(X).
For general shapes we cannot expect to obtain an unbiased estimator by ( 13).A suitable choice of λ i for r > s will strongly depend on the choice of the pairs (B i , W i ).In the sequel we propose one method to choose the weight functions λ i and show how the asymptotic relative worst case error can be determined in this case.Suppose we can determine coefficients µ i ≥ 0 such that for all z ∈ [0, 1] we have then by ( 12) we obtain for all ψ ∈ [0, π] that m i=1 µ i g i (r, ψ) ≈ f (r) −1 , where f (r) = max{1, r/s}.The function f ( a ) −1 is the probability that a is contained in the random 2-slice T s , see (Jensen, 1998).Setting λ i (r) = µ i f (r), we obtain by ( 13) We suggest to chose (µ 1 , . . ., µ m ) within the set S ⊆ [0, ∞) m of all (µ 1 , . . ., µ m ) such that there exists (a 1 , . . ., a m ) ∈ [0, 1] m with m i=1 a i = 1 and µ i = a i η i .This guarantees that the estimator is asymptotically unbiased for sets X ⊆ B(0, s).
In the remainder of this section we show how to determine the asymptotic relative worst case error for given coefficients (µ 1 , . . ., µ m ) ∈ S. It is immediate from (10) This error bound is independent of the size and shape of X.If we know that X ⊆ B(0, R) for some R > 0, then the worst case error is typically smaller and one can determine a bound using the Lipschitz continuity of the function In order to find a Lipschitz constant with respect to ψ for we use partial integration to rewrite the function for r > s as where Let q = s/r and ψ, q, z such that |α ψ,q (z)| ≤ 1.Then we obtain and hence The above expression is non-negative and bounded by 1 − q 2 .Therefore is bounded by hence we suggest to use the upper bound To find a Lipschitz constant for m i=1 λ i (r)g i (r, ψ) with respect to r we first differentiate with respect to r and obtain The first term of the above expression is bounded by 1/r H(arcsin(•)) ∞ .In order to find a bound for the second term we use partial integration to rewrite m i=1 µ i g i (r, ψ) for r > s as in (14).Then (∂/∂r) m i=1 µ i g i (r, ψ) is given by , and The term on the right-hand side of the above equation is bounded in absolute value by s r 2 π 2 .

Coefficients for
there is a hyperplane, which strictly separates B and W .In this section we want to investigate the surface area estimator (9) in the case, where (B i , W i ) runs through the family of all 102 informative 2 × 2 × 2 configurations.These configurations were thoroughly investigated in Ziegel and Kiderlen (2009).In particular the functions (−h(B i + Wi )) + are explicitly given, hence we can numerically determine the functions H i .As in Ziegel and Kiderlen (2009) we classify the informative 2 × 2 × 2 configurations in 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 affinely dependent in the case of type four, and affinely independent in the case of type five.
For configurations of type one, all functions H i are identical and we denote them by H 1 .The function H 1 (arcsin(z)) for z ∈ [0, 1] is shown in Figure 1.All functions H i (arcsin(z)), z ∈ [−1, 1], are symmetric with respect to the origin, which is why we only display them for values z ∈ [0, 1].For configurations of type two, three and four there are two different functions H i occurring per type.We denote them by H 2,1 , H 2,2 , H 3,1 , H 3,2 and H 4,1 , H 4,2 respectively, see Figures 2, 3 and 4. For configurations of type five all functions H i coincide and we denote them by H 5 , which is displayed in Figure 5. Figure 6 shows all functions H i scaled by the number of their occurrence amongst all functions H i induced by informative 2 × 2 × 2-configurations.The coefficients η i for informative 2 × 2 × 2-configurations are given in Table 1.We have seen that the estimator in (9) with is unbiased for S(X) if X is a subset of B(0, s), whenever (µ 1 , . . ., µ 102 ) ∈ S. We wish to choose (a 1 , . . ., a 102 ) ∈ [0, 1] 102 such that ν m 1 , ν M 1 are small.In (Ziegel and Kiderlen, 2009) a one-parameter family of coefficients (µ ′ 1 (u), . . ., µ ′ 102 (u)), u ∈ [0, 1] was derived that minimizes max for each u ∈ [0, 1].From the definition (6) of the functions H i it is clear that ( 17) is an upper bound for ν m 1 , ν M 1 .Adapting the coefficients (µ ′ 1 (u), . . ., µ ′ 102 (u)) to yield an asymptotically unbiased estimator for spherical shapes, we obtain the family of coefficients (µ 1 (u), . . ., µ 102 (u)), u ∈ [0, 1] given in Table 2.It turns out that (µ 1 (0), . . ., µ 102 (0)) ∈ S or in other words We therefore suggest to set µ i = µ i (0) and λ i as in ( 16).With this choice we obtain ν m 1 = 0.0461, ν M 1 = 0.0225, (d/dz)(H(arcsin(z))) 1 = 5.7002 and L = 13.3768.As a simulation example we consider a cylinder with radius 1 and height 2 centred at 0 which is contained in B(0, √ 2).We observe an isotropic slice of thickness s = 1.For a lattice distance of t = 0.055 we obtain a mean estimated surface area of 18.115 with variance 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 variance 1.205 and mean relative error 1.6%.The asymptotic relative mean error for sets X with X ⊆ B(0, √ 2) and ψ X a ∈ [0, π/4] for all a ∈ X is less than 1.2%.We determined this value numerically using the method described in the previous section.With ε = 0.0030 we obtain ν M 2 = 0.0042 and ν m 2 = 0.0082.In summary this simulation example indicates that the bias of the simple linear surface area estimator (9) with weight functions ( 16) 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.
S d−1 we denote the unit sphere in R d .The standard scalar product on R d is •, • .By a k-subspace we mean a k-dimensional linear subspace of R d .Let A, B ⊂ R d .The reflection of A at the origin is denoted by Ǎ = {−x | x ∈ A}, its complement by A c = R d \A and its topological boundary by ∂A.We write A ⊕ B = {a + b | a ∈ A, b ∈ B} for the Minkowski sum of A and B, and

Figure 1 :Figure 2 :
Figure 1: Plot of the function H i (arcsin(•)) for a configuration of type one.There are 16 configurations of type one (not identifying twins).

Figure 3 :Figure 4 :Figure 5 :Figure 6 :
Figure3: Plot of the functions H i (arcsin(•)) that occur for configurations of type three.There are 32 configurations of type three with H i = H 3,1 (left curve) and 16 configurations with H i = H 3,2 (right curve).